1、重復測量資料的統(tǒng)計分析方法簡介在醫(yī)學研究中，一些干預研究和縱向研究都需要對研究對象進行隨訪，每次隨訪進行觀測或測量一些效應指標，考察同一研究對象同一指標的變化情況。同一個對象的多次觀察或測量所獲得的資料稱為重復測量的資料。由于同一對象同一指標的相鄰兩個時間點的效應指標觀測值往往是相關的，也就是重復測量的資料存在不獨立的問題，然而大多數(shù)的醫(yī)學統(tǒng)計方法都要求資料是獨立，所以這些資料的統(tǒng)計分析需要用比較特殊的統(tǒng)計方法進行分析。重復測量資料的統(tǒng)計分析方法可以用重復測量的方差分析，也可以用混合回歸模型（Mixed regression Model），由于重復測量的方差分析要求資料滿足球形對稱性（可以理解

2、為相關資料情況下的方差齊性），而Mixed回歸模型并不要求資料滿足球形對稱，并可以借助計算機統(tǒng)計軟件對未知參數(shù)進行限制的最大似然估計，其他統(tǒng)計分析的思想都是類似的。本節(jié)將主要介紹如何借助統(tǒng)計軟件應用Mixed回歸模型對重復測量資料進行統(tǒng)計分析。為了幫助讀者對重復測量資料分析有一個簡單的了解，本節(jié)將舉一個非常簡單的例子初步說明重復測量資料的統(tǒng)計分析概況。例1 為了比較A藥和B藥在療程為6個月中的持續(xù)減肥的療效，現(xiàn)有10個身高為160cm的女性肥胖者志愿參加這項研究。隨機分成2組，每組各5人。分別考察這2組肥胖者在服藥前、服藥3個月和服藥6個月的體重變化。這2組肥胖者在服用該藥前、服藥3個月和服藥

3、6個月的體重測量值(kg)見表1。表 1 2組肥胖者在服用該藥前、服藥3個月和6個月的體重服藥前 3個月6個月組別和肥胖者編號（t1=0,t2=0） （t1=1,t2=0） （t1=0,t2=1）A藥組1號 A藥組2號 A藥組3號 A藥組4號 A藥組5號 B藥組1號 B藥組2號 B藥組3號 B藥組4號 B藥組5號52 51 50 51 49 51 49 50 49 5249 50 49 49 47 54 47 47 48 5042 46 41 44 40 53 46 44 41 48這是兩組觀察對象的多個測量時間點的重復觀察測量資料，同樣對于同一對象的不同觀察時間點的觀察資料是相關的，但由于需

4、要比較兩個藥的減肥療效，所以采用兩因素方差分析，隨機區(qū)組設計的方差分析或Friedman秩檢驗的統(tǒng)計方法都不適用于本例的數(shù)據(jù)統(tǒng)計分析，但可用Mixed模型對本例資料進行統(tǒng)計分析。設A藥組對象在服藥前體重總體均數(shù)為0；服藥3個月后的總體體重改變量為1，故服藥3個月時的體重總體均數(shù)為0+1；服藥6個月時，體重比服藥前的總體改變量為2，即服藥6個月時的體重總體均數(shù)為0+2；同理，設B藥組對象在服藥前體重總體均數(shù)為10；服藥3個月后的總體體重改變量為11，故服藥3個月時的體重總體均數(shù)為10+11；服藥6個月時，體重比服藥前的總體改變量為12，即服藥6個月時的體重總體均數(shù)為10+12；為了便于兩組比較，

5、引入兩組比較的差異參數(shù)如下：記服藥前的兩組差異為010，即：B組服藥前的總體均數(shù)可以表示為10=0+3；在服藥3個月時A藥和B藥的體重總體改變量分別為1和11 ，記服用B藥和A藥3個月時的體重總體改變量的差異為4=11-1，即在服藥3個月時B藥的體重總體改變量可以表示為1+4；在服藥6個月時A藥和B藥的體重總體改變量分別為2和12 ，記服用B藥和A藥6個月時的體重總體改變量的差異為5=12-2，即在服藥6個月時B藥的體重總體改變量可以表示為2+5，把兩組差異的參數(shù)代入上述表達式，得到下列2組肥胖者在服用該藥前、服藥3個月和6個月的體重總體均數(shù)表達式如表2所示：表2 兩組3個時點總體均數(shù)表達式A

6、組總體均數(shù)(g=0) B組總體均數(shù)(g=1)服藥前(t1=0，t2=0)服藥3個月 (t1=1，t2=0)服藥6個月 (t1=0，t2=1)0 0+1 0+20+3 0+1+3+4 0+2+3+5不難驗證表2的總體均數(shù)表達式可以用下列總體回歸方程(式1)表示：y=0+1t1+2t2+3g+4gt1+5gt2 (1)由于不同對象之間存在個體差異i,同一對象不同時點之間也存在隨機差異it,因此第g組第i個對象第t時刻的體重觀察值可以用式（2）表示為ygit=0+1t1+2t2+3g+4gt1+5gt2+gi+git （2）并且假定giN(0,)，gitN(0,)，稱式（2）為混合線性模型（Mixe

7、d Model）。 若4和5不全為0，則稱兩種藥物與服藥時間對療效有交互作用。兩組在3個時間點的總體均數(shù)差異分別為3，3+4和3+5，因此只需檢驗H0：3=0、H0：3+4=0和 H0：3+5=0就可以推斷兩組總體均數(shù)差異。反之若4和5全為0，則稱兩種藥物與服藥時間對療效無交互作用，并且兩組各個時間點的總體均數(shù)差異均為3，因此只需檢驗H0：30就可以推斷兩組的總體均數(shù)差異。我們同樣借助Stata軟件對上述資料用混合模型進行統(tǒng)計分析，相應的Stata軟件的數(shù)據(jù)格式如下Stata操作命令如下：gen gt1=g*t1 產(chǎn)生交互作用項變量gt1gen gt2=g*t2 產(chǎn)生交互作用項變量gt2由此得

8、到下列統(tǒng)計結論：在服藥前，A組的體重總體均數(shù)0的估計值為50.6(kg)服藥3個月時，A組的體重改變量的總體均數(shù)為1的估計值為-1.8(kg)，P值0.088>0.05，因此沒有足夠證據(jù)推斷在服用A藥3個月時，服用A藥的人群平均體重發(fā)生改變。服A藥6個月時，A組的體重改變量的總體均數(shù)為2的估計值為-8(kg)，P值<0.001，因此可以認為在服用A藥6個月時，服用A藥人群的平均體重已有下降。服藥前兩組平均體重的差異3的估計值為0.4（kg）,相應的P值為=0.804>0.05，差別無統(tǒng)計學意義。服藥3個月時，服藥A藥和服擁B藥的兩個人群平均體重下降的差異4的估計值為0.8，P

9、值=0.591，差異無統(tǒng)計學意義。服藥6個月時，服藥A藥和服擁B藥的兩個人群平均體重下降的差異5的估計值為4.2，P值=0.005<0.05，差異有統(tǒng)計學意義。服用A藥和B藥在各個時間的總體均數(shù)估計見表3。表3 在各個時間點的兩組總體均數(shù)及其估計值個月時A組的總體均數(shù)估計值50.6-1.848.8。服藥后的2個時間點的兩組平均體重比較的Stata統(tǒng)計檢驗命令和輸出結果如下 設=0.05服藥3個月時的兩組平均體重比較的Stata命令和輸出結果如下：test ggt10 （H0：3+40 即服藥3個月時的兩組體重總體均數(shù)相等） 相應的P值0.8041>，差異無統(tǒng)計學意義，故無證據(jù)顯示兩

10、組總體均數(shù)不等。 服藥6個月時的兩組平均體重比較的Stata命令和輸出結果如下： test ggt20 （H0：3+50）相應的P值0.0184<，差異有統(tǒng)計學意義，故可以認為服A藥6個月時的人群平均體重低于服B藥6個月的人群平均體重。（四）隨機效應的Logistic模型在第十章中介紹了常見的兩分類logistic模型為Pr(Y=1|x1,x2, ,xp)=或者可以寫為exp(0+1x1+ +mxm)（3）1+exp(0+1x1+ +mxm)Plogit(P)=log =0+1x1+ +mxm （4） 1-p并由此可知：P越大，logit(P)也越大；P越小，logit(P)也越小，并且

11、兩者一一對應。logistic模型的主要適應范疇為：應變量為兩分類變量，每個觀察記錄是獨立的。為了使讀者能較連貫地理解隨機效應的logistic模型，故通過下面一個簡單例子復習相關內(nèi)容。例2 收集100名HP陽性患者，隨機分成兩組，分別用兩種治療方案進行治療（A方案g=1，B方案g=0），2個月后考察HP治療的療效（陰性Y=1,仍為陽性y=0），比較兩組療效的差別。表4 兩組100名HP陽性患者治療結果A組(g=1)B組(g=0) 陰性(Y=1) 19(38%) 13(26%) 陽性(Y=0) 31(62%) 37(74%) 合計 50 50由于本例的應變量為二分類變量，所以可以用Logist

12、ic模型對本例資料進行統(tǒng)計分析。本例的Logistic模型為Logit(P)=0+1g （5）對于A組，g=1，logit(PA)01；對于B組，g=0，logit（PB）0，因此若1>0，說明logit(PA)> logit（PB），相應的A組轉陰率PA高于B組轉陰率PB；若1<0，說明logit(PA)< logit（PB），相應的A組轉陰率PA低于B組轉陰率PB；若10，說明兩組轉陰率相同。我們?nèi)杂肧tata軟件進行參數(shù)估計和統(tǒng)計檢驗。Stata數(shù)據(jù)格式為（數(shù)據(jù)文件為ex12-5.dta）Stata命令： logit y g freq=w得到0的估計值為-1.04

13、5969，1的估計值為0.5564203，1的P值=0.200，因此對于檢驗水準0.05而言，兩組療效的差異無統(tǒng)計學意義。故沒有充分證據(jù)說明兩組療效有差異。由于HP經(jīng)常呈現(xiàn)反復，這次檢查可能是陰性，但下次檢查可能又呈陽性了，所以往往要考察幾個療程才能評價兩種治療方案的療效。例3 收集100名HP陽性患者，隨機分成兩組，分別用兩種治療方案進行治療（A方案g=1，B方案g=0），2個月為一個療程并檢查一次HP，連續(xù)治療3個療程，考察兩種治療HP方案的療效（陰性y=1,仍為陽性y=0），比較兩組療效的差別。表5 兩種治療方案3個療程的HP轉陰情況一個療程t1=0，t2=0 二個療程 t1=1，t2=

14、0 三個療程 t1=0，t2=1 A組（g=0） 人數(shù) B組(g=1) 人數(shù)陽性(y=0) 陽性(y=0) 陽性(y=0) 陽性(y=0) 陰性(y=1) 陰性(y=1) 陰性(y=1) 陰性(y=1) 陽性(y=0) 陽性(y=0) 陰性(y=1) 陰性(y=1) 陽性(y=0) 陽性(y=0) 陰性(y=1) 陰性(y=1) 陽性(y=0) 陰性(y=1) 陽性(y=0) 陰性(y=1) 陽性(y=0) 陰性(y=1) 陽性(y=0) 陰性(y=1)1 3 7 20 1 5 5 8 1 4 1 15 1 3 4 21由于患者轉陰后還有一定的概率轉陽，因此不能用生存分析的方法進行比較，但本例

15、資料可以認為是重復測量的二分類資料，可以用下列隨機效應的Logistic模型：.Pr(Y=1|x1,x2, ,xm)=exp(0+1x1+ +mxm+u)1+exp(0+1x1+ +mxm+u)（6）其中uN(0,2)該模型的特點是常數(shù)項含有隨機誤差項以及隨機效應與協(xié)變量無關以及與回歸系數(shù)無關。也可以寫成log(Odds)=Logit(P|x1,x2, ,xm)=0+1x1+ +mxm+u（7）0，1，m為參數(shù)，一般采用限制的最大似然估計，由于這種估計方法非常復雜，本書不作介紹，我們僅介紹借助Stata軟件對參數(shù)進行估計和檢驗。本例模型可以表示為（8）式，其中變量定義見表12.9。Logit(

16、P)=0+1t1+2t2+3g+4gt1+5gt2+u（8）其中u服從N(0,2)，0+1t1+2t2+3g+4gt1+5gt2稱為Logit(P)的確定性部分，我們主要關心的是Logit(P)的確定性部分，其與各個變量之間的對應列表如下。表6 兩種治療方案與多個療程的Logit(P)的確定性部分關系方案A（g=1） 方案B（g=0）兩種方案Log(Odds)的差異Odds Ratio一個療程（t1=0,t2=0） 二個療程（t1=1,t2=0） 三個療程（t1=0,t2=1）030 3 exp(3)013401 34 exp(34)023502 35 exp(35)因此兩種治療方案在一個療程

17、治療的轉陰率比較就是檢驗H0：30；在二個療程治療的轉陰率比較就是檢驗H0：340；在三個療程治療的轉陰率比較就是檢驗H0：350。如果兩種治療方案與隨訪時間無交互作用，即40且50，則上述模型可以簡化為 （9）式，相應的兩種治療方案與多個療程的Logit(P)的確定性部分關系見表7。Logit(P)=0+1t1+2t2+3g+u方案A（g=0） 方案B（g=1）一個療程（t1=0,t2=0）0 03二個療程（t1=1,t2=0）01 013（9）表7 兩種治療方案與多個療程的Logit(P)的確定性部分關系三個療程（t1=0,t2=1）02 023兩種方案的差異3 3 3第二療程與第一療程比

18、較：B方案：（013）（03）1 A方案：010=1 第三療程與第一療程比較：B方案：（023）（03）2 A方案：020=2第三療程與第二療程比較：B方案：（023）（013）=2-1 A：方案：02(0+1)=2-1因此各個療程的兩種治療方案轉陰率比較就是檢驗H0：30。本例數(shù)據(jù)和Stata數(shù)據(jù)結構如下本例的Stata操作命令gen gt1=g*t1 產(chǎn)生g×t1項 交互項4的估計值為-0.6822685，相應的P值為0.297，5的估計值為0.0585067，相應的P值為0.929，因此我們需要檢驗H0：450 vs H1：4和5不全為0，=0.05,相應的Stata操作命令如下交互作用logistic模型，相應的Stata命令為公式（12.9）的參數(shù)3的估計值為0.6409473，相應的OR=e0.641=1.898，P值為0.018，對于統(tǒng)計檢驗水準為0.05而言，兩組療效的差別有統(tǒng)計學意義，可以認為B方案的療效優(yōu)于A方案；參數(shù)1的估計值為1.563917，相應的OR=e1.564=4.778，P值<0.001，因此可以認為二個療程的HP轉陰率高于一個療程，并且差別有統(tǒng)