全部代码都在这里了（有代码解释），大家可以自行取用。
# 轨迹分析# LCGAlibrary(lcmm)library(tidyverse)# devtools::install_github("hlennon/LCTMtools")library(LCTMtools)# help(LCTMtools)# 加载数据集data = read.table("ex8.1.dat",col.names = c("y1", "y2", "y3", "y4", 'x1','x2'))data %>% glimpse()data$id = seq(1,500,1) #增加个体编号# 预处理，宽型数据变长型数据newd = data %>% pivot_longer(cols = y1:y4,names_to = 'times',values_to = 'y') %>%select(times,y,id)newd %>% glimpse()# 建模set.seed(123)# 先建立一个基础模型，必须要有的mod1 = hlme(y~times,random=~times,subject = 'id',ng=1,data=newd)# 依次拟合有2~7个潜在类别的增长模型(可以自己写个循环)# 速度比较慢mod2 = hlme(y~times,random=~times,subject = 'id',mixture=~times,#当ng>1时需要设置ng=2,data=newd,B=mod1)mod3 = hlme(y~times,random=~times,subject = 'id',mixture=~times,#当ng>1时需要设置ng=3,data=newd,B=mod1)mod4 = hlme(y~times,random=~times,subject = 'id',mixture=~times,#当ng>1时需要设置ng=4,data=newd,B=mod1)mod5 = hlme(y~times,random=~times,subject = 'id',mixture=~times,#当ng>1时需要设置ng=5,data=newd,B=mod1,nproc=6)mod6 = hlme(y~times,random=~times,subject = 'id',mixture=~times,#当ng>1时需要设置ng=6,data=newd,B=mod1,nproc=6)mod7 = hlme(y~times,random=~times,subject = 'id',mixture=~times,#当ng>1时需要设置ng=7,data=newd,B=mod1,nproc=6)summary(mod1)summary(mod2)LCTMtoolkit(mod2)bic = rbind(mod1$BIC,mod2$BIC,mod3$BIC,mod4$BIC,mod5$BIC,mod6$BIC,mod6$BIC)bic# 上面的方法比较笨，我们可以用循环来实现，下面是个简单的演示# 可供大家参考，就不运行了for (i in 2:7) {modeli = hlme(y~times,random=~times,subject = 'id',mixture=~times,ng=i,data=newd,B=mod1,nproc=6)print(modeli$BIC)}# 上述代码很慢，所以我可能更倾向于用Mplus跑LCGA# 查看各模型的评价指标LCTMtoolkit(mod3)LCTMtoolkit(mod4)LCTMtoolkit(mod5)LCTMtoolkit(mod6)LCTMtoolkit(mod7)# 模型比较LCTMcompare(mod4, mod5)# 个体分配概率mod2$pprob# 各类别比例summarytable(mod2)# 数据出结果不易，保存下save.image()# 下次要用，可以直接导入模型load('.RData')
参考文献：
Lennon H, Kelly S, Sperrin M, Buchan I, Cross AJ, Leitzmann M, Cook MB, Renehan AG. Framework to construct and interpret latent class trajectory modelling. BMJ Open. 2018 Jul 7;8(7):e020683. doi: 10.1136/bmjopen-2017-020683. PMID: 29982203; PMCID: PMC6042544.
正文图片来自上述参考文献，封面是笔者自己画的思维导图。本文仅供学习、分享使用，如有侵权，请联系我们删除，谢谢。
