贝叶斯诊断试验参数估计代码报错求助及敏感性分析需求

问题概述 #

尝试对AGID、P-ELISA、H-ELISA三种诊断试验的真实患病率、灵敏度、特异度进行贝叶斯估计时，模型运行持续失败，错误提示：

Error: unexpected '{' in: "model <- model {"

同时需要开展敏感性分析。

现有代码（含语法错误） #

# 数据定义存在语法错误，缺少data列表结构与关键观测数据
c(199, 89, 68, 78, 69, 74, 85, 500),
  agid_sens <- c(41.4, 59.8),
  pelisa_sens <- c(92.7, 100.0),
  helisa_sens <- c(90.9, 99.4),
  agid_spec <- c(98.4, 99.9),
  pelisa_spec <- c(95.1, 99.0),
  helisa_spec <- c(91.4, 96.6),
  prev <- c(0.01, 0.07)
)

# Define the model
model <-
model {
  # Prior distributions
  prev ~ dunif(0, 1)
  agid_sens ~ dunif(0, 100)
  agid_spec ~ dunif(0, 100)
  pelisa_sens ~ dunif(0, 100)
  pelisa_spec ~ dunif(0, 100)
  helisa_sens ~ dunif(0, 100)
  helisa_spec ~ dunif(0, 100)

  # Likelihood
  for (i in 1:length(n)) {
    agid_pos[i] ~ dbin(p_agid[i], n[i])
    pelisa_pos[i] ~ dbin(p_pelisa[i], n[i])
    helisa_pos[i] ~ dbin(p_helisa[i], n[i])
    p_agid[i] <- prev * agid_sens / (prev * agid_sens + 
                                    (1 - prev) * (1 - agid_spec))
    p_pelisa[i] <- prev * pelisa_sens / (prev * pelisa_sens + 
                                        (1 - prev) * (1 - pelisa_spec))
    p_helisa[i] <- prev * helisa_sens / (prev * helisa_sens + 
                                        (1 - prev) * (1 - helisa_spec))
  }
}

# Run the model
samples <- bugs(
  data = data,
  inits = NULL,
  parameters.to.save = c("prev", "agid_sens", "agid_spec", "pelisa_sens", 
            "pelisa_spec", "helisa_sens", "helisa_spec"),
  model.file = textConnection(model),
  n.chains = 3,
  n.iter = 10000,
  n.burnin = 5000,
  n.thin = 1,
  DIC = FALSE
)

# Summarize the results
summary(samples)

代码错误修正 #

1. 模型定义语法错误
   原代码中model <- model { ... }不符合R语法，需将模型代码以字符串形式赋值给model变量，用引号包裹整个模型块：

   model <- "
   model {
     # 模型内容...
   }
   "
   

2. 数据结构修正
   需将输入数据整理为标准列表，包含样本量n、各试验阳性观测数agid_pos/pelisa_pos/helisa_pos，注意不要将先验范围混入数据：

   data <- list(
     n = c(199, 89, 68, 78, 69, 74, 85, 500),
     agid_pos = c(/* 填入各样本组AGID阳性数 */),
     pelisa_pos = c(/* 填入各样本组P-ELISA阳性数 */),
     helisa_pos = c(/* 填入各样本组H-ELISA阳性数 */)
   )
   

3. 概率单位匹配
   模型中灵敏度/特异度以0-100的百分比表示，但二项分布dbin要求概率为0-1的小数，需在计算阳性预测值时转换单位：

   p_agid[i] <- prev * (agid_sens/100) / (prev * (agid_sens/100) + 
                                   (1 - prev) * (1 - agid_spec/100))
   

4. 变量定义完整性
   确保数据中包含模型用到的所有变量（如n），且观测数据长度与n一致。

敏感性分析方案 #

• 先验分布敏感性：将均匀先验dunif替换为有信息的先验（如基于文献的dbeta分布），对比参数估计结果的波动程度。
• 样本量敏感性：拆分数据集为不同子集，测试样本量对患病率、灵敏度/特异度估计稳定性的影响；或通过模拟不同样本量的数据集验证结果鲁棒性。
• 模型结构敏感性：尝试加入试验间的相关性结构（如分层模型），或调整似然函数设定，观察参数估计的差异。

内容的提问来源于stack exchange，提问作者sunday O. ochai