6
$\begingroup$

For OLS, $\hat{\beta} = (X'X)^{-1}X'y$, and $\text{var}(\hat{\beta}) = (X'X)^{-1} X' \sigma^2 I X (X'X)^{-1}$. I can reproduce these "by hand".

For WLS, with heteroskedastic errors and weights in diagonal $W$, $\hat{\beta} = (X'WX)^{-1}X'Wy$, and

$\text{var}(\hat{\beta}) = (X'WX)^{-1} \hspace{2 pt} X'W \hspace{2 pt} \text{diag}(\sigma^2_i) \hspace{2 pt} WX \hspace{2 pt} (X'WX)^{-1}$. I can also reproduce these "by hand".

Sandwich standard errors act on the variance estimates by substitututing estimates for $\sigma^2_i$. For example for HC0 (Zeiles 2004 JSS) the squared residuals are used. I can also reproduce these "by hand" both for OLS and WLS (see code below).

However, summary(lm_wls) produces an SE estimate of 0.25081. I cannot reproduce this estimate "by hand". See below at the "### <-- HERE" arrow for the spot. What gets plugged in for "diag(sig^2_i)" there?

########################################################### ### DATA GENERATION ####################################### ########################################################### set.seed(1234) # Generate a covariate x <- rnorm(100) # Generate the propensity score ps <- (1 + exp(-(-.5 + .8*x)))^-1 # Generate the exposure (i.e., treatment) variable z <- rbinom(n = 100, size = 1, prob = ps) # Generate the outcome y <- 1.1*x - .6*z + rnorm(100, sd = .5) ### Estimate the average treatment effect by OLS ########################################################### ### ORDINARY LEAST SQUARES REGRESSION ##################### ########################################################### ### OLS via lm() lm_ols <- lm(formula = y ~ x + z) summary(lm_ols) ### Verify OLS 'by hand' X <- cbind(1, x, z) (betas <- solve(t(X) %*% X) %*% t(X) %*% y) (resid_var <- sqrt(sum(lm_ols$residuals^2)/(100 - 3))) (var_betas <- solve(t(X) %*% X) %*% (t(X) %*% diag(resid_var^2, 100, 100) %*% X) %*% solve(t(X) %*% X)) sqrt(diag(var_betas)) # SEs -- Compare to summary(lm_ols) ### Sandwich SEs based on squared residuals ### (see p.4 of Zeiles 2004 JSS) library(sandwich) (vcovHC(x = lm_ols, type = "HC0", sandwich = TRUE)) sqrt(diag(vcovHC(x = lm_ols, type = "HC0", sandwich = TRUE))) ### Verify Sandwich SEs for the OLS model var_betas <- solve(t(X) %*% X) %*% t(X) %*% diag(lm_ols$residuals^2) %*% X %*% solve(t(X) %*% X) sqrt(diag(var_betas)) # SEs -- Compare to sandwich SEs ### Estimate the average treatment effect by propensity ### score weighting ########################################################### ### INVERSE PROPENSITY SCORE WEIGHTING #################### ########################################################### ### Estimate propensity scores glm1 <- glm(z ~ x, family = "binomial") ps_est <- predict(glm1, type = "response") ### Create inverse ps weights wts <- (z/ps_est) + (1-z)/(1-ps_est) ### Estimate average treatment effect via WLS lm_wls <- lm(formula = y ~ z, weights = wts) summary(lm_wls) ### Verify WLS X <- cbind(1, z) W <- diag(wts) (betas2 <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% y) (resid_var2 <- sqrt(sum(wts*(lm_wls$residuals^2))/(100 - 2))) (var_betas2 <- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% "diag(sig^2_i)" %*% t(W) %*% X) %*% ### <-- HERE solve(t(X) %*% W %*% X)) sqrt(diag(var_betas2)) # SEs ### Sandwich SEs with sandwich package library(sandwich) vcovHC(x = lm_wls, type = "HC0", sandwich = TRUE) sqrt(diag(vcovHC(x = lm_wls, type = "HC0", sandwich = TRUE))) ### Verify Sandwich SEs for the WLS model (var_betas3 <- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% diag(lm_wls$resid^2) %*% t(W) %*% X) %*% solve(t(X) %*% W %*% X)) sqrt(diag(var_betas3)) # SEs -- Compare to sandwich SEs
Improve this question
$\endgroup$

1 Answer 1

Reset to default
7
$\begingroup$

Essentially you already computed everything you need. The missing piece is just that the sig_i should be the residual standard error divided by the corresponding square root of the weight. In OLS this isn't necessary because all weights are 1.

sig_i <- resid_var2 / sqrt(wts) var_betas2 <- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% diag(sig_i^2) %*% t(W) %*% X) %*% solve(t(X) %*% W %*% X)

And then you get:

sqrt(diag(var_betas2)) ## z ## 0.1760843 0.2508150

which matches the output of summary() and vcov():

sqrt(diag(vcov(lm_wls))) ## (Intercept) z ## 0.1760843 0.2508150

Even more familiar might be the equations as $\hat \sigma^2 (X^\top W X)^{-1}$ where the terms from the full sandwich (= bread * meat * bread) have already been simplified to just the bread:

var_betas2a <- resid_var2^2 * solve(t(X) %*% W %*% X) sqrt(diag(var_betas2a)) ## z ## 0.1760843 0.2508150
Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ This is exactly what I was looking for. Can you recommend a regression textbook that works out the matrix algebra for WLS? And thanks for your work on sandwich and other useful packages. $\endgroup$ Commented Jun 13, 2017 at 20:38
  • 2
    $\begingroup$ @bsbk Thanks, glad if the packages are useful. As for WLS: This should be covered in all standard econometrics textbooks, e.g., the 2003 edition of Greene's Econometric Analysis covers it in Chapter 11.5. $\endgroup$ Commented Jun 13, 2017 at 21:00

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.