【原】R数据分析：论文中的轨迹的做法，潜增长模型和增长混合模型

CodewarCodewar 2021-08-21

展开全文

好多同学手上有纵向数据，想看轨迹，看人群异质性，咨询做法，今天给大家写两个方法，一个叫潜增长模型Latent Class Growth Analyses (LCGA) ，一个叫增长混合模型Growth Mixture Modeling (GMM)。

这个异质性怎么看呢，就是基于人的不同发展的轨迹的出来的：

They can be used to identify latent subgroups, classes or clusters of individuals based on their common growth trajectories over time.

两个模型都可以看作是增长模型growth mode的拓展：

LCGA can roughly be seen as an extension of a fixed effect growth model, whereas GMM can be seen as an extension of a random effect growth model

LCGA和GMM的区别

这两种方法都是将传统增长模型与潜类别分析相结合的模型，既可以刻画增长趋势又可以考虑群体异质性的目的。

二者的区别主要在于类别组内的发展轨迹是否考虑增长曲线内部的个体扰动，潜增长模型可以看作是增长混合模型的特例，就是说增长混合模型不考虑随机效应的时候就可以认为是潜增长模型：

A special case of GMMs is latent class growth analysis (LCGA)[15],[16] which does not allow for departure from the average trajectory within each latent class。Thus, in contrast to mixed effects models where each subject's intercept and slope are drawn from a normal distribution or GCMMs where they are drawn from a mixture of normal distributions,LCGAs allow only for a limited set of discrete options。

实例操练

我现在手上有500个重复测量的数据集，100个观测，每个观测连续测量5次。数据大概长这样：

这个数据还有一个协变量covar，是一个取0和1的二分类数据，对于这么一个纵向数据，我可以先非常直观地把每个个体的轨迹画出来，假设它是这样的：

其实画出来大体一看，似乎是有两个类别的轨迹出现的，具体是不是呢？我们得使用分析方法验证。

我们的分析的目的就是识别这些轨迹的异质性，从而将人群划分为不同的类别。

先看用潜增长模型如何做，我们需要用到lcmm包中的hlme函数，基本形式如下图：

其中fixed为线性混合模型的固定效应部分，“~”符号左边写因变量右边写自变量，自变量用加号链接。mixture参数只有在类别数大于一的时候才需要设置，我们做1个类别是不需要的；random参数是随机效应部分，因为我们做的是潜增长（没有混合），也没有必要设置这一个参数，subject用来设置嵌套结构的主体，此例中是“ID”；ng是潜类别个数；classmb是逻辑增长中的协变量，所以也不需要设置；最终我们写出代码如下：

lcga1 <-hlme(y ~ time, subject = "ID", ng = 1, data = mydata) 
lcga2 <-gridsearch(rep = 100, maxiter = 10, minit = lcga1,
                   hlme(y ~ time, subject = "ID",
                        ng = 2, data = mydata, mixture = ~ time)) 
lcga3 <-gridsearch(rep = 100, maxiter = 10, minit = lcga1,
                   hlme(y ~ time, subject = "ID",
                        ng = 3, data = mydata, mixture = ~ time))

在上面的代码中因为我们是做潜增长，所以省去了随机效应部分（slope and intercept)，什么意思呢，就是说我们的模型中做出来的每一类都只考虑固定效应，不会考虑每一类中的个体变异了。还有需要注意的是我们跑一类之后的类别的时候是在gridsearch这个函数之中去嵌套了一hlme，这个操作是为了获得全局最优的结果，具体原理就是将每一个hlme函数用不同的起始值跑100遍。此时我们还设置了mixture参数，因为我们是跑大于1个类别了嘛，就是说每一类我都要考虑时间的固定效应。

lcga3的代码的解释也请参考上段。

运行上面的代码之后我们的潜增长模型就跑好了，下面的代码可以方便地比较不同类别数量模型的拟合优度，从而帮助我们判断：

summarytable(lcga1, lcga2, lcga3)

我们还可以用summary很方便地查看具体模型的信息：

summary(lcga2)

运行代码便可以得到我们需要在论文中报告的系数了：

到此，潜增长模型做完。

继续，增长混合模型，增长混合模型的混合又分为两种了，一种是随机截距，另一种是随机斜率，我们分开看

先看随机截距

还是我们之前的数据，我们做随机截距增长混合模型可以写出如下代码：

gmm1 <-hlme(y ~ time, subject = "ID", random=~1, ng = 1, data = mydata)
gmm2 <-gridsearch(rep = 100, maxiter = 10, minit = gmm1, 
                  hlme(y ~ time, subject = "ID", random=~1,
                       ng = 2, data = mydata, mixture = ~ time, nwg=T))
gmm3 <-gridsearch(rep = 100, maxiter = 10, minit = gmm1,
                  hlme(y ~ time, subject = "ID", random=~1,ng = 3, 
                       data = mydata, mixture = ~ time, nwg=T))

可以看到增长混合模型与增长模型唯一的不同就是混合模型多了一个random参数

这个random参数就是用来设定随机效应的，我们只要随机截距所以直接设定为1就行。上面的代码中还有一个参数nwg我们设定为True，意思是随机效应的方差协方差是类别特异的，我们的例子中就是说每个类别的随机截距的方差是不同的：

运行上面的代码，一个带随机截距的增长混合模型就拟合好了，我们看结果：

summarytable(gmm1, gmm2, gmm3)

可以看到随机截距增长混合模型输出和和潜增长模型的差别就在于多了一个随机效应的方差协方差矩阵，在我们的结果中，类别2的截距方差为0.306，类别1的截距方差为0.306*1.12=0.343

以上就是随机截距增长混合模型。

接着看随机斜率增长混合模型

在随机斜率增长混合模型中我们认为，每一个类别中每个人允许有不同的时间效应，就是说每个人的增长斜率可以不一样，具体我们写出如下代码：

gmm1_2 <-hlme(y ~ time, subject = "ID", random=~1 + time, ng = 1, 
              data =mydata)
gmm2_2 <-gridsearch(rep = 100, maxiter = 10, minit = gmm1_2, 
                    hlme(y ~ time, subject = "ID", random=~1 + time,
                         ng = 2, data = mydata, mixture = ~ time, nwg=T))
gmm3_2 <-gridsearch(rep = 100, maxiter = 10, minit = gmm1_2,
                    hlme(y ~ time, subject = "ID", random=~1+time,ng = 3, 
                         data = mydata, mixture = ~ time, nwg=T))