使用plm的predict()时含lag()项的回归预测异常问题

问题描述 #

使用plm包进行混合OLS回归时，若模型包含plm::lag()形式的AR(1)项，会出现模型系数与lm完全一致，但预测结果差异显著的情况：

• lm()预测：每个截面的第一期因无滞后观测值显示为NA，其余结果符合预期
• 提前生成滞后变量再传入plm的模型（plm_model2）：预测结果剔除了NA值，剩余部分与lm一致
• 公式内直接使用plm::lag()的plm模型（plm_model）：预测结果整体后移1期，且为每个截面的最后一期生成额外预测值，完全不符合预期

问题重现代码：

library(plm)

data <- data.frame(country = rep(c("A", "B", "C"), each = 3),
                  year = rep(c(1, 2, 3), 3),
                  population = c(100, 150, 200, 1000, 1200, 1250, 10, 20, 45)) |>
        pdata.frame(index = c('country', 'year'))

data$lag_population <- plm::lag(data$population, k = 1)

lm_model <- lm(population ~ lag_population, data)
plm_model <- plm(population ~ plm::lag(population, k = 1), data, model = 'pooling')
plm_model2 <- plm(population ~ lag_population, data, model = 'pooling')

预测结果差异：

predict_plm <- predict(plm_model, newdata = data, na.fill = FALSE)
predict_plm2 <- predict(plm_model2, newdata = data, na.fill = FALSE)
predict_lm <- predict(lm_model, newdata = data)

> print(predict_lm) 
       A-1        A-2        A-3        B-1        B-2        B-3        C-1        C-2        C-3 
        NA  139.50424  193.43973         NA 1110.34313 1326.08511         NA   42.42035   53.20745 
> print(predict_plm2)
       A-2        A-3        B-2        B-3        C-2        C-3 
 139.50424  193.43973 1110.34313 1326.08511   42.42035   53.20745 
> print(predict_plm)
       A-1        A-2        A-3        B-1        B-2        B-3        C-1        C-2        C-3 
 139.50424  193.43973  247.37522 1110.34313 1326.08511 1380.02060   42.42035   53.20745   80.17519 

核心原因 #

1. 公式内滞后项的处理逻辑差异：plm在拟合公式内包含plm::lag()的模型时，会自动对齐观测值（剔除滞后项缺失的行），但在调用predict()时，会对newdata重新计算滞后——它会将每个截面的最后一个观测值当作下一期的滞后项（即使该下一期不存在），导致预测整体后移并生成额外值。
2. 提前生成滞后变量的情况：plm直接使用已生成的lag_population变量，predict()会自动剔除NA对应的观测，但不会额外生成值；而lm()会保留NA位置，仅输出对应结果为NA。

解决方案 #

方法1：提前生成滞后变量并补全NA（推荐） #

先构造好滞后变量，再拟合模型，最后手动将预测结果补回原数据的NA位置：

# 获取plm_model2的预测结果
predict_plm2 <- predict(plm_model2, newdata = data, na.fill = FALSE)

# 创建与原数据长度一致的向量，初始填充NA
predict_plm2_full <- rep(NA, nrow(data))
names(predict_plm2_full) <- rownames(data)

# 将有效预测值填充到对应位置
predict_plm2_full[names(predict_plm2)] <- predict_plm2

# 此时结果与lm完全一致
print(predict_plm2_full)

方法2：手动计算公式内滞后项的预测值 #

若必须在公式内使用plm::lag()，可绕过plm的predict()函数，手动构建模型矩阵并计算预测值：

# 手动计算newdata的滞后项
newdata_lag <- plm::lag(data$population, k = 1)

# 构建模型矩阵
X <- model.matrix(~ newdata_lag, data = data)

# 用模型系数计算预测值
predict_plm_manual <- X %*% coef(plm_model)

# 为滞后项缺失的位置设置NA
predict_plm_manual[is.na(newdata_lag)] <- NA
rownames(predict_plm_manual) <- rownames(data)

# 转换为向量格式，结果与lm一致
predict_plm_manual <- as.vector(predict_plm_manual)
print(predict_plm_manual)

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