Package {estimatr}

estimatr: Fast Estimators for Design-Based Inference

Description

Fast procedures for a small set of commonly-used, design-appropriate estimators with robust standard errors and confidence intervals. Provides lm_robust, lm_lin, iv_robust, lh_robust, difference_in_means, and horvitz_thompson, with feols-style fixed effects absorption and design-aware Horvitz-Thompson variance. See vignette("estimatr2.0") for what changes in version 2.0 and what does not.

Author(s)

Maintainer: Alexander Coppock acoppock@gmail.com

Authors:

• Alexander Coppock acoppock@gmail.com

• Graeme Blair graeme.blair@gmail.com

• Jasper Cooper jjc2247@columbia.edu

• Luke Sonnet luke.sonnet@gmail.com

Other contributors:

• Macartan Humphreys macartan@gmail.com [contributor]

• Neal Fultz nfultz@gmail.com [contributor]

• Lily Medina lilymiru@gmail.com [contributor]

• Russell Lenth russell-lenth@uiowa.edu [contributor]

• Molly Offer-Westort mollyow@uchicago.edu [contributor]

See Also

Useful links:

• https://declaredesign.org/r/estimatr/

• https://github.com/DeclareDesign/estimatr

• Report bugs at https://github.com/DeclareDesign/estimatr/issues

Augment a Model Object with Fitted Values and Residuals

Description

Returns the model frame with .fitted and .resid columns appended, the form downstream packages expect from broom::augment(). Supplying newdata returns that instead, with .fitted only.

Usage

## S3 method for class 'lm_robust'
augment(x, data = NULL, newdata = NULL, ...)

## S3 method for class 'iv_robust'
augment(x, data = NULL, newdata = NULL, ...)

Arguments

Value

A data.frame.

Examples

set.seed(55)
dat <- data.frame(x = rnorm(50), z = rep(0:1, 25))
dat$y <- dat$x + 0.4 * dat$z + rnorm(50)
fit <- lm_robust(y ~ x + z, data = dat)

head(augment(fit))

# Supplying newdata returns predictions on it, with .fitted only
head(augment(fit, newdata = dat[1:5, ]))

Design-Based Difference-in-Means Estimator

Description

Estimates an average treatment effect as a difference in means, choosing the point estimator, the variance, and the degrees of freedom that match the randomization rather than requiring you to name them. Unit, cluster, blocked, block-cluster, matched-pair, and matched-pair cluster designs are recognised, and the design element of the result reports which case applied.

Usage

difference_in_means(
  formula,
  data,
  blocks,
  clusters,
  weights,
  subset,
  se_type = c("default", "none"),
  condition1 = NULL,
  condition2 = NULL,
  ci = TRUE,
  alpha = 0.05
)

Arguments

Details

Selects the appropriate point estimate, standard errors, and degrees of freedom for unit randomized, cluster randomized, block randomized, block-cluster randomized, matched-pairs, and matched-pair cluster randomized designs.

Blocks of different sizes. For unit randomized blocks, blocks are classified by how many units each arm holds rather than by how large the block is. A block with at least two treated and two control units has its own Neyman variance. A block with a single treated or single control unit has no estimable within-block variance, so the variation across such blocks stands in for it. A design containing both kinds combines the two parts by squared share of the sample, following Pashley and Miratrix (2021). The design element of the returned object reports which case applied: "Blocked", "Matched-pair", "Small blocks", or "Hybrid blocked".

Two designs are refused, because the variance genuinely cannot be estimated: exactly one block with a singleton arm, and a set of different-sized such blocks in which one holds half or more of their units. Both messages suggest merging blocks or using lm_robust() with block fixed effects.

If weights are specified, estimation is handed to lm_robust() with HC2 standard errors.

Blocks of clusters. Pashley and Miratrix treat treatment assigned within blocks, not blocks of clusters, so blocked designs with clusters use the earlier estimators. Every block must have at least two treated and two control clusters, unless the design is matched-pair clustered, where the variance is estimated across blocks. A block with a single treated or control cluster is refused: its within-block variance is not estimable, and estimating it anyway understates the standard error by roughly the block's cluster count.

Value

An object of class "difference_in_means", a list holding coefficients, std.error, df, statistic, p.value, conf.low, conf.high, term, outcome, condition1, condition2, vcov, nobs, alpha, and design, a string naming the case that applied: "Standard", "Blocked", "Small blocks", "Hybrid blocked", "Matched-pair", "Clustered", "Block-clustered", or "Matched-pair clustered".

References

Gerber, Alan S. and Donald P. Green. 2012. Field Experiments: Design, Analysis, and Interpretation. New York: W.W. Norton.

Imai, Kosuke, Gary King, and Clayton Nall. 2009. "The Essential Role of Pair Matching in Cluster-Randomized Experiments." Statistical Science 24(1): 29-53. doi:10.1214/08-STS274.

Pashley, Nicole E. and Luke W. Miratrix. 2021. "Insights on Variance Estimation for Blocked and Matched Pairs Designs." Journal of Educational and Behavioral Statistics 46(3): 271-296. doi:10.3102/1076998620946272.

Examples

set.seed(30)
dat <- data.frame(y = rnorm(100), z = rep(0:1, 50))

# Unblocked, unclustered: the Welch-corrected two-sample difference
fit <- difference_in_means(y ~ z, data = dat)
fit
fit$design

# Blocked designs use the Neyman variance within each block
dat_bl <- data.frame(
  bl = rep(1:10, each = 10),
  z  = rep(rep(0:1, each = 5), times = 10)
)
dat_bl$y <- rnorm(100) + 0.3 * dat_bl$z
difference_in_means(y ~ z, data = dat_bl, blocks = bl)

# Matched pairs are recognised as such
dat_pr <- data.frame(pr = rep(1:50, each = 2), z = rep(c(0, 1), 50))
dat_pr$y <- rnorm(100) + 0.3 * dat_pr$z
difference_in_means(y ~ z, data = dat_pr, blocks = pr)$design

# Blocks of unequal shape, which earlier versions refused, use the
# Pashley and Miratrix (2021) estimators. `design` reports which case
# applied rather than leaving it to be inferred from the block sizes.
dat_hy <- rbind(dat_bl[c("bl", "z", "y")],
                transform(dat_pr[c("pr", "z", "y")], bl = pr + 100)[c("bl", "z", "y")])
difference_in_means(y ~ z, data = dat_hy, blocks = bl)$design

# Clustered assignment
dat_cl <- data.frame(cl = rep(1:20, each = 5))
dat_cl$z <- rep(rep(0:1, each = 5), times = 10)
dat_cl$y <- rnorm(100) + 0.3 * dat_cl$z
difference_in_means(y ~ z, data = dat_cl, clusters = cl)

Functions removed in 2.0

Description

commarobust() and starprep() were helpers for producing robust standard errors outside the package's own estimators and formatting them for stargazer. Both are removed.

Usage

commarobust(...)

starprep(...)

Arguments

Details

They are kept here as names that error rather than deleted outright, so that a script written against estimatr 1.x says what happened and what to do instead of failing with ⁠could not find function⁠.

commarobust() recomputed robust standard errors on a fitted lm. Fit the model with lm_robust() instead, which is what it was reimplementing.

starprep() prepared a list of fits for stargazer, which has not been maintained for years. Table-building now goes through modelsummary, which reads tidy() and glance() and therefore works on every estimator in this package without any adapter.

Value

Never returns; both functions signal an error.

Examples

# Both of these error. The replacements:
set.seed(1)
dat <- data.frame(y = rnorm(20), z = rep(0:1, 10))

# was: commarobust(lm(y ~ z, data = dat))
lm_robust(y ~ z, data = dat)

# was: starprep(fit1, fit2) |> stargazer::stargazer()
# now: modelsummary::modelsummary(list(fit1, fit2))

Tidy an estimatr object

Description

Tidy an estimatr object

Usage

## S3 method for class 'lm_robust'
tidy(x, conf.int = TRUE, conf.level = NULL, ...)

## S3 method for class 'iv_robust'
tidy(x, conf.int = TRUE, conf.level = NULL, ...)

## S3 method for class 'difference_in_means'
tidy(x, conf.int = TRUE, conf.level = NULL, ...)

## S3 method for class 'horvitz_thompson'
tidy(x, conf.int = TRUE, conf.level = NULL, ...)

## S3 method for class 'lh_robust'
tidy(x, conf.int = TRUE, conf.level = NULL, ...)

## S3 method for class 'lh'
tidy(x, conf.int = TRUE, conf.level = NULL, ...)

Arguments

Value

A tibble with one row per term (and per outcome, for a multivariate fit): term, estimate, std.error, statistic, p.value, conf.low, conf.high, df, and outcome. A tibble rather than a plain data frame, as broom's tidiers return; 1.x returned a data frame.

Examples

set.seed(50)
dat <- data.frame(x = rnorm(50), z = rep(0:1, 25))
dat$y <- dat$x + 0.4 * dat$z + rnorm(50)
fit <- lm_robust(y ~ x + z, data = dat)

# One row per term, with the interval the fit was built with
tidy(fit)
tidy(fit, conf.int = FALSE)
tidy(fit, conf.level = 0.9)

# The same shape for every estimator in the package
tidy(difference_in_means(y ~ z, data = dat))

Extract model data for the texreg package

Description

Prepares an lm_robust or iv_robust fit for texreg. Largely a clone of texreg's own extract.lm method.

Usage

extract.lm_robust(
  model,
  include.ci = TRUE,
  include.rsquared = TRUE,
  include.adjrs = TRUE,
  include.nobs = TRUE,
  include.fstatistic = FALSE,
  include.rmse = TRUE,
  include.nclusts = TRUE,
  ...
)

extract.iv_robust(
  model,
  include.ci = TRUE,
  include.rsquared = TRUE,
  include.adjrs = TRUE,
  include.nobs = TRUE,
  include.fstatistic = FALSE,
  include.rmse = TRUE,
  include.nclusts = TRUE,
  ...
)

Arguments

Details

These are exported as plain functions rather than registered with S3method() because that is how texreg finds them: it looks up ⁠extract.<class>⁠ by name in the package namespace rather than dispatching on a generic it owns. Registering them the usual way would leave texreg unable to see them.

texreg is the only consumer. Table building through modelsummary needs nothing here, since it reads tidy() and glance() and so already works on every estimator in this package.

Value

A texreg object.

Examples

set.seed(60)
dat <- data.frame(x = rnorm(50), z = rep(0:1, 25))
dat$y <- dat$x + 0.4 * dat$z + rnorm(50)
fit <- lm_robust(y ~ x + z, data = dat)

if (requireNamespace("texreg", quietly = TRUE)) {
  texreg::screenreg(fit)
}

Horvitz-Thompson Estimator with Inverse Probability Weighting

Description

Estimates treatment effects via inverse probability weighting when treatment assignment probabilities are known. Supports all randomizr designs as well as arbitrary designs supplied via a permutation matrix.

Usage

horvitz_thompson(
  formula,
  data,
  condition_prs = NULL,
  condition1 = NULL,
  condition2 = NULL,
  se_type = "youngs",
  ci = TRUE,
  alpha = 0.05
)

Arguments

Details

With more than two arms, condition1 and condition2 select the contrast, and the estimand remains the average treatment effect over all N units the design covers. The estimator therefore divides by N, not by the number of units landing in the two conditions, and the variance uses the joint assignment probabilities implied by the arm sizes. data must hold one row per unit of the design, in the design's order, including units assigned to arms outside the contrast.

Value

An object of class "horvitz_thompson" with fields coefficients, std.error, statistic, p.value, conf.low, conf.high, df, nobs (the number of units in the design, including any arms outside the contrast), vcov, se_type, condition1, condition2, outcome, and term.

Examples

set.seed(40)
dat <- data.frame(y = rnorm(100), z = rep(0:1, 50))

# A named vector of condition probabilities gives the conservative
# simple-randomization bound, valid for any design
horvitz_thompson(y ~ z, data = dat, condition_prs = c("0" = 0.5, "1" = 0.5))

# Passing the randomization declaration instead is what buys the
# design-aware variance, and it is the recommended form
if (requireNamespace("randomizr", quietly = TRUE)) {
  decl <- randomizr::declare_ra(N = 100, m = 50)
  dat$z2 <- randomizr::conduct_ra(decl)
  print(horvitz_thompson(y ~ z2, data = dat, condition_prs = decl))

  # Blocked and clustered designs need no extra arguments: the declaration
  # already carries the structure
  bl <- rep(1:4, each = 25)
  decl_bl <- randomizr::declare_ra(blocks = bl, prob = 0.4)
  dat$z3 <- randomizr::conduct_ra(decl_bl)
  print(horvitz_thompson(y ~ z3, data = dat, condition_prs = decl_bl))

  # Any two arms of a multi-arm design can be contrasted, with the estimand
  # still defined over all N units
  decl3 <- randomizr::declare_ra(N = 100, conditions = c("a", "b", "c"))
  dat$z4 <- randomizr::conduct_ra(decl3)
  print(horvitz_thompson(y ~ z4, data = dat, condition_prs = decl3,
                         condition1 = "a", condition2 = "c"))
}

Two-Stage Least Squares Instrumental Variables Regression

Description

Fits a two-stage least squares instrumental variables regression and returns heteroskedasticity-robust or cluster-robust standard errors, with optional weak-instrument, Wu-Hausman, and overidentification diagnostics.

Usage

iv_robust(
  formula,
  data,
  weights,
  subset,
  clusters,
  fixed_effects,
  se_type = NULL,
  ci = TRUE,
  alpha = 0.05,
  diagnostics = FALSE,
  return_vcov = TRUE,
  try_cholesky = FALSE
)

Arguments

Value

An object of class "iv_robust", a list holding the estimate table in coefficients, std.error, df, statistic, p.value, conf.low, conf.high, term, and outcome; the fit in fitted.values, residuals, vcov, nobs, k, rank, df.residual, and res_var; the summary statistics r.squared, adj.r.squared, tss, and fstatistic; and se_type, weighted, clustered, fes, alpha, terms, xlevels, and call.

residuals are the structural residuals, ⁠y - X beta⁠, rather than the second-stage ones. ei.iv, terms_regressors, and formula record the two-stage structure. With diagnostics = TRUE the object also holds diagnostic_first_stage_fstatistic, diagnostic_endogeneity_test, and diagnostic_overid_test.

Examples

set.seed(25)
n <- 200
dat <- data.frame(z = rbinom(n, 1, 0.5), cl = rep(1:20, each = 10))
dat$x <- dat$z * rbinom(n, 1, 0.7)
dat$y <- dat$x + rnorm(n)

# Endogenous regressor on the left of the bar, instrument on the right
fit <- iv_robust(y ~ x | z, data = dat)
tidy(fit)

# The same variance menu as lm_robust()
iv_robust(y ~ x | z, data = dat, se_type = "classical")
iv_robust(y ~ x | z, data = dat, clusters = cl)

# Weak-instrument, endogeneity, and overidentification tests
summary(iv_robust(y ~ x | z, data = dat, diagnostics = TRUE))

Linear Hypothesis Test for OLS with Robust Standard Errors

Description

Tests a linear combination of coefficients, or several of them jointly, from a model fitted by lm_robust(). The robust variance and the degrees of freedom of the fit are carried through, so a clustered fit is tested on its cluster-adjusted degrees of freedom rather than on the residual ones.

Usage

lh_robust(..., data, linear_hypothesis)

Arguments

Value

An object of class "lh_robust" with three components: lm_robust, the underlying fit; lh, one row per hypothesis holding coefficients, std.error, statistic, p.value, alpha, conf.low, conf.high, df, term, and outcome; and joint_hypothesis, the Wald F test of all of them at once, as value, numdf, dendf, and p.value. Under se_type = "CR2" each hypothesis's df is its own Satterthwaite approximation, as clubSandwich::linear_contrast() computes it, and dendf is the smallest of them.

Examples

set.seed(35)
dat <- data.frame(x = rnorm(100), z = rbinom(100, 1, 0.5),
                  cl = rep(1:10, each = 10))
dat$y <- dat$x + 0.5 * dat$z + rnorm(100)

# One linear combination of coefficients
fit <- lh_robust(y ~ x + z, data = dat, linear_hypothesis = "z + 2*x = 0")
fit
tidy(fit)

# Degrees of freedom follow the fit, so a clustered model tests against the
# cluster-adjusted df rather than the residual df
lh_robust(y ~ x + z, data = dat, clusters = cl,
          linear_hypothesis = "z + 2*x = 0")

# Several restrictions at once give one joint Wald test as well
joint <- lh_robust(y ~ x + z, data = dat, linear_hypothesis = c("x = 0", "z = 0"))
joint$joint_hypothesis

Linear Regression with Lin (2013) Covariate Adjustment

Description

Estimates an average treatment effect with covariate adjustment following Lin (2013): every covariate is centered, interacted with treatment, and entered alongside it. Centering is what makes the treatment coefficient the effect estimate, and the interactions avoid the bias Freedman (2008) identified in ordinary covariate-adjusted regression.

Usage

lm_lin(
  formula,
  covariates,
  data,
  weights,
  subset,
  clusters,
  se_type = NULL,
  ci = TRUE,
  alpha = 0.05,
  return_vcov = TRUE,
  try_cholesky = FALSE
)

Arguments

Value

An object of class "lm_robust", as returned by lm_robust(), with two additions: scaled_center, the covariate means used for centering (taken after any function in the formula is evaluated), and treatment_levels. The treatment row of coefficients is the estimate of the average treatment effect.

References

Lin, Winston. 2013. "Agnostic Notes on Regression Adjustments to Experimental Data: Reexamining Freedman's Critique." The Annals of Applied Statistics 7(1): 295-318. doi:10.1214/12-AOAS583.

Examples

set.seed(20)
dat <- data.frame(
  x  = rnorm(40, mean = 2.3),
  x2 = rpois(40, lambda = 2),
  x3 = runif(40),
  z  = rep(0:1, 20),
  cl = rep(1:20, each = 2)
)
dat$y <- rnorm(40) + dat$x + 0.35 * dat$z

# lm_robust's interface plus one argument
fit <- lm_lin(y ~ z, covariates = ~ x, data = dat)
tidy(fit)

# Several covariates
lm_lin(y ~ z, covariates = ~ x + x2, data = dat)

# Covariates are centered after any function in the formula is evaluated
fit2 <- lm_lin(y ~ z, covariates = ~ x + log(x3), data = dat)
fit2$scaled_center["log(x3)"]
mean(log(dat$x3))

# Clusters, and multi-valued treatments whether or not they are factors
lm_lin(y ~ z, covariates = ~ x, data = dat, clusters = cl)
dat$z3 <- rep(1:3, length.out = 40)
lm_lin(y ~ z3, covariates = ~ x, data = dat)
lm_lin(y ~ factor(z3), covariates = ~ x, data = dat)

# Dropping the intercept gives the mean outcome under each condition
lm_lin(y ~ z3 - 1, covariates = ~ x, data = dat)

Ordinary Least Squares with Robust Standard Errors

Description

Fits a linear model by ordinary least squares and returns heteroskedasticity-robust or cluster-robust standard errors, with the small-sample corrections used in design-based work. Fixed effects can be absorbed rather than expanded into dummy columns, at no cost in the available standard error types.

Usage

lm_robust(
  formula,
  data,
  weights,
  subset,
  clusters,
  fixed_effects,
  se_type = NULL,
  ci = TRUE,
  alpha = 0.05,
  return_vcov = TRUE,
  try_cholesky = FALSE
)

Arguments

Value

An object of class "lm_robust", a list holding the estimate table in coefficients, std.error, df, statistic, p.value, conf.low, conf.high, term, and outcome; the fit in fitted.values, residuals, vcov, nobs, k, rank, df.residual, and res_var; the summary statistics r.squared, adj.r.squared, tss, and fstatistic; and se_type, weighted, clustered, fes, alpha, terms, xlevels, and call.

Absorbed fits add fixed_effects, felevels (the absorbed levels of each factor), and the within-projection summaries proj_r.squared, proj_adj.r.squared, proj_tss, and proj_fstatistic.

Examples

set.seed(15)
dat <- data.frame(
  y = rpois(40, lambda = 4),
  x = rnorm(40),
  z = rbinom(40, 1, prob = 0.4),
  cl = rep(1:10, each = 4),
  bl = rep(c("A", "B", "C", "D"), each = 10),
  w = runif(40)
)

# HC2 is the default
fit <- lm_robust(y ~ x + z, data = dat)
fit
tidy(fit)
summary(fit)
confint(fit, level = 0.8)

# Other variance estimators, including Stata's
lm_robust(y ~ x + z, data = dat, se_type = "classical")
lm_robust(y ~ x + z, data = dat, se_type = "stata")

# Clustered inference defaults to CR2
lm_robust(y ~ x + z, data = dat, clusters = cl)
lm_robust(y ~ x + z, data = dat, clusters = cl, se_type = "stata")

# Weights and subsets behave as they do in lm()
lm_robust(y ~ x + z, data = dat, weights = w, clusters = cl)
lm_robust(y ~ x, data = dat, subset = z == 1)

# Fixed effects are absorbed rather than expanded into dummies. With a
# single factor the HC2 default is exact and costs nothing extra.
lm_robust(y ~ z, data = dat, fixed_effects = ~ bl)

Internal method that creates linear fits

Description

Internal method that creates linear fits

Usage

lm_robust_fit(
  y,
  X,
  weights,
  cluster,
  ci = TRUE,
  se_type,
  has_int,
  alpha = 0.05,
  return_vcov = TRUE,
  return_fit = TRUE,
  try_cholesky = FALSE,
  iv_stage = list(0),
  fe_rank = 0L,
  fe_leverage = NULL,
  femat = NULL,
  linear_hypothesis = NULL
)

Arguments

Examples

# The fitter behind lm_robust(), exported for packages that have already
# built their own design matrix. Most users want lm_robust().
set.seed(45)
X <- cbind(`(Intercept)` = 1, x = rnorm(50))
y <- X[, "x"] + rnorm(50)

lm_robust_fit(
  y = y, X = X,
  weights = NULL, cluster = NULL,
  se_type = "HC2", has_int = TRUE
)

Predict method for lm_robust object

Description

Produces predicted values, obtained by evaluating the regression function in the frame newdata for fits from lm_robust() and lm_lin(). If se.fit is TRUE, standard errors of the predictions are calculated. Setting interval adds confidence or prediction (tolerance) intervals at the level set by alpha, sometimes called narrow and wide intervals respectively.

Usage

## S3 method for class 'lm_robust'
predict(
  object,
  newdata,
  se.fit = FALSE,
  interval = c("none", "confidence", "prediction"),
  alpha = 0.05,
  na.action = na.pass,
  pred.var = NULL,
  weights,
  ...
)

Arguments

Details

Called without newdata, the method returns the in-sample fitted values, and neither se.fit nor interval is available.

The equation used for the standard error of a prediction given a row of data x is:

\sqrt(x \Sigma x'),

where \Sigma is the estimated variance-covariance matrix from lm_robust().

The prediction intervals are for a single observation at each case in newdata with error variance(s) pred.var. The default is to assume that future observations have the same error variance as those used for fitting, which is taken from the fitted lm_robust() object. If weights is supplied, the inverse of those weights scales the variance. If the fit was weighted, the default is to assume constant prediction variance, with a warning.

Value

A numeric vector of predictions, or a data frame with the predictions and their standard errors and interval bounds when se.fit or interval is set.

See Also

lm_robust(), lm_lin()

Examples

# Set seed
set.seed(42)

# Simulate data
n <- 10
dat <- data.frame(y = rnorm(n), x = rnorm(n))

# Fit lm
lm_out <- lm_robust(y ~ x, data = dat)
# In-sample fitted values
predict(lm_out)
# Get predicted fits
fits <- predict(lm_out, newdata = dat)
# With standard errors and confidence intervals
fits <- predict(lm_out, newdata = dat, se.fit = TRUE, interval = "confidence")

# Use new data as well
new_dat <- data.frame(x = runif(n, 5, 8))
predict(lm_out, newdata = new_dat)

# You can also supply custom variance weights for prediction intervals
new_dat$w <- runif(n)
predict(lm_out, newdata = new_dat, weights = w, interval = "prediction")

# Works for 'lm_lin' models as well
dat$z <- sample(1:3, size = nrow(dat), replace = TRUE)
lmlin_out1 <- lm_lin(y ~ z, covariates = ~ x, data = dat)
predict(lmlin_out1, newdata = dat, interval = "prediction")

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

generics

augment(), glance(), tidy()