# Standard Errors with Weighted Least Squares Regression

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


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

• 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. – bsbk Jun 13 '17 at 20:38
• @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. – Achim Zeileis Jun 13 '17 at 21:00