如何在survregbayes中实现城市层面聚类的稳健标准误？

我使用2012-2022年7903个城市的面板数据，通过生存模型研究城市行政数字化进程：

• 定义t0为区间起始、t2为数字化发生年份或面板结束年份，二元变量treat_v标记事件是否发生
• PRO_COM为城市唯一标识
• 自变量包括总人口（total_ln）、平均收入（Avg_income_ln）、行政软件支出（software_ln）

数据集结构如下：

PRO_COM  year treat_v    t0    t2 total_ln Avg_income_ln software_ln
     <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>         <dbl>       <dbl>
 1    1001  2012       0     0     1     7.88          9.95        13.6
 2    1001  2013       0     1     2     7.91          9.96         0  
...

我需要使用加速失效时间模型（AFT）和空间加速失效时间模型，采用spBayesSurv包中的survregbayes函数，请问如何在此模型中实现按城市聚类的稳健标准误？现有代码如下：

install.packages('spBayesSurv')
library(spBayesSurv)
mcmc <- list(nburn = 5000, nsave = 10000, nskip = 5, ndisplay= 2000)
prior <- list(maxL = 15)

res_prova <- survregbayes(Surv(t0, t2, treat_v) ~ total_ln + Avg_income_ln + software_ln
                      , data=df_prova, 
                      survmodel="AFT", dist="lognormal",InitParamMCMC=T, mcmc=mcmc,prior=prior)
summary.survregbayes(res_prova)

可复现数据集示例：

df_prova <- structure(list(PRO_COM = c(1001, 1001,...),...))

由于spBayesSurv基于贝叶斯框架，频率学派的“聚类稳健标准误”思路不直接适用，需通过以下两种贝叶斯范式的方法处理城市内的相关性：

1. 纳入城市水平随机截距（最直接的聚类稳健方案） #

在模型公式中加入城市标识PRO_COM作为随机截距，通过分层模型捕捉城市间未观测异质性，同时天然控制城市内观测的聚集效应：

# 修改模型公式，添加(1|PRO_COM)引入城市随机截距
res_prova <- survregbayes(
  Surv(t0, t2, treat_v) ~ total_ln + Avg_income_ln + software_ln + (1|PRO_COM),
  data=df_prova, 
  survmodel="AFT", dist="lognormal",
  InitParamMCMC=T, mcmc=mcmc,prior=prior
)
summary.survregbayes(res_prova)

此时模型输出的后验标准差已经考虑了城市内的聚类相关性，无需额外调整。

2. 基于MCMC样本手动计算聚类稳健后验标准误 #

如果需要贴近频率学派的聚类稳健逻辑，可以提取MCMC样本后，按城市分组计算调整后的标准误：

# 提取固定效应的MCMC样本
mcmc_samples <- as.matrix(res_prova$mcmc$beta)

# 按城市分组计算得分矩阵
score_list <- by(df_prova, df_prova$PRO_COM, function(sub_df) {
  X <- model.matrix(~ total_ln + Avg_income_ln + software_ln, data = sub_df)
  t(X) %*% X %*% mcmc_samples
})

# 聚合得分并计算聚类稳健协方差矩阵
cluster_scores <- do.call(cbind, score_list)
robust_cov <- t(cluster_scores) %*% cluster_scores / nrow(df_prova)
robust_se <- sqrt(diag(robust_cov))

# 输出结果
cbind(
  Posterior_Mean = colMeans(mcmc_samples),
  Robust_SE = robust_se
)

注意：该方法需确保MCMC样本收敛良好，仅适用于固定效应部分的标准误调整。

3. 空间AFT模型的聚类处理 #

若使用空间AFT模型，需先构建城市空间权重矩阵W（基于地理邻接/距离），再结合城市随机截距同时处理空间依赖和聚类相关性：

# 假设已构建空间权重矩阵W
prior_spatial <- list(maxL = 15, tau = list(prior="gamma", shape=1, rate=0.01))

res_spatial <- survregbayes(
  Surv(t0, t2, treat_v) ~ total_ln + Avg_income_ln + software_ln + (1|PRO_COM),
  data=df_prova, 
  survmodel="SAFT", dist="lognormal",
  W=W, # 传入空间权重矩阵
  InitParamMCMC=T, mcmc=mcmc,prior=prior_spatial
)
summary.survregbayes(res_spatial)

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