R 中计算异方差稳健标准误差的有效方法
Efficient way to compute Heteroscedasticity Robust standard errors in R
我正在尝试计算 R 中的稳健标准误差。我知道有两种解决方案可以满足我的要求,但速度非常慢。因此,我的问题是是否有一种更有效的方法。例如。已在 Rcpp 中编码的内容。
我的背景是我正在拟合一个具有大量变量(固定效应)的模型。然而,我对这些系数不感兴趣,我只关心推断一个系数(下例中的 X 系数)。
快速解决方案
???
慢解1
library(sandwich)
lmfe<-lm(Y ~ X + factor(strata_ids))
coeftest(lmfe, vcov = vcovHC(lmfe, "HC1"))
慢解2
我从here得到的手动解决方案是:
summaryw <- function(model) {
s <- summary(model)
X <- model.matrix(model)
u2 <- residuals(model)^2
XDX <- 0
## Here one needs to calculate X'DX. But due to the fact that
## D is huge (NxN), it is better to do it with a cycle.
for(i in 1:nrow(X)) {
XDX <- XDX + u2[i]*X[i,]%*%t(X[i,])
}
# inverse(X'X)
XX1 <- solve(t(X)%*%X)
# Variance calculation (Bread x meat x Bread)
varcovar <- XX1 %*% XDX %*% XX1
# degrees of freedom adjustment
dfc <- sqrt(nrow(X))/sqrt(nrow(X)-ncol(X))
# Standard errors of the coefficient estimates are the
# square roots of the diagonal elements
stdh <- dfc*sqrt(diag(varcovar))
t <- model$coefficients/stdh
p <- 2*pnorm(-abs(t))
results <- cbind(model$coefficients, stdh, t, p)
dimnames(results) <- dimnames(s$coefficients)
results
}
这个问题已经有一个很好的答案(即使用lfe::felm()
)。
要获得更快的方法,请尝试新的 fixest 包。使用 OP 示例,
library(fixest)
mod = feols(Y ~ X | strata_ids, data = dat)
## SEs are automatically clustered by the strata_ids FE
mod
## We can compute other SEs on the fly with summary.fixest(), e.g.
summary(mod, se = 'standard') ## vanilla
summary(mod, se = 'white') ## HC
# etc
更一般的教训是避免建模 fixed-effects 作为 R 中的因素......或任何其他语言 TBH。这相当于 DV 方法,并且总是很慢。相反,您需要使用利用 FWL 或其他优化估计方法的 purpose-built 包。
我正在尝试计算 R 中的稳健标准误差。我知道有两种解决方案可以满足我的要求,但速度非常慢。因此,我的问题是是否有一种更有效的方法。例如。已在 Rcpp 中编码的内容。
我的背景是我正在拟合一个具有大量变量(固定效应)的模型。然而,我对这些系数不感兴趣,我只关心推断一个系数(下例中的 X 系数)。
快速解决方案
???
慢解1
library(sandwich)
lmfe<-lm(Y ~ X + factor(strata_ids))
coeftest(lmfe, vcov = vcovHC(lmfe, "HC1"))
慢解2
我从here得到的手动解决方案是:
summaryw <- function(model) {
s <- summary(model)
X <- model.matrix(model)
u2 <- residuals(model)^2
XDX <- 0
## Here one needs to calculate X'DX. But due to the fact that
## D is huge (NxN), it is better to do it with a cycle.
for(i in 1:nrow(X)) {
XDX <- XDX + u2[i]*X[i,]%*%t(X[i,])
}
# inverse(X'X)
XX1 <- solve(t(X)%*%X)
# Variance calculation (Bread x meat x Bread)
varcovar <- XX1 %*% XDX %*% XX1
# degrees of freedom adjustment
dfc <- sqrt(nrow(X))/sqrt(nrow(X)-ncol(X))
# Standard errors of the coefficient estimates are the
# square roots of the diagonal elements
stdh <- dfc*sqrt(diag(varcovar))
t <- model$coefficients/stdh
p <- 2*pnorm(-abs(t))
results <- cbind(model$coefficients, stdh, t, p)
dimnames(results) <- dimnames(s$coefficients)
results
}
这个问题已经有一个很好的答案(即使用lfe::felm()
)。
要获得更快的方法,请尝试新的 fixest 包。使用 OP 示例,
library(fixest)
mod = feols(Y ~ X | strata_ids, data = dat)
## SEs are automatically clustered by the strata_ids FE
mod
## We can compute other SEs on the fly with summary.fixest(), e.g.
summary(mod, se = 'standard') ## vanilla
summary(mod, se = 'white') ## HC
# etc
更一般的教训是避免建模 fixed-effects 作为 R 中的因素......或任何其他语言 TBH。这相当于 DV 方法,并且总是很慢。相反,您需要使用利用 FWL 或其他优化估计方法的 purpose-built 包。