R语言Causal包相对效应计算异常问题求助

使用R语言CausalImpact包结合bsts模型计算相对效应时，出现结果远超预期、不符合实际的负值问题。以下是原始代码及结果图，结合常见问题点给出排查和修复方案：

原始代码 #

# Model 3
ss3 <- list()
# Semi Local trend, weekly-seasonal
ss3 <- AddSemilocalLinearTrend(ss3, ts_dash_wk)
# Add weekly seasonal
ss3 <- AddSeasonal(ss3, ts_dash_wk, nseasons = 52)
model3 <- bsts(ts_dash_wk,
               state.specification = ss3,
               niter = 1500,
               burn = 500)
plot(model3, main = "Model 3")
plot(model3, "components")
# Causal impact of dash and Covid-19 
pre.period <- as.Date(c("2019-01-01", "2020-03-11"))
post.period <- as.Date(c("2020-03-11", "2020-12-31"))

# Obtain post period data
dat_dash_causal_post <- data_dash_wk %>% 
  filter(week >= as.Date("2020-03-11"))

# Use model 3 for causal impact
impact <- CausalImpact(bsts.model = model3,
                       post.period.response = dat_dash_causal_post$dash, alpha = 0.05)
plot(impact)

summary(impact)

结果图 #

排查与修复方案 #

1. 时间序列与模型成分匹配检查 #

• 周季节性参数验证：nseasons=52仅适用于年度周度频率的时间序列，需确认ts_dash_wk的时间索引是严格的每周频率（可通过frequency(ts_dash_wk)查看，应为52）。若数据频率错误，季节性成分拟合会完全偏离，导致预测值异常。
• 趋势成分调整：半局部线性趋势（SemilocalLinearTrend）适合趋势缓慢变化的序列，若你的数据存在突变趋势（如疫情前后的大幅波动），建议替换为局部线性趋势（LocalLinearTrend），或调整趋势的先验方差（如AddSemilocalLinearTrend(..., level.sigma.prior = ExponentialPrior(0.1))），让模型更适应趋势变化。

2. 前后周期划分修正 #

• 避免周期重叠：当前pre.period和post.period的起始/结束日期均为2020-03-11，会导致数据点被重复计算，引发模型混淆。建议调整为无重叠的周期：

  pre.period <- as.Date(c("2019-01-01", "2020-03-08")) # 疫情前最后一周
  post.period <- as.Date(c("2020-03-15", "2020-12-31")) # 疫情开始后第一周
  

• 时间索引对齐：确保pre.period和post.period的日期与ts_dash_wk的时间戳完全匹配，避免因日期偏移导致模型错误映射观测值与预测值。

3. 模型收敛性优化 #

• 增加MCMC迭代次数：当前niter=1500、burn=500的迭代次数可能不足以让模型收敛，建议提升至：

  model3 <- bsts(ts_dash_wk,
                 state.specification = ss3,
                 niter = 3000,
                 burn = 1000)
  

• 收敛性诊断：运行plot(model3, "diagnostics")查看参数迹图，确保所有参数的链稳定无漂移；运行plot(model3, "residuals")检查残差是否为白噪声，若残差存在未解释的趋势/季节性，需补充模型成分（如节假日效应）。

4. CausalImpact调用方式优化 #

• 避免手动提取后周期数据：直接使用CausalImpact的原生参数传递数据，减少手动筛选出错概率：

  impact <- CausalImpact(data = ts_dash_wk,
                         pre.period = pre.period,
                         post.period = post.period,
                         model.args = list(state.specification = ss3,
                                           niter = 3000,
                                           burn = 1000),
                         alpha = 0.05)
  

• 非负响应变量适配：若dash是计数类非负变量，在bsts中指定分布族，避免预测出负值：

  model3 <- bsts(ts_dash_wk,
                 state.specification = ss3,
                 niter = 3000,
                 burn = 1000,
                 family = "poisson") # 或"negbin"处理过度离散
  

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