有没有一种简洁的方法可以获得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
更新（2020-04-02）：这个答案是由Grant提供的，使用新软件包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的输出中。
lm2 <- lm(data = iris, Sepal.Length ~ Sepal.Width + factor(Species)) predict(lm2, newdata = data.frame(Sepal.Width = 3, Species = "virginica")) - pbaylis虽然有些晚，但新的fixest包(链接)具有可预测方法。它使用与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
> 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).
如果您非常需要使用这个包来进行预测，那么您可以通过使用对象中可用的信息编写自己简化版本的功能。例如，OLS回归系数可以通过model2$coefficients获得。
如果您希望忽略预测中的群组效应、为新 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)
}
}
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
我认为你可能需要的是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
你可能需要尝试一下来指定你想要的特定效果，但是这个软件包有很好的文档，所以不应该是问题。
Error: sum(nb) == q is not TRUE。 - kennyB
felm()函数（以及它调用的函数）进行相当多的重新设计，因为当前的实现并没有存储固定效应系数，甚至似乎也没有截距 - 参见这个答案，这是一个与此问题至少非常相似的问题。 - duckmayr