R语言中含固定效应的模型使用ggpredict报错的解决方法

问题场景 #

在使用ggpredict生成含固定效应（如国家-年份面板数据的个体/时间固定效应）的模型预测值时，出现如下错误：

Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels) : factor as.factor(FE) has new level 1.

以iris数据集为例，当模型包含Species作为固定效应时，执行ggpredict会触发该错误，移除固定效应则运行正常。

错误原因 #

模型为避免完全共线性，会自动丢弃固定效应因子的一个水平（参考水平）。但如果在模型公式中直接用as.factor()转换变量，ggpredict生成预测用的新数据时，会尝试包含原始数据的所有因子水平，导致与模型中实际保留的水平不匹配，从而报错。

解决方案 #

1. 提前转换因子变量并统一水平定义

不要在模型公式内临时转换因子，先将数据中的固定效应变量提前转为因子（可手动指定参考水平），确保模型和数据的因子水平一致：

data(iris)
# 提前将Species转为因子，默认参考水平为第一个类别（setosa）
iris$Species <- factor(iris$Species)
# 若需指定其他参考水平，用relevel
# iris$Species <- relevel(factor(iris$Species), ref = "versicolor")

# 拟合模型，直接使用转换后的因子变量
mod1 <- lm(Sepal.Length ~ Sepal.Width + Petal.Width + Species, data = iris)

# 正常生成预测值
object <- ggpredict(mod1, terms = c("Sepal.Width"), vcov.fun = "vcovHAC", ci.lvl = .95)

2. 用condition参数固定固定效应水平

如果需要针对特定固定效应水平生成预测，可通过condition参数明确指定，避免ggpredict生成模型中不存在的因子水平：

data(iris)
# 即使模型内用as.factor转换，也可通过condition指定水平
mod1 <- lm(Sepal.Length ~ Sepal.Width + Petal.Width + as.factor(Species), data = iris)

# 指定固定效应为setosa水平，生成对应预测
object <- ggpredict(mod1, terms = c("Sepal.Width"), 
                    condition = c(Species = "setosa"),
                    vcov.fun = "vcovHAC", ci.lvl = .95)

3. 面板数据固定效应适配（以plm包为例）

针对国家-年份面板数据，使用plm拟合固定效应模型时，提前将个体/时间变量转为因子，并可通过condition指定具体的固定效应水平：

library(plm)
# 假设data是包含country、year的面板数据
data$country <- factor(data$country)
data$year <- factor(data$year)

# 拟合含固定效应的模型
mod <- plm(y ~ IV + Interaction + country + year, data = data, model = "within")

# 针对指定国家和年份生成预测
object <- ggpredict(mod, terms = c("IV", "Interaction"),
                    condition = c(country = "USA", year = 2020),
                    vcov.fun = "vcovHAC", ci.lvl = .95)

内容的提问来源于stack exchange，提问作者Kajasb