1MTHM017 统计学高级专题
任务
请确保提交的作业是您自己的。这不是一项小组作业，因此不应与其他学生讨论方法和解决
方案。抄袭和与其他学生串通属于学术不端行为，将被举报。有关学术诚信的更多信息，请
点击此处。
这项任务有两个主要部分。A 部分涉及在 JAGS 中拟合不同的模型，以比较不同模型下的后验分
布。B 部分涉及使用不同方法将数据分为两组。
贝叶斯推理 [66 分］
对于贝叶斯模型拟合，您应该使用与模块问题表相同的 R 函数和语法。
1. A 部分的第一个问题涉及评估冈比亚乙型肝炎免疫计划的 有效性。数据集来自一项研究，在这项研
究中，106 名儿童接受了乙型肝炎免疫接种，并对他们的免疫水平进行了两次或三次跟踪调查。免疫
水平是通过一种称为抗乙肝病毒滴度的量来衡量的。在这项研究中，免疫滴度是在接种疫苗时测量的
（这提供了一个基线），并在所有的随访中测量。
一项类似的研究发现，抗-HB 滴度与接种疫苗后时间的倒数成正比，因此我们假设抗-HB 滴度对
数与接种疫苗后时间对数近似为线性关系。因此，我们希望拟合以下（非层次）模型
yij ∼ N (ψij , σ2 )、 (1)
ψij = α + β(tij - t¯) + γ(y0i - y¯0 )、
其中，y0i 和 yij 分别表示个体 i 的对数基线测量值和第 j 次随访对数时间测量值。此外，tij 表示 yij
对应的对数时间。最后，t¯ 和 y¯0 是 tij s 和 y0i s 的经验平均值。
csv 文件 logtitre.csv、baseline.csv 和 logtime.csv 分别包含 yij 、y0i 和 tij 的数据。
a. [7 分] 读入并总结三个数据集。然后绘制对数滴定系数测量值与对数时间的曲线图。您应该
绘制一幅图，每个孩子都有一条单独的曲线（即每个孩子的数据都应该用直线连接起来）。
评论线性假设是否合理。请注意，每当一个人只有两次随访而不是三次时，对数时间值 10 就
会被任意选作第三个时间点。在绘制曲线图时，应将这个任意选择的时间点排除在外（但也
不能完全排除，因为在 b 部分的模型拟合中还需要它）。
b. [15 分] 在 JAGS 中编码 (1) 所定义的模型，以分析对数抗-HB-滴度数据。对模型参数 α、β、
2γ 和 σ2 使用模糊先验。初始化 2 条链，并用这两条链运行模型。您必须决定 n.iter 和
burnin 的适当值。生成链的轨迹图和所有模型参数的摘要。调查所有参数的链是否收敛。
对您的发现发表评论。
请注意，这里的缺失值并不重要。之所以存在这些值，是因为有些人接受了两次随访，而有些人
则接受了三次随访。您可以让 JAGS 对这些缺失数据进行初始化。
3jags.mod <-
function(){ for(i in
1:12){
r[i] ~ dbin(theta[i],n[i])
logit(theta[i]) <-
logit.theta[i] logit.theta[i] ~
dnorm(mu, tau)
}
mu ~ dnorm(0,1.0E-3)
tau <- 1/sigmaˆ2
sigma ~ dunif(0,100)
}
c. [14 分] 由于非层次模型的截距和斜率共享的假设可能不现实，我们还可以拟合另一种模 型
，即允许不同个体有不同的截距和斜率，而这些截距/斜率在某种意义上是相关的。具体而言
，我们将考虑以下分层模型
其中
yij ∼ N (ψij , σ2 )、
ψij = αi + βi (tij - t¯) + γ(y0i - y¯0 )、
αi ∼ N (µα , ω2 )、 βi ∼ N (µβ , ω2 )、α β
µα , ω2 , µβ 和 ω2 赋予适当的先验。α β
使用模糊先验对 µα , ω2 , µβ , ω2 , σ2 和 γ 进行拟合。比较斜率和γ。α β
截取非层次模型和层次模型的参数。哪个模型能解释更大比例的方差？这背后的原因是什么？
2. 在这个问题中，我们将分析不同模型对后验分布的影响。我们将使用 surgical.csv 文件，该
文件包含 12 家医院的数据，其中 n 列给出了 1991 年 4 月至 1995 年 3 月期间各中心为 1 岁以下
儿童实施的心脏手术总数，r 列给出了这些手术中患者在术后 30 天内死亡的数量。
对于这些数据，二项模型似乎是合理的，因此我们将重点比较医院特定死亡率的不同模型。
我们将首先拟合以下形式的逻辑随机效应模型：
ri ∼ Binom(ni , θi
), logit(θi ) ∼ N (µ, σ2 ).
该模型假定各医院的失败率（θi ）在某种程度上是相似的。这种相似性表现在假设各医院的死亡
率来自相同的（母）分布。
a. [9 分] 在 JAGS 中编码该模型以分析手术数据。对共享正态分布的参数使用模糊先验。特别
是应使用以下模型定义：
用 2 个链对模型进行 10,000 次迭代，放弃前 5,000 次作为 "预演"。生成模型链的轨迹图和拟
合参数的摘要。评论所有参数的链是否收敛。
b. [8 分] 我们想知道上述随机效应模型的假设是否合理，或者是否有证据表明一家或多家医院的死
亡率并非来自同一母体分布。我们将通过比较每家医院的观察死亡人数和预测死亡人数来评估这
一点。
4i编辑您的 JAGS 模型，使其 (i) 求出下列情况中死亡人数的后验预测值 (r )pred
并 (ii) 计算每家医院的后验概率：
预言 1 预测
pi = P (ri > ri ) + 2 P (ri = ri ）。
5上述概率称为 "中间 p 值"。这使我们能够总结离散量之间的冲突。
检查中间 p 值，是否有任何医院的死亡率似乎不寻常？制作各医院预测死亡人数的核密度图
，并将其与各医院观察到的死亡人数进行直观比较，以确认您的结论。
c. [9 分] 作为替代模型，我们将考虑一个离群模型，该模型对 1-10 号和 12 号医院具有随机效应；
对 11 号医院具有独立效应：
logit(pi ) ∼ N (µ, σ2 )、 i = 1, . . ., 10, 12
logit(p11 ) ∼ N (0, 1000)
拟合该替代模型，比较随机效应模型和离群模型下医院死亡率的后验分布。解释两者的差异。(
提示：比较的最简单方法是绘制每个模型下死亡率 θ 的箱形图）。
d. [4 分] 计算离群值模型的中间 p 值，并将其与随机效应模型的中间 p 值进行比较。离群值模型是
否能更好地解释异常死亡率？
B.分类 [34 分］
下图显示了数据集 Classification.csv 中的信息--它显示了与两个解释变量相对应的两个不同组。
这是模拟数据--目的是从各种可能的方 法 中找到一种合适的方法，将 1000 个数据点分为两组。
4
2
组别
0 0
1
-2
0 2 4
X1
1. [5 分] 用变量 X1 和 X2 分别概括两组情况，并为 每组绘制这两个变量的分布图。描述你的发
现。
2. [4 分] 考虑到观测数据图、密度图和数字摘要，你认为以下哪种分类方法适合对这些数据进 行分
X2
I
6类？对于每种方法，你应解释你的答案背后的原因（即为什么给定的方法是合适的/为什么它可能
不合适）？
a. 线性判别分析
7b. 二次判别分析
c. 逻辑回归。
d. K-近邻
e. 支持向量机
3. [1 分] 选择 70% 的数据作为训练集，其余 30% 用于测试/评估。
4. [22 分] 从 Q2 中列出的方法中选择四种可能适合对数据进行分类的方法。使用这些方法进行分类
。(如果您使用的方法超过四种，则只考虑前四种方法进行评分）。在每种情况下，简要描述给定
的分类方法如何对数据进行分类，提出模型拟合度的评估结果（突出模型性能的不同方面）并描
述您的发现。在适当的情况下，优化方法的（超）参数。
5. [2 分] 比较所选四种方法的结果，选择你认为最适合本案例的分类方法，并说明理由。
试卷总分 = 100 分
提交截止日期为 7 月 31 日中午 12 点。请注意，逾期提交将受到处罚。
您应通过 eBart 提交包含问题答案（包括所有代码和相关输出/图表！）的 PDF 文件。在 A
部分，您应使用 R 编程语言，但在 B 部分，您可以选择使用 R 或 Python（或两者兼用）。