R语言动态回归：如何在回归模型中纳入时间维度

数据预处理 #

首先得把数据整理成面板时间序列格式，确保每个地点（UGA）的时间连续，生成促销变量的滞后项（对应t-1、t-2期的促销活动）。

# 加载依赖包
library(dplyr)
library(lubridate)
library(plm)
library(dynlm)
library(ggplot2)

# 数据清洗与转换
df <- df %>%
  # 把日期转成月度粒度，方便按月份分组
  mutate(Date = floor_date(Date, "month")) %>%
  # 按地点、时间排序，保证时间序列连续性
  arrange(UGA, Date) %>%
  # 为每个地点生成促销变量的滞后1期、滞后2期
  group_by(UGA) %>%
  mutate(
    Nb_Mails_lag1 = lag(Nb_Mails, 1),
    Nb_Mails_lag2 = lag(Nb_Mails, 2),
    Nb_Visits_lag1 = lag(Nb_Visits, 1),
    Nb_Visits_lag2 = lag(Nb_Visits, 2),
    Nb_Events_lag1 = lag(Nb_Events, 1),
    Nb_Events_lag2 = lag(Nb_Events, 2)
  ) %>%
  ungroup() %>%
  # 删掉滞后操作产生的缺失行
  drop_na()

核心模型实现 #

1. 静态面板分布滞后模型 #

这是最直接的方案：把当期、滞后1期、滞后2期的促销变量全部纳入回归，同时控制地点固定效应（消除不同地点的固有销量差异）和月度固定效应（消除季节性影响）。

# 拟合固定效应分布滞后模型
dlm_model <- plm(
  Nb_Sales ~ Nb_Mails + Nb_Mails_lag1 + Nb_Mails_lag2 +
             Nb_Visits + Nb_Visits_lag1 + Nb_Visits_lag2 +
             Nb_Events + Nb_Events_lag1 + Nb_Events_lag2 +
             factor(Date),  # 月度时间固定效应
  data = df,
  index = c("UGA", "Date"),  # 面板数据标识：个体（UGA）+时间（Date）
  model = "within"  # 固定效应模型，自动控制地点个体差异
)

# 查看回归结果
summary(dlm_model)

# 计算某类促销的累计效应（当期+滞后1+滞后2）
mails_total_effect <- coef(dlm_model)["Nb_Mails"] + coef(dlm_model)["Nb_Mails_lag1"] + coef(dlm_model)["Nb_Mails_lag2"]
cat("邮件促销3期累计效应：", round(mails_total_effect, 4), "\n")

2. 动态面板GMM模型 #

如果销量存在自相关（比如本月销量受上月影响），或者促销和销量存在双向因果（比如销量好的区域会加做促销），用系统GMM模型可以解决内生性问题，同时纳入滞后销量和滞后促销项。

# 拟合两步系统GMM模型
gmm_model <- pgmm(
  Nb_Sales ~ lag(Nb_Sales, 1) +  # 滞后1期销量，控制自相关
             Nb_Mails + Nb_Mails_lag1 + Nb_Mails_lag2 +
             Nb_Visits + Nb_Visits_lag1 + Nb_Visits_lag2 +
             Nb_Events + Nb_Events_lag1 + Nb_Events_lag2 +
             factor(Date) |
             lag(Nb_Sales, 2:9),  # 用滞后2-9期的销量做工具变量
  data = df,
  index = c("UGA", "Date"),
  model = "twostep"
)

summary(gmm_model)

3. 多项式分布滞后模型（PDL） #

如果想简化滞后项的参数数量，或者平滑时间效应，可以用PDL把多个滞后项转换成多项式项，适合滞后阶数较多的场景。

# 拟合多项式分布滞后模型，这里设置滞后3期（当期+滞后1+滞后2）
pdl_model <- dynlm(
  Nb_Sales ~ factor(UGA) + factor(Date) +
             PDL(Nb_Mails, 3) + PDL(Nb_Visits, 3) + PDL(Nb_Events, 3),
  data = df
)

summary(pdl_model)

结果可视化与诊断 #

• 绘制滞后系数的柱状图，直观展示促销效应随时间的变化：

# 提取邮件促销的各期系数
mails_coef_df <- data.frame(
  time_lag = c("当期", "滞后1月", "滞后2月"),
  coefficient = c(
    coef(dlm_model)["Nb_Mails"],
    coef(dlm_model)["Nb_Mails_lag1"],
    coef(dlm_model)["Nb_Mails_lag2"]
  )
)

ggplot(mails_coef_df, aes(x = time_lag, y = coefficient)) +
  geom_col(fill = "#3498db", width = 0.6) +
  labs(title = "邮件促销的分时间点效应", x = "时间滞后", y = "回归系数") +
  theme_minimal()

• 用plmtest包做异方差、序列相关检验，确保模型稳健性：

library(plmtest)
# 异方差检验
plmtest(dlm_model, type = "bp")
# 序列相关检验
pbgtest(dlm_model)

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