R软件用潜在类别混合模型LCM分析老年人抑郁数据轨迹多变量建模研究

全文链接： tecdat.cn/?p=40283

潜在类别混合模型假设总体具有异质性，由 GG 个潜在类别组成。在多变量的情况下，潜在类别是根据 KK 个纵向结果来定义的，从而形成 GG 个组，每个组的特征由 KK 个轨迹均值轮廓集表示。

多变量潜在类别混合模型

潜在类别成员身份由一个离散随机变量 cici 定义，如果个体 ii 属于潜在类别 gg（g=1,…,Gg=1,…,G），则 cici 等于 gg。变量 cici 是潜在的，其概率使用多项逻辑模型根据协变量 XciXci 进行描述：

其中，ξ0gξ0g 是类别 gg 的截距，ξ1gξ1g 是与时间无关的协变量 XciXci 的 q1q1 维向量相关的类别特定参数的 q1q1 维向量。为了可识别性，ξ0G=0ξ0G=0 且 ξ1G=0ξ1G=0。当没有协变量预测潜在类别成员身份时，该模型简化为类别特定概率。
对于 k=1,…,Kk=1,…,K 的每个 YkYk 的轨迹是在潜在类别条件下定义的。对于高斯结果，在类别 gg 的条件下，模型是一个线性混合模型，对于个体 ii 在时刻 jj 定义为：

其中，X2kijX2kij、X3kijX3kij 和 ZkijZkij 分别是与类别间公共固定效应 βkβk、类别特定固定效应 γkgγkg 以及个体随机效应 bki|ci=gbki|ci=g（称为 bkigbkig，其分布现在是类别特定的）相关的协变量向量。X2kX2k 和 X3kX3k 不能有共同的变量。
随机效应和误差测量在不同结果之间不相关。因此，在潜在类别条件下，KK 个结果是独立的。
对于曲线结果，我们使用由以下公式定义的潜在过程模型：

其中，HkHk 是由 ηkηk 参数化的链接函数。H−1H−1 可以属于线性函数族、重新缩放的 Beta 累积分布函数或二次 I - 样条函数族。然而，需要注意的是，mpjlcmm 函数目前仅支持连续结果，因此 IRT 模型暂时不可用。

后验分类

在涉及潜在类别的模型中，可以对每个潜在类别中的个体进行后验分类。它基于类成员概率的后验计算，用于描述个体的分类情况以及评估模型的拟合优度。
后验类成员概率使用贝叶斯定理计算，即给定所收集的全部信息，个体属于某个潜在类别的概率。在纵向模型中，对于个体 ii 和潜在类别 gg，其定义为：

其中，θ^Gθ^G 是在 GG 个潜在类别模型中估计的参数向量。

双变量示例

在这个示例中，我们同时研究了一组老年人（入组时年龄在 65 岁及以上）在长达 15 年的随访期间，认知标志物（MMSE）和抑郁量表（CESD）的轨迹。这两个结果都具有偏态分布，因此我们将使用具有 I - 样条链接函数的潜在过程模型。

考虑的模型

我们考虑不进行进一步调整的、与年龄相关的类别特定线性轨迹。对于类别 gg、个体 ii 和重复测量 jj，模型如下：

其中，ukig∼N(0,Bkg)ukig∼N(0,Bkg) 且 ϵij∼N(0,σ2)ϵij∼N(0,σ2)。

估计单类别模型（G=1G=1）

为了估计多变量模型，我们首先使用适当的函数定义每个单变量子模型。由于这里使用的是潜在过程模型，我们使用 lcmm 函数。
我们从 MMSE 开始：

# 定义 MMSE 的单变量模型
mmse_model <- lcmm(MMSE ~ I((age - 65) / 10), random = ~I((age - 65) / 10), subject = 'ID', data = paquid,
 link = "5 - quant - splines", verbose = FALSE)
# 查看模型摘要
summary(mmse_model)

 

 

在摘要中，我们发现一个 I - 样条参数非常小。为了避免在更复杂的模型中出现数值问题，我们将这个参数固定为零：

# 获取初始参数
init_params <- mmse_model$best
# 将第 7 个参数固定为 0
init_params[7] <- 0
# 重新定义 MMSE 模型
mmse_model_fixed <- lcmm(MMSE ~ I((age - 65) / 10), random = ~I((age - 65) / 10), subject = 'ID', data = paquid,
 link = "5 - quant - splines", verbose = FALSE, B = init_params, posfix = 7)
# 查看重新定义后的模型摘要
summary(mmse_model_fixed)

 

 

我们对 CESD 使用相同的规范。不过，需要注意的是，不同结果的规范可能会有所不同。

基于这两个子模型，我们可以估计多变量模型

由于两个结果在潜在类别条件下是独立的，并且我们只有一个类别，因此多变量模型等同于两个分离的子模型。我们可以检查估计值是否相同：

# 比较多变量模型和单变量模型的参数估计
cbind(multi_model_1$best, c(mmse_model_fixed$best, cesd_model$best))

 

估计多类别模型（G>1G>1）

对于多个类别，我们同样从单变量模型开始，然后估计多变量模型。这里不需要对单变量模型进行优化，因此我们使用 maxiter = 0 选项。我们还像之前一样将一个 I - 样条参数固定为零。

上述模型没有指定初始值，这些值将从单变量模型中提取。或者，我们可以使用单类别模型作为起点，或者进行网格搜索。

# 进行网格搜索估计双类别多变量模型
multi_m2_c <- gridsearch

我们接下来估计具有 3、4 和 5 个类别的相同模型：

 

我们为每个潜在类别数量保留最佳模型，并绘制一些用于选择最佳类别数量的统计标准：

在这个应用中，我们将选择 2 类别模型，因为它具有更好的 BIC。
我们可以使用桑基图来表示类别和个体流动之间的关系。目前，该包中还没有自动完成此操作的功能。

# 合并 5 个模型的分类结果
class_1_2 <- merge(multi_model_1$pprob[, 1:2], multi_model_2_b$pprob[, 1:2], by = "ID")
colnames(class_1_2) <- c("ID", "mm1", "mm2_b")
class_1_2_3 <- merge(class_

二分类模型的描述

模型总结

通过上述代码，得到二分类多变量模型的各项统计信息，包括数据集信息、模型收敛情况、拟合优度统计量以及各参数的估计值和检验结果等。从这些结果中可以了解到模型对数据的拟合效果，以及不同类别下各变量的固定效应和随机效应等情况。

 

 

更新单变量模型

函数返回用于指定模型的KK个单变量模型，并更新其输出。模型参数及其估计方差会被多变量框架中优化后的结果所替代。

# 更新二分类多变量模型对应的单变量模型
updated_models <- update(multi_model_2_b)
mmse_bivariate <- updated_models[[1]]
cesd_bivariate <- updated_models[[2]]

轨迹预测

只要数据框中包含模型指定的所有协变量，就可以为数据框中的任何数据计算类别特定的预测值。在接下来的代码中，通过生成一个年龄值在65到95之间的向量来创建这样一个数据框。使用predictY计算预测值，并使用相关的plot函数进行绘图。

# 绘制预测轨迹图
par(mfrow = c(1, 2))
plot(mmse_pred, lwd = 3, ylab = "MMSE", main = "MMSE预测轨迹", legend.loc = "bottomleft")
plot(cesd_pred, lwd = 3, ylab = "CESD", main = "CESD预测轨迹", legend = NULL)

若要逐点添加置信区间，可以在调用predictY时使用draws = TRUE选项。

模型评估

残差图

通过观察残差图，可以评估模型对数据的拟合程度。若残差呈现随机分布，则说明模型对数据的拟合效果较好；若残差存在明显的趋势或规律，则可能暗示模型存在改进的空间。

预测值与观测值对比图

为了评估所选模型的拟合效果，我们同时绘制每个潜在类别中观测值和预测值的对比图。

从对比图中，可以直观地看出模型预测值与实际观测值的接近程度，进一步判断模型的拟合优度。

分类

 

结果显示，“Class 1”由107名个体组成（占比21.4%），而393名个体属于“Class 2” 。同时，模型具有良好的区分能力，“Class 1”的平均后验概率为0.8280，“Class 2”的平均后验概率为0.8926。并且，有相当比例的个体后验概率高于0.7。

多变量联合模型

函数还允许使用比例风险模型对事件发生时间结果进行建模，并且可能涉及竞争风险的情况。

通过上述代码，实现了在潜在类别混合模型框架下，综合考虑纵向数据和事件发生时间数据的建模，进一步拓展了模型的应用范围，能够更全面地分析相关数据中的复杂关系。