如何在R中获取异方差GLS模型的预测置信区间

需要为nlme::gls（带varIdent权重的异方差模型）计算预测置信区间（而非系数置信区间），且无法使用tidymodels工作流。下面提供两种可行解决方案：

方法1：基于模型渐近理论的手动计算 #

该方法利用GLS模型的方差-协方差矩阵推导拟合值标准误，进而计算置信区间，适用于样本量较大、渐近假设成立的场景。

代码实现 #

首先拟合目标GLS模型：

library(nlme)
library(ggplot2)

# 拟合带组间异方差的GLS模型
gls_model <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), 
                 weights = varIdent(form = ~ 1|Mare),
                 data = Ovary)

计算预测置信区间：

# 提取模型拟合值
fitted_vals <- fitted(gls_model)
# 构造对应原数据的设计矩阵
X <- model.matrix(gls_model)
# 提取系数的方差-协方差矩阵
vcov_mat <- vcov(gls_model)
# 计算每个拟合值的标准误
se_fit <- sqrt(diag(X %*% vcov_mat %*% t(X)))
# 获取残差自由度
df_resid <- gls_model$df.residual
# 计算95%置信区间上下限
lower_ci <- fitted_vals - qt(0.975, df = df_resid) * se_fit
upper_ci <- fitted_vals + qt(0.975, df = df_resid) * se_fit

# 整理结果数据框
gls_ci_df <- data.frame(
  Time = Ovary$Time,
  follicles = Ovary$follicles,
  fitted = fitted_vals,
  lowerCI = lower_ci,
  upperCI = upper_ci
)

可视化结果：

ggplot(gls_ci_df, aes(x = Time)) +
  geom_point(aes(y = follicles), alpha = 0.6) +
  geom_line(aes(y = fitted), color = "#2E86AB", linewidth = 1) +
  geom_ribbon(aes(ymin = lowerCI, ymax = upperCI), fill = "#2E86AB", alpha = 0.3) +
  labs(title = "GLS模型预测值与95%置信区间", x = "时间", y = "卵泡数") +
  theme_minimal()

方法2：自助法（Bootstrap） #

若担心渐近假设不成立，可使用按组重抽样的自助法，保留原数据的组内结构，结果更稳健，但计算成本更高。

代码实现 #

set.seed(123) # 设置随机种子保证可复现
n_boot <- 1000 # 自助抽样次数
boot_fitted <- matrix(nrow = nrow(Ovary), ncol = n_boot)

# 获取所有分组标识
unique_mare <- unique(Ovary$Mare)

# 循环执行自助抽样与模型拟合
for (i in 1:n_boot) {
  # 按组有放回重抽样
  sampled_mare <- sample(unique_mare, replace = TRUE)
  boot_data <- do.call(rbind, lapply(sampled_mare, function(m) subset(Ovary, Mare == m)))
  
  # 捕获拟合失败的情况，避免中断循环
  tryCatch({
    boot_model <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), 
                      weights = varIdent(form = ~ 1|Mare),
                      data = boot_data)
    # 对原数据的自变量进行预测
    boot_fitted[, i] <- predict(boot_model, newdata = Ovary)
  }, error = function(e) {
    boot_fitted[, i] <- NA
  })
}

# 计算每个观测值的95%自助置信区间（剔除拟合失败的样本）
lower_boot_ci <- apply(boot_fitted, 1, function(x) quantile(x, 0.025, na.rm = TRUE))
upper_boot_ci <- apply(boot_fitted, 1, function(x) quantile(x, 0.975, na.rm = TRUE))

# 整理结果数据框
gls_boot_df <- data.frame(
  Time = Ovary$Time,
  follicles = Ovary$follicles,
  fitted = fitted(gls_model),
  lowerCI = lower_boot_ci,
  upperCI = upper_boot_ci
)

可视化自助法结果：

ggplot(gls_boot_df, aes(x = Time)) +
  geom_point(aes(y = follicles), alpha = 0.6) +
  geom_line(aes(y = fitted), color = "#D81E5B", linewidth = 1) +
  geom_ribbon(aes(ymin = lowerCI, ymax = upperCI), fill = "#D81E5B", alpha = 0.3) +
  labs(title = "GLS模型自助法预测95%置信区间", x = "时间", y = "卵泡数") +
  theme_minimal()

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