R包官网:
https://cecileproust-lima.github.io/lcmm/index.html
本次教程地址:
https://cecileproust-lima.github.io/lcmm/index.html
https://cecileproust-lima.github.io/lcmm/articles/lcmm.html
简介
利用潜在分类和潜在过程(Latent Classes and Latent Processes)的扩展混合模型LCMM
lcmm致力于对混合模型进行多种形式的扩展。它可以处理具有重复测量和潜在分类的连续(正态或非正态)和有序结果。同时还可以在比例风险模型中处理生存结果。所有这些模型都是在最大似然框架下使用了改良的Marquardt-Levenberg算法所估算出的。该软件包还包括多个预测、可视化以及用于进行完整的统计分析的其他实用功能。
tips
潜在类别:
在医学数据中,我们常遇到来自不同病人群体的信息,但这些群体并不总是显而易见的。 潜在类别是指这些隐藏在数据背后的不同群体或类别。例如,在研究病人对某种治疗的反应时,虽然所有病人看似接受了同样的治疗,但他们的反应可能因隐藏的不同群体特性而有所不同。 通过潜在分类分析,我们可以识别出这些不同的病人群体,并更好地理解不同群体对治疗的响应方式。
潜在过程:
在医学研究中,我们关注的某些关键变量(如病情进展)可能不是直接观测到的,而是通过其他可观测数据间接反映出来。 潜在过程指的是这些影响观测数据的隐藏或不可见的动态变化。例如,一个病人的健康恢复过程可能是一个潜在过程,它影响着临床上观察到的各种健康指标,但本身并不直接被记录。 通过分析这些潜在过程,医生和研究者可以更深入地了解病情的发展和恢复过程,从而为病人提供更有针对性的治疗方案。
关于此软件包的详细介绍已在《统计软件杂志》发表; Proust-Lima C, Philipps V, Liquet B. 利用潜在类别和潜在过程估计扩展混合模型:R软件包lcmm。《统计软件杂志》,文章。2017;78(2):1-56. https://doi.org/10.18637/jss.v078.i02
作者在各种统计论文中还描述了特定的统计模型的估计方法。
安装R包
lcmm需要R软件的3.5或更新版本。要安装CRAN发布的版本,请使用以下命令:
install.packages("lcmm")
要获取最新更新的版本,请从github安装:
remotes::install_github("CecileProust-Lima/lcmm")
lcmm依赖于其他一些R软件包,具体如下:
survival (>=2.37-2) 用于处理生存分析结果 parallel 和 doParallel 用于并行计算一些耗时函数 mvtnorm 用于生成随机参数 randtoolbox 用于准蒙特卡罗序列(quasi Monte Carlo sequences) marqLevAlg (>2.0) 用于数值优化 numDeriv 用于计算海森矩阵(Hessian)
为了运行后续的示例,还需要以下软件包:
lattice NormPsy ggplot2 ggpubr dplyr splines gridExtra ggalluvial
让我们荡起双桨
软件包的功能
lcmm基于线性混合模型理论,提供了一系列函数,用于估算统计模型。包括以下:
hlme
:混合模型以及潜在类别混合模型,用于高斯纵向结果lcmm
:曲线型有序单变量纵向结果multlcmm
:曲线型多变量结果jlcmm
和mpjlcmm
:联合潜在类别混合模型,用于(高斯和/或曲线)纵向结果和可能左截断右删失并在竞争环境中定义的生存结果。
tips:上文部分概念的解释
混合模型(Mixed Model):这种模型结合了固定效应(整体趋势)和随机效应(个体差异),从而能够同时考虑群体水平和个体水平的变化。 高斯纵向结果(Gaussian longitudinal outcomes ):指的是数据遵循高斯(正态)分布的纵向(随时间)测量结果。 曲线型(Curvilinear):数据随时间呈现非直线性变化。 有序(Ordinal):数据是按顺序排列的,但不是连续的。 左截断右删失的生存结果:指的是事件发生的时间可能不完整地被观测到,例如,研究开始前已经发生的事件(左截断)或研究结束前尚未发生的事件(右删失)。 在竞争环境中定义:指的是可能有多个不同的事件或结果相互竞争,影响着最终观测到的时间至事件。
lcmm包括5个主要函数:lcmm
、hlme
、multlcmm
、jlcmm
和mpjlcmm
。每个函数在不同的vignette中有所描述。这些函数适用于纵向标记,即使是那些计量属性有限的标记,但函数也可以用于横截面情境。
对于每个模型,我们使用改进的Marquardt算法获得最大似然估计值。这个算法有严格的收敛标准,这些标准基于参数的稳定性、似然函数的稳定性,以及二阶导数的负值。简单来说,我们通过这些标准确保模型的精确性和可靠性。这个算法在R软件包marqLevAlg中实现,并且支持并行计算,这意味着可以更快地处理复杂的计算任务。
tips:上文部分概念的解释
最大似然估计:这是一种统计方法,用于在给定观测数据的情况下,估计模型参数的最优值。简单来说,它帮助我们找到使得观测到的数据最有可能发生的参数值。 改进的Marquardt算法:这是一种用于求解非线性最优化问题的算法。在这里,它被用来优化最大似然估计,即寻找最能解释数据的模型参数。 严格的收敛标准:这意味着算法在找到最优解之前会严格检查一系列条件,比如参数的稳定性和似然函数的稳定性,以及二阶导数的负值。这些都是确保找到的解既准确又可靠的重要条件。
此外,lcmm还提供了各种拟合后功能,包括拟合优度分析、分类、图形展示、预测轨迹、事件的个体动态预测以及预测准确性评估。这些功能对于理解模型的效果和应用模型到实际数据中非常有用。
如果您在使用过程中有任何疑问或问题,可以在lcmm的GitHub页面上提出:https://github.com/CecileProust-Lima/lcmm/issues。这里是一个共享问题和寻找解答的平台,可以帮助您更好地利用这个软件包。
函数调用
每个函数在配套论文(Proust-Lima, JSS 2017 - https://doi.org/10.18637/jss.v078.i02)中都有详细介绍。
hlme
此函数用于标准线性混合模型(standard linear mixed models ),以及扩展到多元轨迹轮廓 (multiple profiles of trajectory)的潜在类别线性混合模型(the latent class linear mixed models)。
tips:什么是多元轨迹轮廓
多元轨迹轮廓Multiple profiles of trajectory指的是数据中存在多种不同的发展轨迹或模式。在统计和医学研究中,这通常意味着研究对象(如病人)随时间的变化模式不是单一的,而是可以被分为几个不同的群体,每个群体有其独特的变化轨迹。 例如,在研究一种疾病的进展时,不同的病人可能表现出不同的病程轨迹——有的可能迅速恶化,有的可能长期稳定,还有的可能有所改善。这些不同的轨迹(或“轨迹轮廓”)就是多元轨迹轮廓的实例。
hlme(fixed, mixture, random, subject, classmb, ng = 1, idiag = FALSE, nwg = FALSE, cor = NULL, data, B, convB = 0.0001, convL = 0.0001, convG = 0.0001, prior, maxiter = 500, subset = NULL, na.action = 1, posfix = NULL)
参数解释
fixed: 指定模型的固定效应部分。此处需要定义模型的主要公式,包括响应变量和预测变量。 mixture: 用于指定混合效应的模型部分,即在不同的潜在类别中可能有所不同的模型参数。 random: 指定模型中的随机效应部分。这涉及到对每个观察单位或个体特有的随机变异的建模。 subject: 指明用于区分不同观察单位(如个体)的变量。 classmb: 用于潜在类别模型的成员资格(membership)预测的公式。 ng: 指定潜在类别的数量,默认为1。 idiag: 一个逻辑值,决定是否对随机效应的协方差矩阵使用对角化简化,默认为FALSE。 nwg: 一个逻辑值,用于决定是否在不同的潜在类别之间对随机效应的协方差结构进行区分,默认为FALSE。 cor: 指定模型中的残差相关结构,NULL意味着不考虑残差相关。 data: 用于分析的数据集。 B: 初始参数估计。 convB, convL, convG: 分别指定参数估计、似然值和梯度的收敛标准。 prior: 先验信息,用于贝叶斯估计。 maxiter: 最大迭代次数,默认为500。 subset: 用于分析的数据子集。 na.action: 缺失值处理方法,默认为1。 posfix: 用于固定某些参数的位置。
lcmm
这个函数用于估算结果不遵循高斯分布(正态分布)的混合模型。该方法被称为“潜在过程模型”。目前,它能够使用连续的链接函数处理曲线或非高斯连续标记;同时也采用分段常数链接函数的Probit框架来处理二元和有序标记。与hlme
函数一样,这种方法也能处理在“潜在过程和潜在类别混合模型”中的不同轨迹轮廓。
tips:上文部分内容和概念的解释
处理曲线形或非高斯连续标记,使用连续的链接函数
曲线型标记指的是随时间变化呈现曲线形状(非线性)的数据 非高斯连续标记则是指那些连续但不符合高斯(正态)分布的数据 在这两种情况下,函数使用连续的链接函数来将模型的线性预测值转换为与这些特殊数据类型相符的形式 链接函数是统计模型中用于连接模型的线性部分与预测变量之间关系的工具。
同时也处理二元和有序标记,采用分段常数链接函数的Probit框架
二元标记通常指的是只有两个可能结果的数据(例如,是/否,成功/失败) 有序标记则指的是有明确顺序但不一定等距的分类数据(例如,评级为低、中、高) 在这两种情况下,函数使用Probit框架中的分段常数链接函数 Probit模型是一种处理二元或有序响应变量的常用统计模型 分段常数链接函数则允许这些分类数据在不同区间内呈现不同的趋势或模式。
lcmm(fixed, mixture, random, subject, classmb, ng = 1, idiag = FALSE, nwg = FALSE, link = "linear", intnodes = NULL, epsY = 0.5, cor = NULL, data, B, convB = 1e-04, convL = 1e-04, convG = 1e-04, maxiter = 100, nsim = 100, prior, range = NULL, subset = NULL, na.action = 1, posfix = NULL, partialH = FALSE)
参数解释
link: 链接函数类型,决定了模型中自变量与因变量之间的关系。 intnodes: 用于分段常数链接函数的内部节点,影响模型处理非线性数据的方式。 epsY: 设置估计过程中使用的精度,尤其对连续标记的处理重要。 nsim: 模拟研究中的模拟次数,用于某些统计方法。 range: 定义连续协变量的取值范围。 partialH: 是否使用部分海森矩阵来提高计算效率。
其余参数,在之前的函数中已经解释。
multlcmm
扩展lcmm函数:
当有多个标记测量同一潜在结构(例如,多个测试用于评估同一心理特质)时,multlcmm
函数可以对lcmm
函数的功能进行拓展。
处理不同性质的标记:
与lcmm
一样,multlcmm
也能处理曲线或非高斯连续标记(使用连续链接函数),以及二元和有序标记(使用Probit框架下的分段常数链接函数)。
二元/有序标记和连续时间纵向IRT模型:
当只包含二元或有序标记时,multlcmm
定义了一个时间连续的纵向IRT(项目反应理论)模型。IRT模型是用于是一种用于设计、分析和评分测验的统计框架,常用于分析问卷或测试数据,尤其适用于心理学和教育测量领域。
处理同质或异质的轨迹轮廓:
类似于其他函数,multlcmm
能够处理不同类型的轨迹轮廓,无论是均质的(所有个体或群体遵循相同的轨迹模式)还是异质的(不同个体或群体有不同的轨迹模式)。
tips:变量、标记、结果这三个词的区别
变量(Variable):就像问题中的不同选项。在一项研究或调查中,变量是你收集的不同类型的信息,比如人们的年龄、性别或喜欢的颜色。
标记(Marker):特别的提示。在研究中,标记是一种特殊的变量,它可以帮助你理解或预测某些情况,比如用血液中的某种物质水平来预测疾病。
结果(Outcome):最终想要了解的事情。在一个实验或研究中,结果是你最关心的部分,比如你想知道某种药物是否能治愈病人。
multlcmm(fixed, mixture, random, subject, classmb, ng = 1, idiag = FALSE, nwg = FALSE, randomY = FALSE, link = "linear", intnodes = NULL, epsY = 0.5, cor = NULL, data, B, convB = 1e-04, convL = 1e-04, convG = 1e-04, maxiter = 100, nsim = 100, prior, range = NULL, subset = NULL, na.action = 1, posfix = NULL, partialH = FALSE)
参数解释
在multlcmm
函数中,未在前文提到的参数有:
randomY: 逻辑值,决定是否为每个标记(结果)添加随机效应。这对于模型的灵活性很重要,特别是在处理多个标记测量同一潜在构造时。
intnodes: 用于分段常数链接函数的内部节点,这有助于定义处理非线性关系时的关键转折点。
epsY: 设置连续标记的测量误差精度。较小的值表示更高的测量精度。
nsim: 指定模拟研究中的模拟次数。在使用蒙特卡罗方法估计模型参数时,这个值决定了模拟的次数,影响模型估计的准确性和可靠性。
range: 为连续协变量定义一个特定的范围,这对于确保模型预测值保持在合理的范围内非常有用。
partialH: 逻辑值,指示是否使用部分海森矩阵以提高计算效率。在处理复杂模型时,这可以显著减少计算时间。
部分参数在前文已经提及。
jlcmm
该函数扩展了hlme
、lcmm
和multlcmm
函数的功能,对竞争环境下的生存结果进行联合分析。该函数特别之处在于它将纵向结果(如连续观测数据)和生存解雇(如疾病复发或死亡时间)通过潜在类别结构联系起来。
Jointlcmm(fixed, mixture, random, subject, classmb, ng = 1, idiag = FALSE, nwg = FALSE, survival, hazard = "Weibull", hazardtype = "Specific", hazardnodes = NULL, TimeDepVar = NULL, link = NULL, intnodes = NULL, epsY = 0.5, range = NULL, cor = NULL, data, B, convB = 1e-4, convL = 1e-4, convG = 1e-4, maxiter = 100, nsim = 100, prior, logscale = FALSE, subset = NULL, na.action = 1, posfix = NULL, partialH = FALSE)
参数解释
survival: 用于指定与时间至事件数据相关的模型部分。这是生存分析的关键组成部分,用于建模时间至事件数据。 hazard: 指定风险函数的类型,例如"Weibull"。风险函数是生存分析中用来描述事件发生率随时间变化的函数。 hazardtype: 指定风险类型,比如"Specific",表示每个潜在类别都有其特定的风险函数。 hazardnodes: 用于定义风险函数中的内部节点,这对于处理风险随时间变化的非线性模式很重要。 TimeDepVar: 用于指定时间依赖的协变量。 link: 链接函数的类型,用于生存模型部分。 logscale: 逻辑值,指示是否在对数尺度上进行风险函数的估计。
部分参数上文已经提及。
mpjlcmm
此函数是jlcmm
的多变量版本:它可以建模多变量纵向标记(这些纵向标记最后可能被归纳为几个潜在过程),并可以对生存结果进行联合分析(在竞争环境下)。潜在过程(即纵向结果)和生存结果都被潜在类别所影响。
mpjlcmm <- function(longitudinal,subject,classmb,ng,survival,
hazard="Weibull",hazardtype="Specific",hazardnodes=NULL,TimeDepVar=NULL,
data,B,convB=0.0001,convL=0.0001,convG=0.0001,maxiter=100,nsim=100,
prior,logscale=FALSE,subset=NULL,na.action=1,posfix=NULL,
partialH=FALSE,verbose=TRUE,nproc=1,clustertype=NULL)
参数解释
longitudinal: 定义纵向数据模型的部分,特别指定如何建模随时间变化的数据。 hazardnodes: 用于定义风险函数中的内部节点,这有助于在处理风险随时间变化的非线性模式时定义关键转折点。 TimeDepVar: 指定时间依赖的协变量,这对于包含时间相关效应的模型而言非常重要。 verbose: 逻辑值,决定函数运行时是否输出详细信息,有助于调试和模型拟合过程中的监控。 nproc: 指定用于并行计算的处理器数量。这在处理大型数据集或进行计算密集型模型估计时非常有用。 clustertype: 在并行计算时使用的集群类型,它定义了并行处理的具体方式,这对于优化计算性能和资源使用非常重要。
部分参数上文已经提及。
下面是栗子
paquid子样本
在lcmm
包的wiki文档中,使用了一个名为paquid的数据集作为示例。这个数据集是从法国Paquid前瞻性研究的原始数据中选取了500名受试者的数据作为子样本。需要特别指出的是,这个子样本并不完全代表原始研究队列的整体特征,特别是因为痴呆症(dementia)病例在这个子样本中被过度采样,因此这份数据集不适宜用于进行流行病学研究。
该数据集的格式为纵向数据,涵盖了包括三项心理测试(简易精神状态检查MMSE
、本顿视觉记忆测试BVRT
和智力测验IST
)、抑郁症状量表CESD
在内的多个变量。此外,还包括年龄age
、直到痴呆症发作的年龄agedem
、是否患有痴呆症的标志dem
(布尔型变量,1代表患病)、入组时的年龄ageinit
、教育水平CEP
以及性别male
等信息。
数据可视化(仅显示表格的头部信息)
head(paquid)
ID MMSE BVRT IST HIER CESD age agedem dem age_init CEP male
1 1 26 10 37 2 11 68.50630 68.5063 0 67.4167 1 1
2 2 26 13 25 1 10 66.99540 85.6167 1 65.9167 1 0
3 2 28 13 28 1 15 69.09530 85.6167 1 65.9167 1 0
4 2 25 12 23 1 18 73.80720 85.6167 1 65.9167 1 0
5 2 24 13 16 3 22 84.14237 85.6167 1 65.9167 1 0
6 2 22 9 15 3 NA 87.09103 85.6167 1 65.9167 1 0
在这个数据集中,不同的标记是在不同的时间点收集的。而在这个数据集里,时间的度量是以年龄为标准的(而不是按照实际的日历日期或研究开始后的时间长度)。
快速获得数据的summary
summary(paquid)
ID MMSE BVRT IST
Min. : 1.0 Min. : 0.00 Min. : 0.00 Min. : 5.00
1st Qu.:132.2 1st Qu.:25.00 1st Qu.: 9.00 1st Qu.:22.00
Median :263.0 Median :27.00 Median :11.00 Median :27.00
Mean :256.2 Mean :25.99 Mean :10.78 Mean :26.52
3rd Qu.:376.0 3rd Qu.:29.00 3rd Qu.:13.00 3rd Qu.:31.00
Max. :500.0 Max. :30.00 Max. :15.00 Max. :40.00
NA's :36 NA's :300 NA's :198
HIER CESD age agedem
Min. :0.000 Min. : 0.000 Min. : 66.28 Min. :66.87
1st Qu.:1.000 1st Qu.: 2.000 1st Qu.: 75.09 1st Qu.:82.03
Median :1.000 Median : 6.000 Median : 80.57 Median :86.61
Mean :1.272 Mean : 8.488 Mean : 80.53 Mean :85.89
3rd Qu.:2.000 3rd Qu.:12.000 3rd Qu.: 85.77 3rd Qu.:89.94
Max. :3.000 Max. :52.000 Max. :102.16 Max. :99.49
NA's :46 NA's :146
dem age_init CEP male
Min. :0.0000 Min. :65.25 Min. :0.0000 Min. :0.0000
1st Qu.:0.0000 1st Qu.:68.00 1st Qu.:0.0000 1st Qu.:0.0000
Median :0.0000 Median :71.92 Median :1.0000 Median :0.0000
Mean :0.3329 Mean :73.01 Mean :0.7373 Mean :0.3853
3rd Qu.:1.0000 3rd Qu.:76.96 3rd Qu.:1.0000 3rd Qu.:1.0000
Max. :1.0000 Max. :92.33 Max. :1.0000 Max. :1.0000
在lcmm中使用的某些变量可能包含缺失值,不过问题不大,因为lcmm会系统地移除缺失值。
MMSE结果
在lcmm包的示例中,MMSE
(简易精神状态检查)通常被视为研究的结果。MMSE
是一种非常常见的神经心理测试,用于衡量老年人的全面认知功能。由于MMSE
分布非常不对称,因此在应用于高斯变量的方法时通常需要进行标准化处理。这可以通过在NormPsy包中提供的专门用于MMSE
的预标准化函数来完成。
library(NormPsy)
# 标准化MMSE数据:
paquid$normMMSE <- normMMSE(paquid$MMSE)
# 该代码使用normMMSE函数对paquid数据集中的MMSE数据进行标准化处理
# 注意,normMMSE函数在这个R包中专门用来处理MMSE数据集,使其正态化
# 处理后的数据集储存在paquid数据集的新列normMMSE上
# 设置绘图参数
par(mfrow=c(1,2))
# par(mfrow=c(1,2)):设置图形参数,使得接下来的两个直方图并排显示
# mfrow=c(1,2)表示一个1行2列的图形布局,即两个图形并排显示
# 绘制原始和标准化后的MMSE数据的直方图:
hist(paquid$MMSE, cex.main=0.8, cex.lab=0.8)
# 这行代码绘制了原始MMSE数据的直方图
# cex.main和cex.lab参数用于调整图形标题和标签的字体大小
hist(paquid$normMMSE, cex.main=0.8, cex.lab=0.8)
# 这行代码绘制了标准化后的MMSE数据的直方图
要建模的MMSE
(简易精神状态检查)的个体重复测量数据包括:
library(lattice)
color <- paquid$ID
xyplot(normMMSE ~ age, paquid, groups = ID, col=color, lwd=2, type="l")