机器算法验证 - 加权最小二乘回归的标准误差 - 吾爱随笔录

加权最小二乘回归的标准误差

机器算法验证回归标准错误稳健标准错误三明治 wls

2022-03-27 11:27:59

对于 OLS，和。我可以“手动”复制这些。 $\hat{\beta} = (X'X)^{-1}X'y$ $\text{var}(\hat{\beta}) = (X'X)^{-1} X' \sigma^2 I X (X'X)^{-1}$

对于 WLS，在对角线中具有异方差误差和权重，和 $W$ $\hat{\beta} = (X'WX)^{-1}X'Wy$

$\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}$ 。我也可以“手动”复制这些。

替换估计来作用于方差估计。例如对于 HC0 (Zeiles 2004 JSS)，使用平方残差。我还可以为 OLS 和 WLS 复制这些“手动”（参见下面的代码）。 $\sigma^2_i$

但是，summary(lm_wls) 产生的 SE 估计值为 0.25081。我无法“手动”复制这个估计。请参阅下面的“### <-- HERE”箭头处的位置。那里的“diag(sig^2_i)”插入了什么？

###########################################################
### 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

1个回答

基本上你已经计算了你需要的一切。缺失的部分只是sig_i应该是残差标准误差除以相应的权重平方根。在 OLS 中，这不是必需的，因为所有权重都是 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)

然后你得到：

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

summary()与and的输出相匹配vcov()：

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

更熟悉的可能是等式，其中完整三明治（= 面包 * 肉 * 面包）的术语已经简化为面包： $\hat \sigma^2 (X^\top W X)^{-1}$

var_betas2a <- resid_var2^2 * solve(t(X) %*% W %*% X)
sqrt(diag(var_betas2a))
##                   z 
## 0.1760843 0.2508150

其它你可能感兴趣的问题

上一篇说 t 分布为正态分布提供“调整”或“估计”是什么意思？下一篇置换也置换依赖吗？