R数据分析：潜类别轨迹模型LCTM的做法，实例解析
最近看了好多潜类别轨迹latent class trajectory models的文章，发现这个方法和我之前常用的横断面数据的潜类别和潜剖面分析完全不是一个东西，做纵向轨迹的正宗流派还是这个方法，当然了这个方法和潜增长和增长曲线模型在做法并没有实际区别，都是用的hlme这个函数。但是文献中的叫法和花样就比较多了。
像本文写的latent class trajectory models，之前写的潜类别增长模型LCGA和增长曲线模型GMM都是潜类别线性混合模型latent class linear mixed models (LCLMM)的分支。
The major difference between LCGA and GMM is that LCGA does not allow within-class variation whereas GMM does allow within-class variation
像这一类的模型都是用hlme这个函数跑，这篇文章也可以看作是作为之前的潜增长和增长曲线文章的一个实际应用的延续。
应用背景
很多的同学关心某个变量的纵向发展轨迹，并且还感兴趣不同轨迹对某个结局的影响如何。如果你的研究也涉及到这样的问题，你就可以考虑用潜类别轨迹模型了，参考文献也甩给大家，大家感兴趣可以去瞅瞅下面这个文章：
Mirza, S. S., Wolters, F. J., Swanson, S. A., Koudstaal, P. J., Hofman, A., Tiemeier, H., & Ikram, M. A. (2016). 10-year trajectories of depressive symptoms and risk of dementia: a population-based study. The Lancet Psychiatry, 3(7), 628-635.
文章作者通过潜类别轨迹模型将人群抑郁症状发展轨迹分成了5类，最终发现只有特定类轨迹才和随后的痴呆有关系，这对痴呆的干预和抑郁痴呆的关系的认识都是有重要意义的。
今天我就仿照这篇文章给大家写写如何做潜类别轨迹模型。
潜类别轨迹的报告内容
做之前我们还是看看这篇文献中是如何介绍这个方法的。
We used latent class trajectory models to identify trajectories of depressive symptoms over time. This is a specialised form of finite mixture modelling, and is designed to identify latent classes of individuals following similar progressions of a determinant over time or with age.
可以看到这个方法在重要作用是识别那些随着时间或者年龄拥有相似症状和疾病进程的人群类别。比如做抑郁的潜类别轨迹就是要识别出人群中可能的抑郁进程亚组。文章中也说明了这个模型就是一个特殊的混合模型specialised form of finite mixture modelling。
在拟合症状随着时间或者年龄变化的时候，我们允许或者说我们需要去考虑症状和时间的曲线关系的，就是说不能简单第认为某个症状的纵向变化一定是线性的，意思就是我们要考虑时间变量的高次项，一般来讲二次就够了。作者的论文中也是加了时间的二次项的。然后根据BIC确定最优类别数，同时确保后验概率大于0.7，类别人数大于0.02.
最终作者出图如下：
作者根据这个轨迹走势，还给每个类别进行了命名，打上标签，包括Low symptoms，Decreasing symptoms，Remitting symptoms，Increasing symptoms， High symptoms，然后将轨迹标签作为预测变量进行了后续的生存分析。
对于轨迹部分结果的报告，因为这个文章轨迹只是一部分而并非主要目的，所以报告很少,只有每一个轨迹的人群数量和占比。
接下来给大家分享如何做这么一个潜类别轨迹。
潜类别轨迹的做法
潜类别轨迹有专门的R包可以做，感兴趣的同学可以去看这篇文章：
Lennon H, Kelly S, Sperrin M, et al Framework to construct and interpret latent class trajectory modelling BMJ Open 2018;8:e020683. doi: 10.1136/bmjopen-2017-020683
上面这篇文章给出了潜类别轨迹模型的做法框架，共8步：
本文的绝大部分步骤也都是参考的上面的文章。
现在我手上的数据长这样：
一个纵向的长型数据，包括每个人不同年龄段测得的bmi，我现在就想看看随着年龄的增长，人群bmi轨迹是不是存在异质性亚组。接下来我就用潜类别轨迹模型回答这个问题，并且出图，并得到每个轨迹类别的人数和占比。
首先我们写出轨迹类别数量为1时候的潜类别轨迹的代码：
m.1 <- hlme(bmi ~ 1+ age + I(age^2),
random = ~ 1 + age,
ng = 1,
data = data.frame(bmi), subject = "id")
关于hlme之前给大家写过各个参数的意思，上面的代码就是要拟合bmi随着年龄变化的轨迹，同时考虑年龄的随机效应（截距+斜率），并声明嵌套的高水平subject = "id"。运行上面的代码我们就拟合了一个轨迹类别为1的模型。
之后我们还需要拟合2到7个类别的模型，这个7是上面文献推荐的哈，我们可以写个循环语句，一次搞定（为什么不从模型1循环到7呢？是因为ng参数为1时我们并不需要设定mixture参数，所以2到7写了循环，1单独做）：
lin <- c(m.1$ng, m.1$BIC)
for (i in 2:7) {
mi <- hlme(fixed = bmi ~ 1+ age + I(age^2),
mixture = ~ 1 + age + I(age^2),
random = ~ 1 + age,
ng = i, nwg = TRUE,
data = data.frame(bmi), subject = "id")
lin <- rbind(lin, c(i, mi$BIC))
}
7个模型跑完，我们需要对比每个模型的BIC（这个也是参考的The Lancet Psychiatry那篇文章的做法），所以我们对模型和相应的BIC进行展示，然后以BIC最小为标准确定轨迹数量：
从上图就可以看得出我们轨迹数量确定为5个时，模型的BIC最小，由此可以确定轨迹数量为5，就是说人群bmi发展轨迹存在5个异质性亚组。
按照论文报告的要求我们需要出图，根据图中每条轨迹的走势确定轨迹类别标签，还有每个轨迹类别的人群数量和占比，具体方法如下：
首先，进行图形的绘制，我们解决这类问题（包括机器学习模型）的基本思路依然是通过自我数据得到模型，通过模型拟合新数据出图，代码如下：
plotpred <- predict(m5, datnew, var.time ="age", draws = TRUE)
plot(plotpred, lty=2,xlab="Age", ylab="BMI", legend.loc = "topleft", cex=0.75)
上面的代码中m5为我们拟合的5个轨迹类别的模型对象。
运行代码就可以得到5个轨迹类别模型的轨迹图如下：
然后我们就可以根据图的走势给每个潜类别轨迹打上有临床意义的标签了，有了标签就可以进行后续的建模了。
然后我们还需要报告每个轨迹类别的人数和占比，方法如下：
m5$pprob
运行代码即可得到，每个个案到底属于哪个轨迹类别，以及其属于每个轨迹类别的概率，如下图：
到这儿我们的模型基本跑完了，论文中有提到后验概率大于0.7，类别人数大于0.02，这个在模型总结中也是可以调出来的：
基本上掌握了上面的方法，柳叶刀精神病学那篇文章的前半部分统计分析就完成了，如果你想将某个症状或者疾病发展的轨迹作为自变量的相关研究都可以进行了，好了今天要给大家分享的潜类别轨迹模型的做法就是这样。
后记：写完这篇文章，我越来越觉得潜类别轨迹模型和增长混合模型就是一个东西，只不过不同的学者用词不一样，不知各位看官怎么看，可以私信我交流。
另外，还特别建议大家好好去看lcmm包的说明文档，相信大家看完之后还会有更大收获。