R语言LCMM多维度潜在类别模型流行病学研究：LCA、MM方法分析纵向数据

全文代码数据：https://tecdat.cn/?p=39710

在数据分析领域，当我们面对一组数据时，通常会有已知的分组情况，比如不同的治疗组、性别组或种族组等（点击文末“阅读原文”获取完整代码数据）。

然而，数据中还可能存在未被观测到的分组，例如素食者与非素食者、经常锻炼者与不锻炼者，或者有某种疾病家族史与无家族史的人群（假设这些数据未被收集）。针对这种情况，有相应的分析方法，本文将帮助客户围绕潜在类别分析（Latent Class Analysis，LCA）和混合建模（Mixture Modeling，MM，有时也指有限混合建模Finite MM）展开讨论，并通过实例展示其在纵向数据中的应用。

模型与实现

数据处理挑战

在使用ProcTraj和Mplus环境时，发现处理并非所有人都在完全相同时间点进行观测的数据较为棘手。在设计严格的研究中，所有受试者在基线、6个月和12个月时进行测量，这些软件使用起来较为直接。但在观察性流行病学研究中，人们的诊断年龄不同，有人离开研究（区域）或去世导致数据收集时间不同，或者选择某个时间点分析数据时，受试者可用数据量不同，此时这些软件处理这类数据的方式并不直观。

实例分析

数据集介绍

我们有如下数据集：

# 展示数据集结构> str(dat)

# 展示数据集头部

 

其中，x为自变量（可视为年龄），y为因变量（如胆固醇水平），id用于标识数据来自哪个受试者，rf1和rf2是两个风险因素，后续可能会纳入分析，totobs表示每个受试者的观测次数。

数据可视化

数据量较大，有超过2000个受试者，每个受试者最多有80个时间点。首先绘制所有数据，结果如图1所示。

图1 所有数据图
有人可能认为第一个图（左上角）没有展示出有效信息，但实际上可以明显看出右上角部分的数据密度远高于左上角，即更多人最终的值高于起始值。接下来的三个图展示了300个随机选择且至少有5个时间点的受试者数据，以避免数据过于杂乱。右上角是原始数据图，左下角是为每个人拟合的直线，右下角是使用ggplot2软件包中geom_smooth()函数默认设置为每个人拟合的平滑曲线（loess）。进一步去除杂乱数据后，更能突出图右侧与左侧相比数据分布的增加。

---

点击标题查阅往期内容

R语言用潜类别混合效应模型(Latent Class Mixed Model ,LCMM)分析老年痴呆年龄数据

左右滑动查看更多

01

02

03

04

平滑算法展示

这里简单直观地展示loess平滑算法的工作原理。图2展示了3个有较多数据的受试者，包括他们的原始（锯齿状）图、平滑曲线以及二者结合，突出该算法如何准确描绘趋势并去除噪声，使我们能轻松识别趋势路径。

图2 loess平滑算法示例

模型构建与分析

用户能够指定希望算法找到的潜在类别（未观测到的组）数量。以下代码用于在数据中寻找2个潜在类别：

# 寻找2个潜在类别result2 <- lcmm(y_variable ~ x

这里使用自变量x_variable对因变量y_variable进行建模，允许x_variable因人而异（随机斜率），并使用（线性）x_variable作为混合项中的变量。目前先保持简单，后续会考虑更合适的模型。在输出结果中，从摘要的主要部分可以看到，“x 1”和“x 2”行表明每个类别有不同的斜率，一个略为正，一个略为负，且截距（见“intercept 2”）也有很大差异。
同时，模型还提供了拟合优度统计量。其中，组隶属的后验概率很重要，每个受试者被分配到每个（2个）类别的概率如下：

# 查看后验概率

 

可以看到一个类别中有62%的受试者，另一个有38%。属于类别1的受试者，其属于类别1的平均后验概率为0.93，属于类别2的平均后验概率为0.09。类似地，属于类别2的受试者，其属于类别2的平均后验概率为0.90，属于类别1的平均后验概率为0.07。通过以下代码进一步查看，最小值接近0.5，说明一些受试者的轨迹可能确实属于任意一个类别，两个下四分位数都大于0.80。

# 进一步查看后验概率统计

 

 

由于处理的是趋势和轨迹，我们可能希望可视化其对数据的意义。通过以下代码，提取受试者所属类别，并根据类别对数据进行绘图：

结果如图3所示。

图3 2个潜在类别数据图
从图中可以看到预期结果，首先，对一个类别内所有数据进行平滑处理时，线相交，表明类别确实产生了实际影响。如果看到两条平行线，可能意味着每个人都有相同“类型”的轨迹，那么类别可能只是算法为了“找到两个类别”而产生的人为结果。这里特别的是，可以看到类别2中的一组受试者在x达到200时，y值迅速上升，而类别1中的受试者y值通常保持稳定。
尝试3个类别时，结果如图4所示。

图4 3个潜在类别数据图
尝试4个类别时，结果如图5所示。

图5 4个潜在类别数据图
在选择合适的类别数量时，普遍接受的方法是使用贝叶斯信息准则（BIC），而不是赤池信息准则（AIC），BIC值越低，模型越好。对于此数据，4个类别的BIC为32941，3个类别的BIC为32972，2个类别的BIC为33084，这表明4个类别比3个类别拟合得稍好。实际上，5个类别能使BIC更低，但如同所有建模一样，需要人为（或到2035年非常智能的AI机器人）干预。在4个类别中，第4个类别仅包含10%的受试者，在本研究的数据子集中，即25人。需要根据数据情况决定是否合适，这里认为不应低于10%。

模型优化

如前所述，线性效应可能不是对该数据建模的最佳方式。这里选择3个类别，并引入二次项，结果如图6所示。与之前3个潜在类别的图相比，有一些细微差异。

图6 含二次项3个潜在类别数据图
在考虑后验概率和类别成员百分比时，仅含线性项的模型和含二次项的模型有很大差异。似乎二次项在确定每个受试者属于哪个类别时更困难。具体如下：
仅含线性项

含二次项

然而，BIC实际上从32972降至32594，这一差异比改变类别数量时看到的差异大得多。在这种情况下，一方面需要进一步探索模型；另一方面，如果必须做出选择，会倾向于仅基于线性项的模型，因为其在确定受试者所属类别时，后验概率表现优于基于BIC的模型 。

结论

未来的研究可以在此基础上，进一步探索不同的模型设定和参数选择，以更好地适应各种复杂的数据结构和研究问题。例如，可以尝试结合其他变量或考虑不同的混合分布形式，以提高模型的解释能力和预测精度。此外，对于处理纵向数据中缺失值的方法，也可以进行更深入的研究和比较，以确保分析结果的可靠性和稳定性。同时，随着数据量的不断增大和数据维度的增加，如何提高算法的计算效率和可扩展性也是值得关注的问题。

本文中分析的数据、代码分享到会员群，扫描下面二维码即可加群！ 

---

资料获取

在公众号后台回复“领资料”，可免费获取数据分析、机器学习、深度学习等学习资料。

点击文末“阅读原文”

获取全文完整代码数据资料。

本文选自《R语言LCMM多维度潜在类别模型流行病学研究：LCA、MM方法分析纵向数据》。

点击标题查阅往期内容

R语言贝叶斯广义线性混合（多层次/水平/嵌套）模型GLMM、逻辑回归分析教育留级影响因素数据

R语言估计多元标记的潜过程混合效应模型（lcmm）分析心理测试的认知过程

R语言因子实验设计nlme拟合非线性混合模型分析有机农业施氮水平

R语言非线性混合效应 NLME模型(固定效应&随机效应)对抗哮喘药物茶碱动力学研究

R语言用线性混合效应（多水平/层次/嵌套）模型分析声调高低与礼貌态度的关系

R语言LME4混合效应模型研究教师的受欢迎程度

R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model分析藻类数据实例

R语言混合线性模型、多层次模型、回归模型分析学生平均成绩GPA和可视化

R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例

R语言用lme4多层次（混合效应）广义线性模型（GLM），逻辑回归分析教育留级调查数据

R语言 线性混合效应模型实战案例

R语言混合效应逻辑回归（mixed effects logistic）模型分析肺癌数据

R语言如何用潜类别混合效应模型（LCMM）分析抑郁症状

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言建立和可视化混合效应模型mixed effect model

R语言LME4混合效应模型研究教师的受欢迎程度

R语言 线性混合效应模型实战案例

R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言如何解决线性混合模型中畸形拟合(Singular fit)的问题

基于R语言的lmer混合线性回归模型

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

R语言分层线性模型案例

R语言用WinBUGS 软件对学术能力测验（SAT）建立分层模型

使用SAS，Stata，HLM，R，SPSS和Mplus的分层线性模型HLM

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

SPSS中的多层（等级）线性模型Multilevel linear models研究整容手术数据

用SPSS估计HLM多层（层次）线性模型模型