lcmm(2)-如何使用hlme函数估计潜在类别混合模型-2

教程地址：https://cecileproust-lima.github.io/lcmm/articles/latent_class_model_with_hlme.html

小心初始值！

用户预先指定的值

在以下示例中，初始值由用户预先指定：方差协方差的参数取线性混合模型的估计值，并为类别特定的轨迹尝试了任意的初始值。

m2b <- hlme(normMMSE ~ age65+I(age65^2)+CEP, random =~ age65+I(age65^2), subject = 'ID',data = paquid, ng = 2, mixture =~ age65+I(age65^2), B = c(0, 60, 40, 0, -4, 0, -10, 10, 212.869397, -216.421323,456.229910, 55.713775, -145.715516, 59.351000, 10.072221))

这个代码使用的是 hlme 函数，该函数通常用于纵向数据的混合效应模型。具体来说，它定义了一个具有潜在类别的混合效应线性模型 (Hidden Mixture Model)。以下是代码各部分的解释：

m2b <- hlme(normMMSE ~ age65+I(age65^2)+CEP, 
            random =~ age65+I(age65^2), 
            subject = 'ID',
            data = paquid, 
            ng = 2, 
            mixture =~ age65+I(age65^2), 
            B = c(0, 60, 40, 0, -4, 0, -10, 10, 212.869397, -216.421323,456.229910, 55.713775, -145.715516, 59.351000, 10.072221))

代码解读

1. normMMSE ~ age65 + I(age65^2) + CEP:

• 这是模型的固定效应部分，表示被解释变量 normMMSE（可能是一个标准化的认知评分）被解释变量 age65、age65^2（年龄的平方项）和 CEP 解释。
• I(age65^2) 用于在公式中指定 age65 的平方项，用于捕捉年龄的非线性影响。

• 指定随机效应结构。在这个模型中，age65 和 age65^2（年龄的平方项）有随机效应。这意味着每个受试者的这些变量的影响可以不同，从而捕捉个体的变化。

• 定义个体的标识符，这里使用 ID 列表示每个受试者的唯一标识符。这有助于识别每个个体的纵向数据。

• 使用的数据集为 paquid，这是一个包含变量 normMMSE、age65、CEP 和 ID 的数据框。

• 定义混合模型中的潜在类别（子组）的数量，这里设为 2。这意味着模型会尝试将数据划分为两个潜在类别（例如，两种不同的认知衰退模式）。

• 定义混合效应中类别特定的部分，表示 age65 和 age65^2 的效应可以在不同的潜在类别中有所不同。

• 开头的几个数字（如 0, 60, 40, 0, -4 等）通常代表固定效应的初始值。
• 后面的数值（如 212.869397, -216.421323 等）可能是随机效应和协方差的初始值，用于控制模型的复杂性和适应性。
• 这是模型的初始参数值向量。B 向量中的数值是参数的初始估计值，用于帮助模型收敛到最佳解。这些值通常需要根据经验或预估来设定，以提高模型拟合的效率和准确性。
• 具体数值含义：

值随机生成

另一种方法是根据单类别模型（这里是 m1）的估计值的渐近分布随机生成初始值：

m2c <- hlme(normMMSE ~ age65+I(age65^2)+CEP, random =~ age65+I(age65^2),subject = 'ID', data = paquid, ng = 2, mixture =~ age65+I(age65^2),  B = random(m1))

这段代码使用 hlme 函数定义了一个混合效应线性模型，具有两个潜在类别。和之前的代码类似，但不同之处在于，这里初始值的设定方式不同，具体使用了 B = random(m1)。以下是代码的详细解释：

m2c <- hlme(normMMSE ~ age65 + I(age65^2) + CEP, 
            random =~ age65 + I(age65^2), 
            subject = 'ID', 
            data = paquid, 
            ng = 2, 
            mixture =~ age65 + I(age65^2), 
            B = random(m1))

代码解读

1. normMMSE ~ age65 + I(age65^2) + CEP:

• 固定效应部分，表示 normMMSE（可能是认知测量值）由 age65、age65^2（年龄的平方项）和 CEP 解释。

• 随机效应结构。指定 age65 和 age65^2（年龄的平方项）具有随机效应，允许这些变量的效应在个体间有所不同。

• 设定个体标识符 ID，用于识别每个受试者的纵向数据。

• 使用的数据集是 paquid，包含分析所需的变量 normMMSE、age65、CEP 和 ID。

• 定义模型中的潜在类别（子群）数量为 2，表示模型会将数据分成两个潜在类别。

• 指定混合效应中类别特定的部分，age65 和 age65^2 的效应在不同的潜在类别中可以不同。

• 这里使用 random(m1) 生成初始值，m1 是一个单类别模型的对象。
• random(m1) 表示从 m1 的估计值的渐近分布中随机生成初始值，这是一种设置初始值的方式，可以增加模型拟合的稳定性。
• 通过从单类别模型的渐近分布中随机生成初始值，可以帮助模型更快收敛到最佳解，尤其在具有多个潜在类别时。

总结

这段代码和之前的 m2b 模型类似，都是用于拟合一个带有两个潜在类别的混合效应模型。但在 m2c 中，初始值是通过从单类别模型 m1 的估计值的渐近分布中随机生成的。这种方法能更好地探索参数空间，可能会使模型更容易收敛到全局最优解。

网格搜索

最后，可以使用 gridsearch 函数进行自动网格搜索。在下一个包含 ( G = 2 ) 和 ( G = 3 ) 类（分别为 m2d 和 m3b）的示例中，hlme 将从 100 个随机初始值向量开始，最多运行 30 次迭代。然后，仅对在 30 次迭代后提供最佳对数似然值的起始点完成估计过程。

m2d <- gridsearch(hlme(normMMSE ~ age65+I(age65^2)+CEP,  random =~ age65+I(age65^2), subject = 'ID', data=paquid, ng = 2, mixture=~age65+I(age65^2)), rep=100, maxiter=30, minit=m1)

这段代码使用了 gridsearch 函数来自动执行网格搜索，以找到模型的最佳初始值。具体来说，它在混合效应线性模型中使用了多个随机初始值来寻找最优解。以下是代码的详细解释：

m2d <- gridsearch(hlme(normMMSE ~ age65 + I(age65^2) + CEP, 
                       random =~ age65 + I(age65^2), 
                       subject = 'ID', 
                       data = paquid, 
                       ng = 2, 
                       mixture =~ age65 + I(age65^2)), 
                 rep = 100, 
                 maxiter = 30, 
                 minit = m1)

代码解读

1. hlme(normMMSE ~ age65 + I(age65^2) + CEP, random =~ age65 + I(age65^2), subject = 'ID', data = paquid, ng = 2, mixture =~ age65 + I(age65^2)):

• normMMSE ~ age65 + I(age65^2) + CEP：固定效应部分，解释变量包括 age65、age65^2 和 CEP。
• random =~ age65 + I(age65^2)：随机效应结构，表示 age65 和 age65^2 具有随机效应。
• subject = 'ID'：用于指定每个个体的唯一标识符。
• data = paquid：数据集为 paquid。
• ng = 2：定义两个潜在类别。
• mixture =~ age65 + I(age65^2)：指定类别特定的效应，允许 age65 和 age65^2 在不同类别中有不同的影响。
• 这部分代码定义了一个混合效应线性模型 (hlme)，其中包含固定效应、随机效应、个体标识符、数据集、类别数量以及类别特定效应。
• 解释如下：

• gridsearch 函数用于执行网格搜索，以便找到最佳的初始值，从而帮助模型更好地拟合。
• rep = 100：指定要运行 100 组不同的随机初始值。
• maxiter = 30：每组初始值最多运行 30 次迭代。
• minit = m1：表示使用模型 m1 的结果作为初始值的一个起点。m1 是一个先前拟合的模型，可能是一个单类别模型，用来为多类别模型提供稳定的初始参数。

总结

这段代码通过 gridsearch 执行网格搜索，拟合一个具有两个潜在类别的混合效应模型。代码会从 100 个随机生成的初始值向量中进行迭代，每个向量最多运行 30 次迭代。最后，hlme 函数会选择在这些迭代中对数似然值最佳的初始值，从而优化模型的拟合效果。这种方法可以帮助避免因初始值选择不当而导致模型收敛到局部最优解的问题。

m3g <- gridsearch(hlme(normMMSE ~ age65+I(age65^2)+CEP,  random =~ age65+I(age65^2), subject = 'ID', data=paquid, ng = 3, mixture=~age65+I(age65^2)), rep=100, maxiter=30, minit=m1)

这段代码使用了 gridsearch 函数来执行网格搜索，以拟合一个具有三个潜在类别的混合效应线性模型。代码的结构和之前的 m2d 类似，但这里将类别数量设置为 3。以下是代码的详细解释：

m3g <- gridsearch(hlme(normMMSE ~ age65 + I(age65^2) + CEP, 
                       random =~ age65 + I(age65^2), 
                       subject = 'ID', 
                       data = paquid, 
                       ng = 3, 
                       mixture =~ age65 + I(age65^2)), 
                 rep = 100, 
                 maxiter = 30, 
                 minit = m1)

代码解读

1. hlme(normMMSE ~ age65 + I(age65^2) + CEP, random =~ age65 + I(age65^2), subject = 'ID', data = paquid, ng = 3, mixture =~ age65 + I(age65^2)):

• normMMSE ~ age65 + I(age65^2) + CEP：固定效应部分，表示 normMMSE（可能是一个认知测量值）由 age65、age65^2（年龄的平方项）和 CEP 解释。
• random =~ age65 + I(age65^2)：指定随机效应结构，表示 age65 和 age65^2 的效应在个体间可以有所不同。
• subject = 'ID'：定义个体标识符 ID，用于识别每个受试者的纵向数据。
• data = paquid：指定数据集为 paquid。
• ng = 3：定义潜在类别数量为 3，表示模型会将数据划分为三个潜在类别。
• mixture =~ age65 + I(age65^2)：指定类别特定效应部分，允许 age65 和 age65^2 的效应在不同类别中有所不同。
• 这部分代码定义了一个混合效应线性模型 (hlme)，包含固定效应、随机效应、个体标识符、数据集、类别数量和类别特定效应。
• 具体解释如下：

• gridsearch 函数执行网格搜索，以便找到模型的最佳初始值。
• rep = 100：指定运行 100 组不同的随机初始值。
• maxiter = 30：每组初始值最多运行 30 次迭代。
• minit = m1：表示使用模型 m1 的结果作为初始值的一个起点。m1 是一个先前拟合的模型，可能是一个单类别模型，用来为多类别模型提供稳定的初始参数。

选择最佳模型

一组模型的估计过程（通常具有不同数量的潜在类别）可以使用 summarytable 进行汇总，并使用 summaryplot 绘制出来。可以显示不同的信息：

summarytable(m1,m2,m2b,m2c, m2d  , m3g, which = c("G", "loglik", "conv", "npm", "AIC", "BIC", "SABIC", "entropy","ICL", "%class"))
    G    loglik conv npm      AIC      BIC    SABIC   entropy     ICL1     ICL2
m1  1 -8920.623    1  11 17863.25 17909.61 17874.69 1.0000000 17909.61 17909.61
m2  2 -8899.228    1  15 17828.46 17891.67 17844.06 0.5504369 18047.48 18050.54
m2b 2 -8899.228    1  15 17828.46 17891.67 17844.06 0.5504366 18047.48 18050.54
m2c 2 -8899.228    1  15 17828.46 17891.67 17844.06 0.5504369 18047.48 18050.54
m2d 2 -8899.228    1  15 17828.46 17891.67 17844.06 0.5504370 18047.48 18050.54
m3g 3 -8891.351    1  19 17820.70 17900.78 17840.47 0.6240325 18107.30 18089.55
    %class1 %class2 %class3
m1    100.0                
m2     12.4    87.6        
m2b    87.6    12.4        
m2c    12.4    87.6        
m2d    12.4    87.6        
m3g    85.8     4.0    10.2

summaryplot(m1, m2, m3g, which = c("BIC", "entropy","ICL"))

AI解释

这段代码使用 summarytable 和 summaryplot 函数来汇总和可视化一组混合效应模型的估计结果，以帮助选择最佳模型。以下是代码和结果的详细解释。

代码解释

summarytable(m1, m2, m2b, m2c, m2d, m3g, 
             which = c("G", "loglik", "conv", "npm", "AIC", "BIC", "SABIC", "entropy", "ICL", "%class"))

• summarytable 函数：用于生成多个模型的汇总表，以便对不同模型的性能指标进行比较。
• 模型列表：m1、m2、m2b、m2c、m2d 和 m3g 是要比较的模型。
• m1 是单类别模型。
• m2、m2b、m2c 和 m2d 是具有两个潜在类别的模型，使用了不同的初始值设定方式或网格搜索方法。
• m3g 是具有三个潜在类别的模型。
• which 参数：指定要显示的模型统计信息指标。
• G：类别数。
• loglik：对数似然值，值越大越好。
• conv：收敛指示符，1 表示模型成功收敛。
• npm：非零参数数。
• AIC、BIC、SABIC：信息准则，数值越小越好。
• entropy：熵，表示类别划分的清晰度，越接近 1 表示分类越清晰。
• ICL：分类似然信息准则，值越小越好。
• %class：各类别的样本百分比分布。

输出解释

• m1：单类别模型，熵为 1，因为所有数据都属于一个类别。
• m2、m2b、m2c、m2d：双类别模型。所有模型的对数似然值、AIC、BIC、SABIC 和 ICL 指标完全一致，说明它们收敛到了相同的解，仅类别顺序不同（这是一种“标签交换”现象）。熵约为 0.55，表示分类清晰度较低。两个类别的样本分布在 12.4% 和 87.6% 之间。
• m3g：三类别模型，AIC、BIC 和 SABIC 数值略低，可能是更好的选择。熵为 0.62，略高于双类别模型，表明分类更清晰。三类别的分布分别为 85.8%、4.0% 和 10.2%。

总结

• 模型选择：通常根据 BIC、AIC、SABIC 和 ICL 选择最佳模型，值越小越好。根据 BIC 值，m2 是最佳选择，而 AIC 和修正后的 BIC 则偏向于 m3g。
• 类别清晰度：熵值显示三类别模型（m3g）具有更好的辨别力。然而，三类别模型包含一个非常小的类别（4.0%），可能不具实际意义。

可视化代码解释

summaryplot(m1, m2, m3g, which = c("BIC", "entropy", "ICL"))

• **summaryplot**：绘制模型 m1、m2 和 m3g 的图表，展示指定的指标。
• which 参数：显示 BIC、熵和 ICL 指标，以便直观比较这些模型。

通过图表，可以更清晰地看到各模型在 BIC、熵和 ICL 指标上的表现，有助于确定最佳模型。

我们在此总结之前估计的 6 个模型。可以看到，所有的双类别模型都收敛到了同一个估计点，只有类别的顺序不同。这种现象称为“标签交换”，常常会发生。这并不是一个问题，因为它们达到了相同的最优解。

这个示例说明了确定“最佳潜在类别数”的复杂性。实际上，根据推荐的 BIC 值，应该选择双类别模型（因为它提供了最低的 BIC 值）。但是，AIC 和大小调整后的 BIC（它们带有较小的惩罚）都更偏向三类别模型。熵值也偏向三类别模型，认为其具有更好的辨别能力（熵更接近 1）。然而，三类别模型包含一个非常小的类别，这通常不是我们所寻找或感兴趣的异质性。在这个例子中，取决于统计和临床标准，可以选择双类别模型或三类别模型。下文中，我们选择了双类别模型用于最终结果描述。

AI解释-标签交换（label switching）

在上述分析中，标签交换（label switching） 现象发生的原因和意义主要涉及模型中类别的本质和模型的解释方法：

原因

1. 类别的对称性：在潜在类别模型中，不同类别在数值计算上没有固定的顺序。类别 1 和类别 2 只是标签，在模型的统计结构上是等价的。模型没有对类别赋予任何先验顺序，因此在多个模型中，类别可能会被不同地标记。

2. 相同的最优解：所有双类别模型（m2、m2b、m2c、m2d）在同一个估计点收敛，说明它们实际上达到了相同的最优解。尽管类别的标签可能不同（例如，一个模型中类别 1 代表高概率群体，类别 2 代表低概率群体，而另一个模型正好相反），但模型的整体结构和拟合效果没有变化。

3. 初始化的随机性：模型可能使用不同的初始值或搜索方法（如网格搜索）来进行拟合。尽管初始化不同，模型收敛到了同一个解，但由于计算过程的细微差别，标签可能在不同模型中切换。

意义

1. 对模型质量没有影响：标签交换不影响模型的对数似然值、AIC、BIC、SABIC 或 ICL 等指标。这些统计量关注的是模型整体的拟合优度，而不是类别标签本身。因此，标签交换不改变模型的估计质量和对数据的解释能力。

2. 分析和解释时需要一致性：在选择最终模型并报告结果时，确保类别标签的一致性可以避免混淆。例如，如果选择了 m2d 作为最终模型，可以固定类别 1 和类别 2 的含义，以便对类别特征进行描述和解读。

3. 说明了类别的不确定性：标签交换表明了类别的相对性。由于潜在类别模型是无监督的，类别标签本身并没有实质性意义，它们只是对数据特征的某种划分。这种现象提醒我们在解释潜在类别时，应该关注类别的实际特征（如各类别的均值、标准差等），而非类别标签本身。

总结

标签交换是一种常见的现象，主要因为类别标签的相对性和模型的等价性。它不会影响模型的拟合效果，但在解释和报告结果时，需要选择一种类别标签的固定方案，以保持一致性。标签交换的出现也提醒我们，潜在类别的划分是基于统计结构的，类别标签并不具有绝对意义。

最新课程-基于R语言的动态预测模型课程-胖子老师独自授课

开课目的及前言

预测模型作为真实世界研究的重要组成部分，其研究被广泛开展。但是，传统的预测模型利用基线数据对最终的生存结果进行预测，这种模型无法纳入患者在后续随访中可能会动态变化的重要数据（比如肿瘤标记物的动态变化）。 以上情况在统计学中会产生估计偏差情况，也是不符合临床实际的。近年来发展起来的动态预测模型方法，利用患者的多次随访数据，结合患者的基线数据，对最终患者的额生存结果（或类似的time to event事件）进行估计。其发文量呈现快速增长趋势。

在临床实际中，医生会根据患者的动态变化指标做出进一步诊断及治疗的判断。动态预测模型结合患者的纵向数据与最终的生存结果，对于最终结果进行更加准备的预测。由于当前R语言在医学统计工作中占据重要地位，但很多临床大夫、护士因为时间工作关系很难将R语言与临床科研相结合，故开设R语言动态预测模型课程，旨在快速让学员掌握统计工作中常用到的R语言，助力临床科研工作。天企助力(天津)生产力促进有限公司特举办“基于R语言的动态预测模型课程培训班”。

预测模型类文章目前总结起来发展经历了以下三个阶段：

1. 基于传统流行病学的列线图模型（本质都是cox回归及glm回归），简单的统计学分析模型，是模型依赖的方法，临床上实际情况很难满足其前提假设，实际效果不好。

2. 基于机器学习/深度学习的预测模型的构建（在数据上提高了维度，在算法上引入了机器学习），虽然算法上引入了机器学习模型，处理数据更加灵活，模型的假设也更少。但是在使用的数据上还是患者的一次基线数据进行预测，与临床实际不符。

3. 基于纵向数据的动态预测模型（基于纵向多次随访数据，模型应用联合模型等动态预测模型方法），应用患者的多次随访数据对最终的生存结果进行预测，从数据和方法上都更类似于临床实际。

考虑到动态预测模型有以下特点，因此必然是后续高分文章的必备方法：

1. 数据上必须有同一个患者的多次随访数据，相对于既往横断面一次基线数据，数据的收集难度更大，而且动态预测模型需拟合纵向的线性混合模型，因此需要的数据量较大。这就提示我们如果能收集到如上数据更加容易发高分文章。

2. 应用方法学动态预测模型需首先掌握普通生存分析及普通预测模型的方法，并且还需要熟悉纵向数据分析的广义线性混合模型，再次基础上还需要掌握tidyverse语法基础来将自己的数据转换为满足函数要求的纵向数据，另外对于联合模型，模型的结合形式及变量选择也均需要从临床背景及统计学方法考虑。

近期高分文章举例

授课老师

灵活胖子-独自

双一流学校肿瘤学博士毕业，目前就职于国内五大肿瘤中心之一。科研方向为真实世界研究，生物信息学分析及人工智能研究。目前以第一或共同第一作者身份发表SCI论文10余篇，累计IF50+。目前与国内多个院校及医院有科研合作。联合翻译小组同学，在国内第一次将jmbayes2及dynamicLM全文翻译为中文并在公众号发表。

课程目录及安排

授课形式及时间

授课形式：远程在线实时直播授课。

授课时间：2024年12月开课，总课时不少于30小时，每周进行3-5小时的授课，有充分时间学习，预计6-8周完成所有授课内容。

答疑支持：建立课程专属微信群，1年内课程内容免费答疑。

视频回看：3年内免费无限次回看。

课程售价及售后保证

课程售价：总价3000元，报名可先交300元预定即可，开课后2周内交齐即可

对公转账等手续务必提前联系助教

承办公司：天企助力（天津）生产力促进有限公司

奖励政策：学员应用所学内容发表IF 10+文章可退还学费（具体要求及流程需要咨询助教）

报名咨询

可联系我的助教进行咨询

助教联系电话：18502623993

正式通知

pdf版通知可联系助教获取