R语言使用rethinking包从后验分布估算发病率的代码问题求助

• 变量名冗余：预先定义的fev、no.fev变量未被实际使用，属于无效代码，若要复用变量需保证模型传入的变量名和定义一致
• 语法错误：duration <- rnorm(1e3, 7, 3).行末尾多余英文句号，属于语法错误会直接终止代码运行
• 先验设置不符合业务假设：明确假设发热平均病程为7天，但病程模型的mu先验设置为均值14天的正态分布，和假设冲突会拉偏后验结果
• 后验样本计算逻辑错误：extract.samples()返回的是包含所有模型参数的data.frame，直接对两个完整数据框做除法不符合计算逻辑，需要分别提取患病率参数p、病程均值参数mu的样本值参与计算
• 单位未统一：直接按患病率/病程计算得到的是每日发病率，需要乘以365天转换为你需要的每人每年发病率

library(rethinking)

# 发热患病率模型（过去2周调查数据）
fever_cases <- 132
no_fever_cases <- 675
total <- fever_cases + no_fever_cases

fever.p1 <- quap(
  alist(
    fever ~ dbinom(total, p),
    p ~ dbeta(10, 90)
  ), 
  data = list(fever = fever_cases)
)
precis(fever.p1)
# 提取患病率p的后验样本
sample_fev_p <- extract.samples(fever.p1, 1e5)$p

# 发热病程模型
flist <- alist(
  duration ~ dnorm(mu, sigma),
  mu ~ dnorm(7, 3),  # 先验均值调整为和7天病程假设一致
  sigma ~ dunif(0, 7)
)
# 移除多余句号修正语法错误
duration <- rnorm(1e3, 7, 3)
d <- data.frame(duration = duration)

fev.d <- quap(flist, data = d)
precis(fev.d)
# 提取病程均值mu的后验样本
sample_fev_duration_mu <- extract.samples(fev.d, 1e5)$mu

# 计算每人每年发病率：每日发病率*365天
fever_incidence_per_year <- (sample_fev_p / sample_fev_duration_mu) * 365

# 查看后验结果
precis(fever_incidence_per_year)
# 绘制后验分布密度图
dens(fever_incidence_per_year, main = "发热年发病率后验分布", xlab = "每人每年发病率")

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