23

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

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
4

6 回答 6

15

更新（2020-04-02）：下面格兰特使用新包的答案提供了一个更简洁的解决方案。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.

于 2015-12-28T23:16:20.410 回答
9

派对迟到了，但新的fixst包（链接）有一个 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
于 2019-11-08T20:52:57.940 回答
6

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

该对象与“lm”对象有一些相似之处，并且为 lm 设计的一些后处理方法可能会起作用。然而，可能需要强制对象成功完成此操作。

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

lfe 包旨在在非常大的数据集上运行，并努力节省内存：作为直接结果，与felm对象相反，对象不使用/包含 qr 分解lm。不幸的是，该lm predict过程依赖于这些信息来计算预测。因此，强制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.

于 2015-07-05T14:05:26.100 回答
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)
  }
}
于 2016-03-04T23:41:27.237 回答
3

为了扩展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
于 2019-05-26T09:02:29.303 回答
-1

我认为您正在寻找的可能是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

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

于 2015-07-02T18:22:35.047 回答