R实现LM滞后检验代码报错及PCSE使用问题咨询

研究背景与方法要求 #

• 研究主题：检验商业军事行为体（commercial military actors, CMA）的存在与军事/平民死亡数上升的相关性及因果关系，基于R语言开展线性回归分析，需确定模型应纳入的因变量滞后阶数
• 检验方法要求：采用Engle(1984)提出、Beck和Katz(1996)沿用的拉格朗日乘数（LM）检验判定滞后阶数，操作步骤为：
  1. 估计当前设定的基准模型，保存模型残差
  2. 将残差作为因变量，对残差的1期滞后项、原模型所有自变量做回归，若残差滞后项系数统计显著，说明存在剩余序列相关，需要纳入更高阶的因变量滞后项
  3. 从无滞后项的基准模型开始逐次检验，加入滞后项后重复上述步骤，直到残差滞后项系数不再显著为止
• 标准误要求：回归需报告Katz和Bailey提出的面板校正标准误（panel corrected standard errors, PCSE）
• 核心变量说明：
  • 因变量：log_military_cas，国家层面年度军事死亡数的对数转换值
  • 核心自变量：CMA，虚拟变量，国家-年份观测单元存在CMA活动记为1，无活动记为0
  • 滞后变量：lag_md，log_md的1期滞后项
  • 分析数据集：lagr

遇到的问题 #

代码报错情况 #

运行LM检验代码时连续触发两类报错，原始代码与报错信息如下：

# 拟合无滞后项基准模型
lagtest_0a <- lm(log_military_cas ~ CMA + as.factor(country) + as.factor(year), data = lagr)

# 保存残差
lagr$Risid_0 <- resid(lagtest_0)

# 第二步检验（模型设定存在错误）
lagtest_0b <- lm(log_military_cas  ~ CMA + Risid_0 + as.factor(country) + as.factor(year), data = lagr)
summary(lagtest_0b)

# 根据检验结果判定至少需要纳入1期滞后，拟合含1期滞后的模型
lagtest_1a <- lm(log_military_cas ~ CMA + lag_md + as.factor(country) + as.factor(year), data = lagr)

# 尝试保存残差到原数据集，触发第一次报错
lagr$Risid1 <- resid(lagtest_1a)
# 报错信息：
# Error in `$<-.data.frame`(`*tmp*`, Risid1, value = c(`2` = 1.84005148256506,  : 
#   replacement has 2855 rows, data has 2856

# 尝试不将残差存入原数据集，单独存储残差对象
Risid1 <- resid(lagtest_1a)

# 重新运行检验模型，触发第二次报错
lagtest1 <- lm(log_military_cas  ~ CMA + Rs_lagtest_md1 + as.factor(country) + as.factor(year), data = lagr)
# 报错信息：
# Error in model.frame.default(formula = log_military_cas ~ CMA + Rs_lagtest_md1 +  : 
#   variable lengths differ (found for 'Rs_lagtest_md1')

• 初步判断：纳入首年存在NA值的滞后变量lag_md时出现变量长度不匹配，已尝试显式设置na.action = na.omit仍未解决问题

第二个疑问 #

基准回归阶段是否需要纳入PCSE？

问题解答 #

1. 报错原因与修正方案 #

你的判断完全准确，报错核心原因是缺失值导致的模型拟合样本与原数据集行数不匹配，同时你的LM检验第二步模型设定存在错误：

• lm()函数默认会自动删除所有变量中含NA的观测行，lag_md作为1期滞后项，每个国家的第一年观测都是NA，因此lagtest_1a实际仅用了2855个观测拟合，输出的残差长度为2855，而原数据集lagr有2856行，直接往原数据框新增残差列自然会触发长度不匹配错误。
• 单独存储残差对象后，后续模型调用的残差向量长度为2855，和数据集中长度为2856的其他变量长度不一致，同样会触发报错。
• LM检验第二步的模型设定错误：按照操作指引，第二步的因变量应该是第一步模型输出的残差，而非原因变量log_military_cas，你之前的写法得到的检验结果完全无效。

修正后的LM检验代码逻辑如下：

# 提前按国家、年份排序，保证滞后项匹配正确
lagr <- lagr[order(lagr$country, lagr$year), ]

# ------------ 以1期滞后模型的LM检验为例 ------------
# 第一步：先提取模型用到的所有变量，剔除缺失值，保证拟合样本和数据行数一致
model_data <- na.omit(lagr[, c("log_military_cas", "CMA", "lag_md", "country", "year")])
# 拟合当前模型
current_model <- lm(log_military_cas ~ CMA + lag_md + as.factor(country) + as.factor(year), data = model_data)
# 提取残差
model_data$resid <- resid(current_model)
# 生成残差的1期滞后项，注意按国家分组生成滞后，避免跨国家错配
model_data$resid_lag1 <- ave(model_data$resid, model_data$country, FUN = function(x) c(NA, head(x, -1)))
# 剔除残差滞后项为NA的首年观测
lm_test_data <- na.omit(model_data)

# 第二步：做LM检验回归
lm_test <- lm(resid ~ resid_lag1 + CMA + lag_md + as.factor(country) + as.factor(year), data = lm_test_data)
summary(lm_test)
# 若resid_lag1系数显著，说明存在剩余序列相关，需要加入2期滞后项，重复上述步骤即可

2. 基准回归是否需要纳入PCSE #

需要。
PCSE是Beck & Katz针对国家-年份这类时间序列截面（面板）数据设计的标准误校正方法，专门用来处理面板回归中常见的组间异方差、组内序列相关问题，不管是不含滞后项的基准模型，还是后续确定最优滞后阶数后的最终模型，只要你需要报告系数的统计显著性，都应该使用PCSE校正后的标准误，不能直接报告普通lm输出的默认标准误。
实际操作时不需要改动模型的变量设定，拟合完基础线性模型后，调用PCSE计算函数传入模型对象即可得到校正后的标准误与显著性结果。

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