如何用R的marginaleffects包获取ITS分析中预测值与反事实的置信区间

我正在学习edX平台的中断时间序列分析（ITS）课程，使用R语言的nlme包构建了如下GLS模型：

# Fit the GLS regression model
model_p10 <- gls(flow ~ time + level + trend,
    data=data,
    correlation=corARMA(p=10,form=~time),
    method="ML")

目前已能计算干预后25年的预测值（pred）、反事实值（cfac）及其绝对/相对变化，但无法获取这些指标的95%置信区间，以下是针对问题的具体解决方案：

1. 用marginaleffects包获取反事实值（cfac）的95%置信区间 #

反事实值的核心是模拟未发生干预的场景，需将level和trend重置为干预前的基准状态，通过marginaleffects的predictions()函数直接生成带置信区间的结果：

library(marginaleffects)

# 构建反事实数据集：保留干预后25年时间点，将level/trend设为干预前基准值（示例为0）
cfac_data <- transform(data, level = 0, trend = 0)
cfac_data <- subset(cfac_data, time >= intervention_time & time <= intervention_time + 24)

# 生成反事实预测及95%置信区间
cfac_ci <- predictions(model_p10, newdata = cfac_data, conf_level = 0.95)
# 结果包含estimate（反事实值）、conf.low、conf.high等列
head(cfac_ci)

注：若干预变量编码非0/1，需调整cfac_data中level和trend的取值为对应基准状态

2. 通过comparisons()函数获取pred与cfac的绝对/相对变化的95%置信区间 #

comparisons()函数可直接对比干预场景与反事实场景的差异，自动计算置信区间：

绝对变化 #

abs_change <- comparisons(
    model_p10,
    newdata = subset(data, time >= intervention_time & time <= intervention_time + 24),
    variables = list(level = c(0, 1), trend = c(0, 1)),  # 0=反事实，1=真实干预
    conf_level = 0.95
)
# estimate列即为绝对变化值，附带95%置信区间上下限
head(abs_change)

相对变化（百分比） #

设置type = "ratio"即可计算相对变化，再转换为百分比：

rel_change <- comparisons(
    model_p10,
    newdata = subset(data, time >= intervention_time & time <= intervention_time + 24),
    variables = list(level = c(0, 1), trend = c(0, 1)),
    type = "ratio",
    conf_level = 0.95
)
# 转换为百分比变化
rel_change$percent_change <- (rel_change$estimate - 1) * 100
head(rel_change)

注：若干预变量编码非0/1，需对应调整variables中的取值对

3. 将level、trend的绝对系数转为带95%置信区间的相对百分比 #

根据模型类型（线性/对数转换），可通过两种方式实现：

线性模型场景 #

方法1：提取系数后手动转换

# 提取模型系数及95%置信区间
coef_ci <- tidy(model_p10, conf.int = TRUE, conf.level = 0.95)
# 筛选level和trend的系数，转换为百分比变化
coef_ci <- subset(coef_ci, term %in% c("level", "trend"))
coef_ci$percent_change <- coef_ci$estimate * 100
coef_ci$percent_low <- coef_ci$conf.low * 100
coef_ci$percent_high <- coef_ci$conf.high * 100
# 查看结果
coef_ci[, c("term", "percent_change", "percent_low", "percent_high")]

方法2：用hypotheses()直接计算转换后的置信区间

# 计算level系数的百分比变化及置信区间
hyp_level <- hypotheses(model_p10, "level * 100")
# 计算trend系数的百分比变化及置信区间
hyp_trend <- hypotheses(model_p10, "trend * 100")
# 合并结果并重命名列
percent_coef <- rbind(hyp_level, hyp_trend)
colnames(percent_coef) <- c("term", "percent_change", "std.error", "statistic", "p.value", "percent_low", "percent_high")
print(percent_coef)

对数转换响应变量场景（如log(flow)） #

需用指数转换计算百分比变化：

# 计算level系数的百分比变化及置信区间
hyp_level_log <- hypotheses(model_p10, "(exp(level) - 1) * 100")
# 计算trend系数的百分比变化及置信区间
hyp_trend_log <- hypotheses(model_p10, "(exp(trend) - 1) * 100")

内容的提问来源于stack exchange，提问作者M. Yates