glmmTMB中AR(1)相关参数的标准误及置信区间提取问题

核心原理 #

glmmTMB中AR(1)结构的相关系数(ρ)通过逆logit变换从无约束theta参数（记为θ₂）推导而来：ρ = plogis(θ₂)，θ₂的标准误可从模型完整方差协方差矩阵中提取。我们需要用delta方法将θ₂的标准误转换为ρ的标准误，进而计算置信区间。

一、基础AR(1)模型（g2）的处理 #

1. 定位并提取目标theta参数 #

首先查看模型的theta参数集合：

summary(g2)$theta

AR(1)结构包含两个theta参数：第一个是随机效应的对数标准差，第二个是AR(1)相关系数的logit变换值（即θ₂）。

接着提取该theta的方差与标准误：

# 提取完整方差协方差矩阵
vcov_full <- vcov(g2, full = TRUE)
# 定位AR(1)相关系数对应的theta列
theta_col <- grep("theta_time\\+0\\|group\\.2", colnames(vcov_full))
theta_est <- summary(g2)$theta[theta_col]
theta_se <- sqrt(diag(vcov_full)[theta_col])

2. 计算AR(1)相关系数的标准误与置信区间 #

# 相关系数点估计
rho_est <- plogis(theta_est)
# delta方法计算标准误
rho_se <- rho_est * (1 - rho_est) * theta_se

# 95%正态近似置信区间
rho_ci <- rho_est + c(-1, 1) * 1.96 * rho_se

# 输出结果
cat("AR(1)相关系数估计值：", rho_est, "\n")
cat("标准误：", rho_se, "\n")
cat("95%置信区间：", rho_ci[1], "-", rho_ci[2], "\n")

二、含额外随机效应的模型（g3）的处理 #

1. 定位AR(1)相关系数对应的theta参数 #

查看模型theta参数：

summary(g3)$theta

结合你提供的vcov输出，theta_time+0|group.2对应AR(1)相关系数的logit变换值（θ₂），theta_time+0|group.1为随机效应的对数标准差。

2. 提取参数并转换为相关系数指标 #

# 提取完整方差协方差矩阵
vcov_full_g3 <- vcov(g3, full = TRUE)
# 定位目标theta列
theta_col_g3 <- grep("theta_time\\+0\\|group\\.2", colnames(vcov_full_g3))
theta_est_g3 <- summary(g3)$theta[theta_col_g3]
theta_se_g3 <- sqrt(diag(vcov_full_g3)[theta_col_g3])

# 计算相关系数、标准误与置信区间
rho_est_g3 <- plogis(theta_est_g3)
rho_se_g3 <- rho_est_g3 * (1 - rho_est_g3) * theta_se_g3
rho_ci_g3 <- rho_est_g3 + c(-1, 1) * 1.96 * rho_se_g3

# 输出结果
cat("AR(1)相关系数估计值：", rho_est_g3, "\n")
cat("标准误：", rho_se_g3, "\n")
cat("95%置信区间：", rho_ci_g3[1], "-", rho_ci_g3[2], "\n")

3. 稳健的参数定位方法（不依赖列名） #

若模型结构变化导致列名改变，可通过随机效应组件结构定位：

# 获取模型随机效应结构信息
re_struct <- g3$modelInfo$reTrms
# 找到AR(1)结构对应的组件
ar1_comp <- which(sapply(re_struct$cnms, function(x) any(grepl("ar1", x))))
# AR(1)结构的第二个theta参数即为logit(ρ)
theta_idx <- sum(re_struct$thetaDim[1:ar1_comp])
theta_est_g3 <- summary(g3)$theta[theta_idx]
theta_se_g3 <- sqrt(diag(vcov(g3, full=TRUE))[theta_idx])

补充说明 #

• confint()默认仅返回固定效应置信区间，需指定parm = "theta"获取theta参数的置信区间，再用plogis()转换为ρ的置信区间：

  theta_ci <- confint(g3, parm = "theta")
  rho_ci <- plogis(theta_ci[grep("theta_time\\+0\\|group\\.2", rownames(theta_ci)), ])
  

• 样本量较小时，建议用参数自助法提升置信区间稳健性：

  # 加载boot包（需提前安装）
  library(boot)
  
  # 定义自助抽样函数
  boot_fun <- function(data, indices) {
    d_boot <- data[indices,]
    mod <- glmmTMB(value ~ (1|time) + ar1(time + 0 | group), data = d_boot)
    plogis(summary(mod)$theta[grep("theta_time\\+0\\|group\\.2", names(summary(mod)$theta))])
  }
  
  # 运行1000次自助抽样
  boot_res <- boot(data = d, statistic = boot_fun, R = 1000)
  # 获取百分位数法95%置信区间
  boot_ci <- boot.ci(boot_res, type = "perc")
  print(boot_ci)
  

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