从 lfe 包中预测 felm 的方法

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 包 (link) 有一个 predict 方法。它使用与 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

这可能不是您正在寻找的答案，但似乎作者没有向lfe 包添加任何功能，以便使用拟合的felm 模型对外部数据进行预测。主要关注点似乎是对组固定效应的分析。然而，有趣的是，在包的文档中提到了以下内容：

该对象与“lm”对象有一些相似之处，有些 为 lm 设计的后处理方法可能会起作用。有可能 但是有必要强制对象成功。

因此，可以将felm 对象强制转换为lm 对象以获得一些额外的lm 功能（如果对象中存在执行必要计算所需的所有信息）。

lfe 包旨在在非常大的数据集上运行，并努力节省内存：作为直接结果，felm 对象不使用/包含 qr 分解，而不是 @987654327 @ 目的。不幸的是，lmpredict 过程依赖此信息来计算预测。因此，强制 felm 对象并执行 predict 方法将失败：

> model2 <- felm(data = iris, Sepal.Length ~ Sepal.Width | Species)
> class(model2) <- c("lm","felm") # coerce to lm object
> predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica"))
Error in qr.lm(object) : lm object does not have a proper 'qr' component.
 Rank zero or should not have used lm(.., qr=FALSE).

如果你真的必须使用这个包来执行预测，那么你可以使用 felm 对象中提供的信息来编写你自己的简化版本的这个功能。例如，OLS 回归系数可通过model2$coefficients 获得。

Answer 4

为了扩展pbaylis 的答案，我创建了一个稍微冗长的函数，该函数很好地扩展以允许多个固定效果。请注意，您必须手动输入 felm 模型中使用的原始数据集。该函数返回一个包含两项的列表：预测向量和基于 new_data 的数据帧，其中包含预测和固定效果作为列。

predict_felm <- function(model, data, new_data) {

  require(dplyr)

  # Get the names of all the variables
  y <- model$lhs
  x <- rownames(model$beta)
  fe <- names(model$fe)

  # Demean according to fixed effects
  data_demeaned <- demeanlist(data[c(y, x)],
                             as.list(data[fe]),
                             na.rm = T)

  # Create formula for LM and run prediction
  lm_formula <- as.formula(
    paste(y, "~", paste(x, collapse = "+"))
  )

  lm_model <- lm(lm_formula, data = data_demeaned)
  lm_predict <- predict(lm_model,
                        newdata = new_data)

  # Collect coefficients for fe
  fe_coeffs <- getfe(model) %>% 
    select(fixed_effect = effect, fe_type = fe, idx)

  # For each fixed effect, merge estimated fixed effect back into new_data
  new_data_merge <- new_data
  for (i in fe) {

    fe_i <- fe_coeffs %>% filter(fe_type == i)

    by_cols <- c("idx")
    names(by_cols) <- i

    new_data_merge <- left_join(new_data_merge, fe_i, by = by_cols) %>%
      select(-matches("^idx"))

  }

  if (length(lm_predict) != nrow(new_data_merge)) stop("unmatching number of rows")

  # Sum all the fixed effects
  all_fixed_effects <- base::rowSums(select(new_data_merge, matches("^fixed_effect")))

  # Create dataframe with predictions
  new_data_predict <- new_data_merge %>% 
    mutate(lm_predict = lm_predict, 
           felm_predict = all_fixed_effects + lm_predict)

  return(list(predict = new_data_predict$felm_predict,
              data = new_data_predict))

}

model2 <- felm(data = iris, Sepal.Length ~ Sepal.Width | Species)
predict_felm(model = model2, data = iris, new_data = data.frame(Sepal.Width = 3, Species = "virginica"))
# Returns prediction and data frame

Answer 5

这应该适用于您希望忽略预测中的组效应、预测新 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)
  }
}

Answer 6

我认为您正在寻找的可能是lme4 包。我能够使用这个来预测工作：

library(lme4)
data(iris)

model2 <- lmer(data = iris, Sepal.Length ~ (Sepal.Width | Species))
predict(model2, newdata = data.frame(Sepal.Width = 3, Species = "virginica"))
       1 
6.610102 

您可能需要花一点时间来指定您正在寻找的特定效果，但该软件包有据可查，因此应该不是问题。