从lfe包中预测from的方法

Question

有没有人有一个好的、干净的方法来获取felm模型的predict行为？

library(lfe)
model1 <- lm(data = iris, Sepal.Length ~ Sepal.Width + Species)
predict(model1, newdata = data.frame(Sepal.Width = 3, Species = "virginica"))
# Works

model2 <- felm(data = iris, Sepal.Length ~ Sepal.Width | Species)
predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica"))
# Does not work

Answer 1

更新(2020-04-02)：下面来自Grant的answer使用新的包fixest提供了一个更简洁的解决方案。

作为一种解决办法，您可以按如下方式组合felm、getfe和demeanlist：

library(lfe)

lm.model <- lm(data=demeanlist(iris[, 1:2], list(iris$Species)), Sepal.Length ~ Sepal.Width)
fe <- getfe(felm(data = iris, Sepal.Length ~ Sepal.Width | Species))
predict(lm.model, newdata = data.frame(Sepal.Width = 3)) + fe$effect[fe$idx=="virginica"]

其思想是使用demeanlist将变量居中，然后使用lm使用居中的变量估计Sepal.Width上的系数，从而得到一个lm对象，您可以在该对象上运行predict。然后运行felm+getfe以获得固定效果的条件平均值，并将其添加到predict的输出中。

Answer 2

虽然晚了一步，但是新的fixest package (link)有一个预测方法。支持高维固定效果(以及聚类等)使用与lfe非常相似的语法。值得注意的是，对于我测试过的基准测试用例，它也比lfe快得多。

library(fixest)

model_feols <- feols(data = iris, Sepal.Length ~ Sepal.Width | Species)
predict(model_feols, newdata = data.frame(Sepal.Width = 3, Species = "virginica"))
# Works

Answer 3

这应该适用于以下情况:您希望忽略预测中的组效应，正在预测新的X，并且只想要置信区间。它首先查找clustervcv属性，然后查找robustvcv，最后查找vcv。

predict.felm <- function(object, newdata, se.fit = FALSE,
                         interval = "none",
                         level = 0.95){
  if(missing(newdata)){
    stop("predict.felm requires newdata and predicts for all group effects = 0.")
  }

  tt <- terms(object)
  Terms <- delete.response(tt)
  attr(Terms, "intercept") <- 0

  m.mat <- model.matrix(Terms, data = newdata)
  m.coef <- as.numeric(object$coef)
  fit <- as.vector(m.mat %*% object$coef)
  fit <- data.frame(fit = fit)

  if(se.fit | interval != "none"){
    if(!is.null(object$clustervcv)){
      vcov_mat <- object$clustervcv
    } else if (!is.null(object$robustvcv)) {
      vcov_mat <- object$robustvcv
    } else if (!is.null(object$vcv)){
      vcov_mat <- object$vcv
    } else {
      stop("No vcv attached to felm object.")
    }
    se.fit_mat <- sqrt(diag(m.mat %*% vcov_mat %*% t(m.mat)))
  }
  if(interval == "confidence"){
    t_val <- qt((1 - level) / 2 + level, df = object$df.residual)
    fit$lwr <- fit$fit - t_val * se.fit_mat
    fit$upr <- fit$fit + t_val * se.fit_mat
  } else if (interval == "prediction"){
    stop("interval = \"prediction\" not yet implemented")
  }
  if(se.fit){
    return(list(fit=fit, se.fit=se.fit_mat))
  } else {
    return(fit)
  }
}