logo标签列表

使用felm输出和标准误差进行预测

rpredictionlfe
9
有没有办法通过使用felm中的投影方法将固定效应排除来获得lfe::felm的预测行为和标准误差？这个问题非常类似于这里的问题，但是那个问题的答案都不能用于估计标准误差或置信度/预测区间。我知道当前没有predict.felm，但我想知道是否有类似于上述链接的解决方法也可以用于估计预测区间。
library(DAAG)
library(lfe)

model1 <- lm(data = cps1, re74 ~ age + nodeg + marr)
predict(model1, newdata = data.frame(age=40, nodeg = 0, marr=1), se.fit = T, interval="prediction")$fit
# Result:        fit      lwr      upr
# 1 18436.18 2339.335 34533.03

model2 <- felm(data = cps1, re74 ~ age | nodeg + marr)
predict(model2, newdata = data.frame(age=40, nodeg = 0, marr=1), se.fit = T, interval="prediction")$fit
# Does not work

目标是估计yhat的预测区间，为此我认为需要计算完整的方差-协方差矩阵（包括固定效应）。我还没有找到如何做到这一点，并且我想知道这是否在计算上可行。

2个回答
4

经过与多人交流后，我认为直接从felm中获取yhat=Xb（其中X包括协变量和固定效应）的分布估计是不可能的，这就是这个问题的核心。但是可以使用引导法来解决。以下代码并行地执行此操作。尽管还有提高性能的空间，但这提供了一般思路。

注意：在此处，我没有计算完整的预测区间，只需要计算Xb的SE即可，但获取预测区间很简单——只需将sigma^2的根加到SE上即可。

library(DAAG)
library(lfe)
library(parallel)

model1 <- lm(data = cps1, re74 ~ age + nodeg + marr)
yhat_lm <- predict(model1, newdata = data.frame(age=40, nodeg = 0, marr=1), se.fit = T)

set.seed(42)
boot_yhat <- function(b) {
  print(b)
  n <- nrow(cps1)
  boot <- cps1[sample(1:n, n, replace=T),]

  lm.model <- lm(data=demeanlist(boot[, c("re74", "age")], list(factor(boot$nodeg), factor(boot$marr))), 
                 formula = re74 ~ age)
  fe <- getfe(felm(data = boot, re74 ~ age | nodeg + marr))

  bootResult <- predict(lm.model, newdata = data.frame(age = 40)) + 
    fe$effect[fe$fe == "nodeg" & fe$idx==0] + 
    fe$effect[fe$fe == "marr" & fe$idx==1]
  return(bootResult)
}

B = 1000
yhats_boot <- mclapply(1:B, boot_yhat)

plot(density(rnorm(10000, mean=yhat_lm$fit, sd=yhat_lm$se.fit)))
lines(density(yhats), col="red")
3

在你的第一个模型中，predict(.)返回以下结果：

#        fit      lwr      upr
# 1 18436.18 2339.335 34533.03

按照李哲源的方法，我们也可以手动实现这些结果。

beta.hat.1 <- coef(model1)  # save coefficients
# model matrix: age=40, nodeg = 0, marr=1:
X.1 <- cbind(1, matrix(c(40, 0, 1), ncol=3))  
pred.1 <- as.numeric(X.1 %*% beta.hat.1) # prediction
V.1 <- vcov(model1)  # save var-cov matrix
se2.1 <- unname(rowSums((X.1 %*% V.1) * X.1))  # prediction var
alpha.1 <- qt((1-0.95)/2, df = model1$df.residual)  # 5 % level
pred.1 + c(alpha.1, -alpha.1) * sqrt(se2.1)  # 95%-CI
# [1] 18258.18 18614.18
sigma2.1 <- sum(model1$residuals ^ 2) / model1$df.residual  # sigma.sq
PI.1 <- pred.1 + c(alpha.1, -alpha.1) * sqrt(se2.1 + sigma2.1) # prediction interval
matrix(c(pred.1, PI.1), nrow = 1, dimnames = list(1, c("fit", "lwr", "upr")))
#        fit      lwr      upr
# 1 18436.18 2339.335 34533.03

现在，将您链接的示例应用于多个前端，我们得到以下结果：
lm.model <- lm(data=demeanlist(cps1[, c(8, 2)], 
                               list(as.factor(cps1$nodeg), 
                                    as.factor(cps1$marr))), re74 ~ age)
fe <- getfe(model2)
predict(lm.model, newdata = data.frame(age = 40)) + fe$effect[fe$idx=="1"]
# [1] 15091.75 10115.21

第一个值是添加了FE的宽度，第二个是未添加FE的宽度（尝试fe$effect[fe$idx=="1"]）。
现在我们正在按照上面的手动方法进行。
beta.hat <- coef(model2)  # coefficient
x <- 40  # age = 40
pred <- as.numeric(x %*% beta.hat)  # prediction
V <- model2$vcv  # var/cov
se2 <- unname(rowSums((x %*% V) * x))  # prediction var
alpha <- qt((1-0.95)/2, df = model2$df.residual)  # 5% level
pred + c(alpha, -alpha) * sqrt(se2)  # CI
# [1]  9599.733 10630.697
sigma2 <- sum(model2$residuals ^ 2) / model2$df.residual  # sigma^2
PI <- pred + c(alpha, -alpha) * sqrt(se2 + sigma2)  # PI
matrix(c(pred, PI), nrow = 1, dimnames = list(1, c("fit", "lwr", "upr")))  # output
#        fit       lwr      upr
# 1 10115.21 -5988.898 26219.33

如我们所见，拟合与链接示例方法相同，但现在具有预测区间。（免责声明：该方法的逻辑应该很简单，PI的值仍应进行评估，例如使用reghdfe在Stata中。）
编辑：如果您想从felm（）实现与线性模型1中predict.lm（）生成的完全相同的输出，则只需在模型中再次“包含”固定效应（请参见下面的model3）。然后按照相同的方法进行即可。为了更方便，您可以将其轻松包装成一个函数。
library(DAAG)
library(lfe)

model3 <- felm(data = cps1, re74 ~ age + nodeg + marr)

pv <- c(40, 0, 1)  # prediction x-values

predict0.felm <- function(mod, pv.=pv) {
  beta.hat <- coef(mod)  # coefficient
  x <- cbind(1, matrix(pv., ncol=3))  # prediction vector
  pred <- as.numeric(x %*% beta.hat)  # prediction
  V <- mod[['vcv'] ] # var/cov
  se2 <- unname(rowSums((x %*% V) * x))  # prediction var
  alpha <- qt((1-0.95)/2, df = mod[['df.residual']])  # 5% level
  CI <- structure(pred + c(alpha, -alpha) * sqrt(se2), 
                  names=c("CI lwr", "CI upr"))  # CI
  sigma2 <- sum(mod[['residuals']] ^ 2) / mod[['df.residual']] # sigma^2
  PI <- pred + c(alpha, -alpha) * sqrt(se2 + sigma2)  # PI
  mx <- matrix(c(pred, PI), nrow = 1, 
               dimnames = list(1, c("PI fit", "PI lwr", "PI upr")))  # output
  list(CI, mx)
}

predict0.felm(model3)[[2]]
#     PI fit   PI lwr   PI upr
# 1 18436.18 2339.335 34533.03

通过使用 felm()，您可以实现与 predict.lm() 相同的预测区间。
1
谢谢你的回答，这有助于我澄清思路，也帮我找到了一些 MWE 中的错误。但我认为这并没有回答我所提出的问题。我想找到一种技术，可以从 felm 得到与使用 lm 和 predict.lm 获得的相同的预测区间。您给出的答案仅使用了 felm 中非固定效应组件的 var/cov 矩阵。我认为这里的根本问题是，在使用 felm 时是否可能估计所有协变量（包括固定效应）的完整 var/cov 矩阵。 - pbaylis
感谢您的评论，我已相应地编辑了我的答案。 - jay.sf
1
不幸的是，在 VCV 中包括固定效应的问题在于我正在运行 felm 而不是 lm，因为在我的实际问题中，固定效应有太多的因素需要在内存中估计 - 这个工作示例是在一个小数据集上进行的，但我想象中的数据集有数百万或数十亿条观测值和成千上万个固定效应。因此，该问题的假设是必须通过 felm 将固定效应清除。 - pbaylis
也许你会发现 lfe::fevcov() 很有用，它可以抛出 FE vcv 矩阵。 - jay.sf
1
我查看了那个函数，但它只包括固定效应的vcov，而我需要完整的方差-协方差矩阵（包括协变量和固定效应之间的协方差）才能得到正确的预测区间。 - pbaylis
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接