好多同學(xué)手上有縱向數(shù)據(jù)，想看軌跡洽沟，看人群異質(zhì)性，咨詢做法蜗细，今天給大家寫兩個(gè)方法裆操，一個(gè)叫潛增長(zhǎng)模型Latent Class Growth Analyses (LCGA) ，一個(gè)叫增長(zhǎng)混合模型Growth Mixture Modeling (GMM)炉媒。

這個(gè)異質(zhì)性怎么看呢踪区，就是基于人的不同發(fā)展的軌跡的出來的：

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

兩個(gè)模型都可以看作是增長(zhǎng)模型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的區(qū)別

這兩種方法都是將傳統(tǒng)增長(zhǎng)模型與潛類別分析相結(jié)合的模型，既可以刻畫增長(zhǎng)趨勢(shì)又可以考慮群體異質(zhì)性的目的吊骤。

二者的區(qū)別主要在于類別組內(nèi)的發(fā)展軌跡是否考慮增長(zhǎng)曲線內(nèi)部的個(gè)體擾動(dòng)缎岗，潛增長(zhǎng)模型可以看作是增長(zhǎng)混合模型的特例，就是說增長(zhǎng)混合模型不考慮隨機(jī)效應(yīng)的時(shí)候就可以認(rèn)為是潛增長(zhǎng)模型：

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传泊。

實(shí)例操練

我現(xiàn)在手上有500個(gè)重復(fù)測(cè)量的數(shù)據(jù)集，100個(gè)觀測(cè)鸭巴，每個(gè)觀測(cè)連續(xù)測(cè)量5次眷细。數(shù)據(jù)大概長(zhǎng)這樣：

這個(gè)數(shù)據(jù)還有一個(gè)協(xié)變量covar，是一個(gè)取0和1的二分類數(shù)據(jù)鹃祖，對(duì)于這么一個(gè)縱向數(shù)據(jù)溪椎，我可以先非常直觀地把每個(gè)個(gè)體的軌跡畫出來，假設(shè)它是這樣的：

其實(shí)畫出來大體一看恬口，似乎是有兩個(gè)類別的軌跡出現(xiàn)的校读，具體是不是呢？我們得使用分析方法驗(yàn)證祖能。

我們的分析的目的就是識(shí)別這些軌跡的異質(zhì)性歉秫，從而將人群劃分為不同的類別。

先看用潛增長(zhǎng)模型如何做芯杀，我們需要用到lcmm包中的hlme函數(shù)端考，基本形式如下圖：

其中fixed為線性混合模型的固定效應(yīng)部分雅潭，“~”符號(hào)左邊寫因變量右邊寫自變量，自變量用加號(hào)鏈接却特。mixture參數(shù)只有在類別數(shù)大于一的時(shí)候才需要設(shè)置扶供，我們做1個(gè)類別是不需要的；random參數(shù)是隨機(jī)效應(yīng)部分裂明，因?yàn)槲覀冏龅氖菨撛鲩L(zhǎng)（沒有混合）椿浓，也沒有必要設(shè)置這一個(gè)參數(shù)，subject用來設(shè)置嵌套結(jié)構(gòu)的主體闽晦，此例中是“ID”扳碍；ng是潛類別個(gè)數(shù)；classmb是邏輯增長(zhǎng)中的協(xié)變量仙蛉，所以也不需要設(shè)置笋敞；最終我們寫出代碼如下：

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))

在上面的代碼中因?yàn)槲覀兪亲鰸撛鲩L(zhǎng)，所以省去了隨機(jī)效應(yīng)部分（slope and intercept)荠瘪，什么意思呢夯巷，就是說我們的模型中做出來的每一類都只考慮固定效應(yīng)，不會(huì)考慮每一類中的個(gè)體變異了哀墓。還有需要注意的是我們跑一類之后的類別的時(shí)候是在gridsearch這個(gè)函數(shù)之中去嵌套了一hlme趁餐，這個(gè)操作是為了獲得全局最優(yōu)的結(jié)果，具體原理就是將每一個(gè)hlme函數(shù)用不同的起始值跑100遍篮绰。此時(shí)我們還設(shè)置了mixture參數(shù)后雷，因?yàn)槲覀兪桥艽笥?個(gè)類別了嘛，就是說每一類我都要考慮時(shí)間的固定效應(yīng)吠各。

lcga3的代碼的解釋也請(qǐng)參考上段臀突。

運(yùn)行上面的代碼之后我們的潛增長(zhǎng)模型就跑好了，下面的代碼可以方便地比較不同類別數(shù)量模型的擬合優(yōu)度走孽，從而幫助我們判斷：

summarytable(lcga1, lcga2, lcga3)

我們還可以用summary很方便地查看具體模型的信息：

summary(lcga2)

運(yùn)行代碼便可以得到我們需要在論文中報(bào)告的系數(shù)了：

到此惧辈，潛增長(zhǎng)模型做完。

繼續(xù)磕瓷，增長(zhǎng)混合模型盒齿，增長(zhǎng)混合模型的混合又分為兩種了，一種是隨機(jī)截距困食，另一種是隨機(jī)斜率边翁，我們分開看

先看隨機(jī)截距

還是我們之前的數(shù)據(jù)，我們做隨機(jī)截距增長(zhǎng)混合模型可以寫出如下代碼：

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))

可以看到增長(zhǎng)混合模型與增長(zhǎng)模型唯一的不同就是混合模型多了一個(gè)random參數(shù)

這個(gè)random參數(shù)就是用來設(shè)定隨機(jī)效應(yīng)的硕盹，我們只要隨機(jī)截距所以直接設(shè)定為1就行符匾。上面的代碼中還有一個(gè)參數(shù)nwg我們?cè)O(shè)定為True，意思是隨機(jī)效應(yīng)的方差協(xié)方差是類別特異的瘩例，我們的例子中就是說每個(gè)類別的隨機(jī)截距的方差是不同的：

運(yùn)行上面的代碼啊胶，一個(gè)帶隨機(jī)截距的增長(zhǎng)混合模型就擬合好了甸各，我們看結(jié)果：

summarytable(gmm1, gmm2, gmm3)

可以看到隨機(jī)截距增長(zhǎng)混合模型輸出和和潛增長(zhǎng)模型的差別就在于多了一個(gè)隨機(jī)效應(yīng)的方差協(xié)方差矩陣，在我們的結(jié)果中焰坪，類別2的截距方差為0.306趣倾，類別1的截距方差為0.306*1.12=0.343

以上就是隨機(jī)截距增長(zhǎng)混合模型。

接著看隨機(jī)斜率增長(zhǎng)混合模型

在隨機(jī)斜率增長(zhǎng)混合模型中我們認(rèn)為某饰，每一個(gè)類別中每個(gè)人允許有不同的時(shí)間效應(yīng)儒恋，就是說每個(gè)人的增長(zhǎng)斜率可以不一樣，具體我們寫出如下代碼：

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))

可以看到黔漂，代碼的不同之處就是random參數(shù)的設(shè)定不一樣了诫尽，其余都是一樣的，運(yùn)行代碼后我們的隨機(jī)斜率增長(zhǎng)混合模型就出來了

summarytable(gmm1_2, gmm2_2, gmm3_2)

可以看到模型二表現(xiàn)得最好炬守，增長(zhǎng)混合模型驗(yàn)證了我們數(shù)據(jù)分為兩類的假設(shè)牧嫉，我們具體看看模型二的結(jié)果輸出：

輸出的結(jié)果中展示了兩個(gè)類別時(shí)的每一類中時(shí)間的固定效應(yīng)和和第二類隨機(jī)效應(yīng)，同樣的類別1的隨機(jī)效應(yīng)可以通過矩陣和比例系數(shù)相乘得到劳较。

以上就是所有的潛增長(zhǎng)和增長(zhǎng)混合模型的做法介紹驹止。

如何在論文中報(bào)告結(jié)果

看了一些些水水的論文哈，基本都是報(bào)告幾個(gè)擬合優(yōu)度指數(shù)观蜗，還有軌跡的時(shí)間系數(shù)，然后討論一番：

上面的結(jié)果描述來自于一篇中國(guó)青少年抑郁癥狀軌跡的研究衣洁。大家可以參考墓捻。在我們的結(jié)果輸出中都是有相應(yīng)的系數(shù)的。

小結(jié)

今天給大家寫了常見軌跡的做法坊夫，希望對(duì)大家有用砖第，感謝大家耐心看完，自己的文章都寫的很細(xì)环凿，代碼都在原文中梧兼，希望大家都可以自己做一做，請(qǐng)轉(zhuǎn)發(fā)本文到朋友圈后私信回復(fù)“數(shù)據(jù)鏈接”獲取所有數(shù)據(jù)和本人收集的學(xué)習(xí)資料智听。如果對(duì)您有用請(qǐng)先收藏羽杰，再點(diǎn)贊分享。

也歡迎大家的意見和建議到推，大家想了解什么統(tǒng)計(jì)方法都可以在文章下留言考赛，說不定我看見了就會(huì)給你寫教程哦，另歡迎私信莉测。