基于Nimble的贝叶斯GLM模型：假设检验与残差获取问询

一、模型假设检验 #

贝叶斯框架下的假设检验核心基于后验分布推断，针对你带人口偏移项的计数GLM模型，可按以下方式操作：

1. 解释变量系数显著性推断 #

• 可信区间判断：若变量系数的95%后验可信区间不包含0，说明该变量对计数结果有显著影响。Nimble中可通过summary()提取参数分位数：

  mcmc_output <- runMCMC(model, niter = 10000, nburnin = 2000)
  summary(mcmc_output)$quantiles  # 查看各参数的分位数及95%可信区间
  

• 后验概率检验：计算系数大于/小于0的后验概率，若概率接近1或0，可认为变量显著。以Var1为例：

  var1_post_samples <- as.matrix(mcmc_output)[, "beta_var1"]
  mean(var1_post_samples > 0)  # 系数大于0的后验概率
  

2. 模型拟合优度检验 #

用**后验预测检验（PPC）**对比观测数据与后验预测数据的一致性：

• 先在Nimble模型代码中添加预测节点（以泊松模型为例）：

  modelCode <- nimbleCode({
    # 先验设定
    beta0 ~ dnorm(0, sd = 10)
    beta_var1 ~ dnorm(0, sd = 10)
    beta_var2 ~ dnorm(0, sd = 10)
    
    # 线性预测与观测模型
    for(i in 1:N){
      log(mu[i]) <- beta0 + beta_var1*Var1[i] + beta_var2*Var2[i] + log(Pop[i])
      Counts[i] ~ dpois(mu[i])
      # 添加预测节点
      Counts_pred[i] ~ dpois(mu[i])
    }
  })
  

• 运行MCMC时监控预测节点，再计算拟合统计量的后验分布：

  mcmc_output <- runMCMC(model, 
                         niter = 10000, 
                         nburnin = 2000,
                         monitors = c("beta0", "beta_var1", "beta_var2", "Counts_pred"))
  
  # 提取观测与预测数据
  obs_counts <- model$Counts
  pred_counts <- as.matrix(mcmc_output)[, grep("Counts_pred", colnames(mcmc_output))]
  
  # 计算卡方统计量并获取PPC p值
  chi_obs <- sum((obs_counts - colMeans(pred_counts))^2 / colMeans(pred_counts))
  chi_pred <- apply(pred_counts, 1, function(x) sum((obs_counts - x)^2 / x))
  ppc_p <- mean(chi_pred > chi_obs)
  

  若ppc_p落在0.05-0.95区间内，说明模型拟合效果较好。

二、获取模型残差 #

贝叶斯GLM有多种残差类型，可按需计算：

1. 原始残差 #

直接用观测值减去预测均值：

# 方式1：若模型监控了mu，直接提取后验均值
mu_post_mean <- colMeans(as.matrix(mcmc_output)[, grep("mu", colnames(mcmc_output))])
# 方式2：用参数后验均值重新计算mu
beta0_mean <- mean(as.matrix(mcmc_output)[, "beta0"])
beta_var1_mean <- mean(as.matrix(mcmc_output)[, "beta_var1"])
beta_var2_mean <- mean(as.matrix(mcmc_output)[, "beta_var2"])
mu_post_mean <- exp(beta0_mean + beta_var1_mean*Var1 + beta_var2_mean*Var2 + log(Pop))

# 计算原始残差
raw_residuals <- obs_counts - mu_post_mean

2. 皮尔逊残差 #

针对计数模型调整方差后的残差：

pearson_residuals <- (obs_counts - mu_post_mean) / sqrt(mu_post_mean)

3. 离差残差 #

基于模型离差计算，适配泊松计数模型：

deviance_residuals <- sign(obs_counts - mu_post_mean) * sqrt(2*(obs_counts*log(obs_counts/mu_post_mean) - (obs_counts - mu_post_mean)))
# 处理观测值为0的情况
deviance_residuals[obs_counts == 0] <- -sqrt(2*mu_post_mean[obs_counts == 0])

4. 后验预测残差 #

体现残差的后验不确定性：

# 计算每个MCMC样本对应的残差
post_pred_residuals <- apply(pred_counts, 1, function(x) obs_counts - x)
# 提取残差的后验均值与95%可信区间
post_pred_res_mean <- colMeans(post_pred_residuals)
post_pred_res_ci <- apply(post_pred_residuals, 2, quantile, c(0.025, 0.975))

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