全文链接：https://tecdat.cn/?p=36112
原文出处：拓端数据部落公众号
淋巴结疾病作为一类复杂的健康问题，其风险预测一直是临床和公共卫生领域的研究热点。随着统计学的进步和计算能力的提升，广义线性模型（GLM）成为了分析这类数据的有力工具。特别是当数据呈现比例特性时，beta二项分布作为广义线性模型的一个特例，为我们提供了一种灵活且强大的方法来建模和预测淋巴结疾病的风险。
在统计分析中，beta二项分布模型因其能够处理比例数据且能够同时考虑数据的均值和变异程度而备受关注。特别地，在医疗领域，当数据表现为成功次数与总次数的比例时，beta二项分布模型能够为我们提供一个强大的工具来分析这些比例数据的内在规律。
本文旨在帮助客户探讨beta二项分布模型在淋巴结疾病风险预测中的应用，并通过R语言实现数据的可视化分析。通过这一研究，我们期望能够为淋巴结疾病的风险预测提供新的视角和方法，并为临床医生和患者提供更有价值的决策支持。
数据集
本项研究使用的数据集包含17,659个观测值，涉及六个变量。展示了数据集的前六行，以便对数据的结构和内容有一个初步的了解。
head(data)数据集前六行概览
CS.lymph.nodes表示淋巴结计数，Regional.nodes.examined和Regional.nodes.positive分别表示被检查的淋巴结数量和阳性淋巴结数量。此外，数据集还包括种族(Race.recode)、年龄(Age.recode.with.<1.year.olds)和性别(Sex)等分类变量。
数据结构描述
使用str()函数对数据结构进行进一步分析
数据结构描述
CS.lymph.nodes：数值型变量，表示淋巴结计数，共有17,659个观测值。Regional.nodes.examined：数值型变量，表示被检查的淋巴结数量，共有17,659个观测值。Regional.nodes.positive：数值型变量，表示阳性淋巴结数量，共有17,659个观测值。Race.recode：字符型变量，表示种族信息，共有17,659个观测值。Age.recode.with.<1.year.olds：字符型变量，表示年龄信息，共有17,659个观测值。Sex：字符型变量，表示性别信息，共有17,659个观测值。
模型拟合
为了探究淋巴结数量及其阳性比例与可能的影响因素之间的关系，我们使用了beta二项分布模型进行拟合。beta二项分布模型是一种广义线性模型，特别适用于比例数据的建模。具体地，我们使用vglm()函数进行模型拟合，将Regional.nodes.examined和Regional.nodes.positive作为响应变量，betabinomial作为分布函数，data作为数据集，并设置trace = TRUE以显示迭代过程。模型拟合的代码如下：
egional.nodes.positive) ~ 1, betabinomial,在模型拟合完成后，我们将对模型结果进行详细的分析和解释。
这表明模型在第五次迭代时收敛，并且对应的对数似然值为-33041.2181。这一数值将作为后续模型评估的重要依据。
模型结果分析
为了获取模型参数的估计值，我们使用了coef()函数，并设置了matrix = TRUE以矩阵形式展示结果：
这里的
logit(mu)和logit(rho)分别表示beta二项分布中均值（mu）和分散参数（rho）的logit变换后的估计值。为了更容易地解释这些值，我们进一步将这些logit变换后的系数转换回原始尺度：
转换后的系数值如下：
其中，mu的估计值表示在给定条件下，阳性淋巴结数量的期望比例，而rho的估计值则描述了这一比例的变异程度。在本例中，mu的估计值约为0.67，表明在平均情况下，大约67%的被检查淋巴结会被诊断为阳性。同时，rho的估计值约为0.79，说明在数据集中，阳性淋巴结比例的变化相对较小，这可能与数据集的特性或潜在的影响因素有关。
在模型拟合完成后，为了进一步理解模型中的数据和权重，我们结合了模型的响应变量（fit@y）和先验权重（weights(fit, type = "prior")）进行了初步的数据查看。以下是合并后数据的头部展示：
其中，第一列（[,1]）表示模型响应变量的观测值（例如，阳性淋巴结的比例），而第二列（[,2]）则代表与这些观测值相对应的先验权重。这些权重在模型拟合过程中用于调整不同观测值对模型参数估计的影响。
模型结果汇总
接下来，我们对模型的整体拟合结果进行了汇总。以下是使用summary(fit)函数得到的输出：
模型系数估计
summary(fit)模型估计了两个参数：logit(mu)和logit(rho)的截距项。mu代表阳性淋巴结的期望比例，而rho则描述了这一比例的变异程度。从上述结果可以看出，两个截距项的估计值分别为0.71364和1.29980，且均具有较高的统计显著性（p值远小于0.05）。
模型包含两个线性预测器，分别对应于logit(mu)和logit(rho)。模型的对数似然值为-33041.22，在32764个自由度下计算得出。此外，模型在5次迭代后收敛。
数据可视化
为了更直观地展示模型结果，我们使用ggplot2包绘制了散点图，并通过geom_point和geom_smooth函数添加了数据点和拟合曲线。在散点图中，我们以Regional.nodes.positive为x轴，以模型响应变量fit@y（即模型预测的阳性淋巴结比例）为y轴。通过这一可视化手段，我们可以清晰地看到数据点的分布情况以及模型拟合的效果。
geom_point(aes(text = paste("fit@y:", fit@y)))
淋巴结疾病数据
在本研究中，我们来拟合一个beta二项分布模型，以分析淋巴结疾病数据。首先，我们构建了一个仅包含截距项的模型，用于估计阳性淋巴结（R）与阴性淋巴结（N-R）之间的比例关系。模型构建如下：
在模型迭代过程中，我们观察到了对数似然值（log-likelihood）的变化，最终模型在第3次迭代时收敛，并给出了稳定的对数似然值-122.69531。这表示模型已找到了最优的参数估计值。
模型的参数估计结果如下：
这里，logit(mu)和logit(rho)分别表示阳性淋巴结比例的期望（mu）和变异程度（rho）的对数几率。通过反变换，我们可以得到这些参数的原始值：
Coef(fit)接下来，为了进一步提高模型的预测能力，我们引入了额外的解释变量。具体地，我们将患者分组（grp）和血红蛋白水平（hb）作为预测因子加入到模型中，并设置beta二项分布的零参数为2，以适应数据的特定特征。新的模型构建如下：
fgrp + hb 在模型迭代过程中，我们观察到了对数似然值的稳定下降，最终模型在第8次迭代时收敛，并给出了稳定的对数似然值-92.907851。这表明，通过引入额外的解释变量，我们得到了一个更加拟合数据的模型。
模型参数估计的解读
经过在beta二项分布模型中引入患者分组（fgrp）和血红蛋白水平（hb）作为预测因子，我们得到了新的模型fit2。通过执行coef(fit2, matrix = TRUE)，我们得到了模型参数的估计值，具体如下：
从上述结果中，我们可以看到不同患者分组（fgrp2、fgrp3、fgrp4）相对于基准组（fgrp1）对阳性淋巴结比例（mu）的对数几率的影响。特别地，fgrp2、fgrp3、和fgrp4的系数均为负数，表明这些组的患者相较于基准组具有较低的阳性淋巴结比例。此外，血红蛋白水平（hb）的系数为负，说明随着血红蛋白水平的降低，阳性淋巴结的比例也随之降低。
对于变异程度（rho），我们注意到除了截距项外，其他预测因子的系数均为0。这可能是由于在模型构建过程中，我们没有为rho引入额外的预测因子，或者这些因子对rho的影响不显著。
数据可视化分析
为了更直观地理解模型结果，我们对数据进行了可视化分析。首先，我们使用plot函数绘制了血红蛋白水平（hb）与死亡比例（R / N）之间的散点图，其中不同的患者分组用不同的颜色和符号表示。图表的标题为“拟合值（线条）”。
smalldf = with(lirat, lirat[N > 1, ])接下来，为了在每个患者分组中绘制模型的拟合线，我们首先筛选出N > 1的数据子集smalldf。然后，我们遍历每个患者分组，提取出该组的血红蛋白水平（hb）和模型拟合值（fitted(fit2)），并使用lines函数在散点图上绘制出拟合线。
通过这一步骤，我们可以清楚地看到在每个患者分组中，模型如何根据血红蛋白水平的变化来预测阳性淋巴结的比例。这些拟合线不仅有助于验证模型的准确性，还能为我们提供关于患者分组和血红蛋白水平对疾病严重程度影响的深入理解。