This formula fits a linear model, provides a variety of options for robust standard errors, and conducts coefficient tests
lm_robust( formula, data, weights, subset, clusters, fixed_effects, se_type = NULL, ci = TRUE, alpha = 0.05, return_vcov = TRUE, try_cholesky = FALSE )
Arguments
| formula | an object of class formula, as in |
|---|---|
| data | A |
| weights | the bare (unquoted) names of the weights variable in the supplied data. |
| subset | An optional bare (unquoted) expression specifying a subset of observations to be used. |
| clusters | An optional bare (unquoted) name of the variable that corresponds to the clusters in the data. |
| fixed_effects | An optional right-sided formula containing the fixed
effects that will be projected out of the data, such as |
| se_type | The sort of standard error sought. If |
| ci | logical. Whether to compute and return p-values and confidence intervals, TRUE by default. |
| alpha | The significance level, 0.05 by default. |
| return_vcov | logical. Whether to return the variance-covariance matrix for later usage, TRUE by default. |
| try_cholesky | logical. Whether to try using a Cholesky decomposition to solve least squares instead of a QR decomposition, FALSE by default. Using a Cholesky decomposition may result in speed gains, but should only be used if users are sure their model is full-rank (i.e., there is no perfect multi-collinearity) |
Value
An object of class "lm_robust".
The post-estimation commands functions summary and tidy
return results in a data.frame. To get useful data out of the return,
you can use these data frames, you can use the resulting list directly, or
you can use the generic accessor functions coef, vcov,
confint, and predict. Marginal effects and uncertainty about
them can be gotten by passing this object to
margins from the margins,
or to emmeans in the emmeans package.
Users who want to print the results in TeX of HTML can use the
extract function and the texreg package.
If users specify a multivariate linear regression model (multiple outcomes), then some of the below components will be of higher dimension to accommodate the additional models.
An object of class "lm_robust" is a list containing at least the
following components:
the estimated coefficients
the estimated standard errors
the t-statistic
the estimated degrees of freedom
the p-values from a two-sided t-test using coefficients, std.error, and df
the lower bound of the 1 - alpha percent confidence interval
the upper bound of the 1 - alpha percent confidence interval
a character vector of coefficient names
the significance level specified by the user
the standard error type specified by the user
the residual variance
the number of observations used
the number of columns in the design matrix (includes linearly dependent columns!)
the rank of the fitted model
the fitted variance covariance matrix
The \(R^2\), $$R^2 = 1 - Sum(e[i]^2) / Sum((y[i] - y^*)^2),$$ where \(y^*\) is the mean of \(y[i]\) if there is an intercept and zero otherwise, and \(e[i]\) is the ith residual.
The \(R^2\) but penalized for having more parameters, rank
a vector with the value of the F-statistic with the numerator and denominator degrees of freedom
whether or not weights were applied
the original function call
the matrix of predicted means
Details
This function performs linear regression and provides a variety of standard
errors. It takes a formula and data much in the same was as lm
does, and all auxiliary variables, such as clusters and weights, can be
passed either as quoted names of columns, as bare column names, or
as a self-contained vector. Examples of usage can be seen below and in the
Getting Started vignette.
The mathematical notes in
this vignette
specify the exact estimators used by this function.
The default variance estimators have been chosen largely in accordance with the
procedures in
this manual.
The default for the case
without clusters is the HC2 estimator and the default with clusters is the
analogous CR2 estimator. Users can easily replicate Stata standard errors in
the clustered or non-clustered case by setting `se_type` = "stata".
The function estimates the coefficients and standard errors in C++, using
the RcppEigen package. By default, we estimate the coefficients
using Column-Pivoting QR decomposition from the Eigen C++ library, although
users could get faster solutions by setting `try_cholesky` = TRUE to
use a Cholesky decomposition instead. This will likely result in quicker
solutions, but the algorithm does not reliably detect when there are linear
dependencies in the model and may fail silently if they exist.
If `fixed_effects` are specified, both the outcome and design matrix
are centered using the method of alternating projections (Halperin 1962; Gaure 2013). Specifying
fixed effects in this way will result in large speed gains with standard error
estimators that do not need to invert the matrix of fixed effects. This means using
"classical", "HC0", "HC1", "CR0", or "stata" standard errors will be faster than other
standard error estimators. Be wary when specifying fixed effects that may result
in perfect fits for some observations or if there are intersecting groups across
multiple fixed effect variables (e.g. if you specify both "year" and "country" fixed effects
with an unbalanced panel where one year you only have data for one country).
As with `lm()`, multivariate regression (multiple outcomes) will only admit
observations into the estimation that have no missingness on any outcome.
References
Abadie, Alberto, Susan Athey, Guido W Imbens, and Jeffrey Wooldridge. 2017. "A Class of Unbiased Estimators of the Average Treatment Effect in Randomized Experiments." arXiv Pre-Print. https://arxiv.org/abs/1710.02926v2.
Bell, Robert M, and Daniel F McCaffrey. 2002. "Bias Reduction in Standard Errors for Linear Regression with Multi-Stage Samples." Survey Methodology 28 (2): 169-82.
Gaure, Simon. 2013. "OLS with multiple high dimensional category variables." Computational Statistics \& Data Analysis 66: 8-1. doi: 10.1016/j.csda.2013.03.024
Halperin, I. 1962. "The product of projection operators." Acta Scientiarum Mathematicarum (Szeged) 23(1-2): 96-99.
MacKinnon, James, and Halbert White. 1985. "Some Heteroskedasticity-Consistent Covariance Matrix Estimators with Improved Finite Sample Properties." Journal of Econometrics 29 (3): 305-25. doi: 10.1016/0304-4076(85)90158-7 .
Pustejovsky, James E, and Elizabeth Tipton. 2016. "Small Sample Methods for Cluster-Robust Variance Estimation and Hypothesis Testing in Fixed Effects Models." Journal of Business & Economic Statistics. Taylor & Francis. doi: 10.1080/07350015.2016.1247004 .
Samii, Cyrus, and Peter M Aronow. 2012. "On Equivalencies Between Design-Based and Regression-Based Variance Estimators for Randomized Experiments." Statistics and Probability Letters 82 (2). doi: 10.1016/j.spl.2011.10.024 .
Examples
library(fabricatr) dat <- fabricate( N = 40, y = rpois(N, lambda = 4), x = rnorm(N), z = rbinom(N, 1, prob = 0.4) ) # Default variance estimator is HC2 robust standard errors lmro <- lm_robust(y ~ x + z, data = dat) # Can tidy() the data in to a data.frame tidy(lmro)#> term estimate std.error statistic p.value conf.low conf.high #> 1 (Intercept) 3.9737898 0.4485498 8.8591939 1.118544e-10 3.064942 4.8826380 #> 2 x -0.5307197 0.4882497 -1.0869841 2.840745e-01 -1.520008 0.4585682 #> 3 z 0.2195192 0.9326624 0.2353683 8.152208e-01 -1.670234 2.1092726 #> df outcome #> 1 37 y #> 2 37 y #> 3 37 y#> #> Call: #> lm_robust(formula = y ~ x + z, data = dat) #> #> Standard error type: HC2 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF #> (Intercept) 3.9738 0.4485 8.8592 1.119e-10 3.065 4.8826 37 #> x -0.5307 0.4882 -1.0870 2.841e-01 -1.520 0.4586 37 #> z 0.2195 0.9327 0.2354 8.152e-01 -1.670 2.1093 37 #> #> Multiple R-squared: 0.04158 , Adjusted R-squared: -0.01022 #> F-statistic: 0.654 on 2 and 37 DF, p-value: 0.5259# Can also get coefficients three ways lmro$coefficients#> (Intercept) x z #> 3.9737898 -0.5307197 0.2195192#> (Intercept) x z #> 3.9737898 -0.5307197 0.2195192#> [1] 3.9737898 -0.5307197 0.2195192# Can also get confidence intervals from object or with new 1 - `alpha` lmro$conf.low#> (Intercept) x z #> 3.064942 -1.520008 -1.670234#> 10 % 90 % #> (Intercept) 3.3884976 4.5590820 #> x -1.1678144 0.1063751 #> z -0.9974694 1.4365078# Can recover classical standard errors lmclassic <- lm_robust(y ~ x + z, data = dat, se_type = "classical") tidy(lmclassic)#> term estimate std.error statistic p.value conf.low conf.high #> 1 (Intercept) 3.9737898 0.4908831 8.0951850 1.039569e-09 2.979166 4.9684135 #> 2 x -0.5307197 0.4368420 -1.2149007 2.321003e-01 -1.415846 0.3544063 #> 3 z 0.2195192 0.8975742 0.2445694 8.081411e-01 -1.599139 2.0381773 #> df outcome #> 1 37 y #> 2 37 y #> 3 37 y# Can easily match Stata's robust standard errors lmstata <- lm_robust(y ~ x + z, data = dat, se_type = "stata") tidy(lmstata)#> term estimate std.error statistic p.value conf.low conf.high #> 1 (Intercept) 3.9737898 0.4517381 8.7966665 1.338896e-10 3.058481 4.8890982 #> 2 x -0.5307197 0.4798608 -1.1059867 2.758718e-01 -1.503010 0.4415707 #> 3 z 0.2195192 0.9172312 0.2393281 8.121720e-01 -1.638968 2.0780060 #> df outcome #> 1 37 y #> 2 37 y #> 3 37 y# Easy to specify clusters for cluster-robust inference dat$clusterID <- sample(1:10, size = 40, replace = TRUE) lmclust <- lm_robust(y ~ x + z, data = dat, clusters = clusterID) tidy(lmclust)#> term estimate std.error statistic p.value conf.low conf.high #> 1 (Intercept) 3.9737898 0.6045112 6.5735587 0.0004326423 2.520725 5.4268549 #> 2 x -0.5307197 0.3732881 -1.4217430 0.2009822588 -1.426117 0.3646777 #> 3 z 0.2195192 0.7957097 0.2758785 0.7913574413 -1.699550 2.1385881 #> df outcome #> 1 6.479238 y #> 2 6.541055 y #> 3 6.382863 y# Can also match Stata's clustered standard errors lm_robust( y ~ x + z, data = dat, clusters = clusterID, se_type = "stata" )#> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF #> (Intercept) 3.9737898 0.6145043 6.4666589 0.000115865 2.583684 5.3638951 9 #> x -0.5307197 0.3689754 -1.4383606 0.184182152 -1.365400 0.3039607 9 #> z 0.2195192 0.7855876 0.2794331 0.786227898 -1.557603 1.9966418 9# Works just as LM does with functions in the formula dat$blockID <- rep(c("A", "B", "C", "D"), each = 10) lm_robust(y ~ x + z + factor(blockID), data = dat)#> Estimate Std. Error t value Pr(>|t|) CI Lower #> (Intercept) 4.1595712 0.9428286 4.4117999 9.804957e-05 2.243513 #> x -0.5860799 0.5340797 -1.0973641 2.801935e-01 -1.671460 #> z 0.3870498 1.0611313 0.3647521 7.175555e-01 -1.769428 #> factor(blockID)B -0.3962175 1.2706935 -0.3118120 7.570875e-01 -2.978577 #> factor(blockID)C 0.2792504 1.3812823 0.2021675 8.409906e-01 -2.527853 #> factor(blockID)D -0.8210846 1.1856909 -0.6924946 4.933297e-01 -3.230698 #> CI Upper DF #> (Intercept) 6.0756295 34 #> x 0.4993006 34 #> z 2.5435280 34 #> factor(blockID)B 2.1861423 34 #> factor(blockID)C 3.0863538 34 #> factor(blockID)D 1.5885293 34# Weights are also easily specified dat$w <- runif(40) lm_robust( y ~ x + z, data = dat, weights = w, clusters = clusterID )#> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper #> (Intercept) 4.0120339 0.6209751 6.4608607 0.0004885974 2.517552 5.5065158 #> x -0.4938056 0.4225456 -1.1686443 0.2893483806 -1.542660 0.5550485 #> z 0.3181113 1.0831727 0.2936847 0.7816242950 -2.528463 3.1646859 #> DF #> (Intercept) 6.443394 #> x 5.666692 #> z 4.660184# Subsetting works just as in `lm()` lm_robust(y ~ x, data = dat, subset = z == 1)#> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF #> (Intercept) 3.972133 0.7444527 5.335641 0.0003302924 2.313389 5.6308765 10 #> x -2.601287 1.4320438 -1.816486 0.0993453585 -5.792079 0.5895056 10# One can also choose to set the significance level for different CIs lm_robust(y ~ x + z, data = dat, alpha = 0.1)#> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper #> (Intercept) 3.9737898 0.4485498 8.8591939 1.118544e-10 3.064942 4.8826380 #> x -0.5307197 0.4882497 -1.0869841 2.840745e-01 -1.520008 0.4585682 #> z 0.2195192 0.9326624 0.2353683 8.152208e-01 -1.670234 2.1092726 #> DF #> (Intercept) 37 #> x 37 #> z 37# We can also specify fixed effects # Speed gains with fixed effects are greatests with "stata" or "HC1" std.errors tidy(lm_robust(y ~ x + z, data = dat, fixed_effects = ~ blockID, se_type = "HC1"))#> term estimate std.error statistic p.value conf.low conf.high df outcome #> 1 x -0.5860799 0.5156647 -1.1365522 0.2636730 -1.634037 0.4618769 34 y #> 2 z 0.3870498 1.0278841 0.3765501 0.7088478 -1.701862 2.4759618 34 yif (FALSE) { # Can also use 'margins' or 'emmeans' package if you have them installed # to get marginal effects library(margins) lmrout <- lm_robust(y ~ x + z, data = dat) summary(margins(lmrout)) # Can output results using 'texreg' library(texreg) texreg(lmrout) # Using emmeans to obtain covariate-adjusted means library(emmeans) fiber.rlm <- lm_robust(strength ~ diameter + machine, data = fiber) emmeans(fiber.rlm, "machine") }