4.11 重復測量設計一元定量資料方差分析
4.11.1 問題與數據
【例4-12】某課題組研究益髓生血顆粒治療_SEA/ααCS基因型患者血紅蛋白H（Hemoglobin H，HbH）病的臨床療效，選取13名_SEA/ααCS基因型HbH病患者，給以益髓生血顆粒治療，療程為3個月。分別記錄患者治療前、服藥1個月、服藥2個月及服藥3個月時血紅蛋白的含量，結果如表4-14所示。請進行合適的統計分析。
表4-14 益髓生血顆粒治療前后_SEA/ααCS基因型HbH患者血紅蛋白含量
【例4-13】某研究者欲觀察隔日禁食結合耐饑食餌對小鼠生理指標的影響，將20只成年ICR小鼠隨機地等分為對照組（自由攝食標準飼料）、食餌組（自由攝食耐饑食餌）、隔日禁食組（隔日禁食標準飼料）、辟谷食餌組（隔日禁食耐饑食餌），實驗8周，自由飲水，觀察4組小鼠體重的變化，如表4-15所示。請進行適當的統計分析。
表4-15 4組小鼠體重的變化
【例4-14】某研究者用賁門癌患者的標本制成液體，在3種不同處理條件下觀察雞胚背根神經節與雞胚交感神經節中長出突起的神經節比例。現有賁門癌患者10例，將每人的標本均分成3份，分別給予3種不同的處理（因素A），即A1（加入100ng/ml神經生長因子）、A2（加入200ng/ml神經生長因子）和A3（單用賁門癌培養液）；并對每種處理后的標本中的2種類型的神經節（因素B），即B1（背根神經節）與B2（交感神經節），觀測長出突起的神經節的比例（Y）。實驗結果如表4-16所示。請進行適當的統計分析。
表4-16 賁門癌患者的標本經3種處理后2種神經節中長出突起的神經節的比例
4.11.2 對數據結構的分析
在例4-12中，對每一個患者在4個時間點上分別測量其血紅蛋白含量，說明“時間”因素是一個重復測量因素，并且所有患者均接受同一種治療方法—服用益髓生血顆粒，因而這是具有一個重復測量的單因素設計定量資料。
在例4-13中，對每一只小鼠在4個時間點上分別測量其體重，說明“時間”因素是一個重復測量的因素。此外，還有2個實驗因素“攝食方式”（自由攝食或隔日禁食）和“飼料種類”（標準飼料或耐饑食餌），因而這是具有一個重復測量的三因素設計定量資料。
在例4-14中，涉及2個實驗因素，即“處理”和“神經節類型”，前者有3個水平，后者有2個水平。由于是對同一個患者的標本采用不同的處理，同時在不同類型神經節中重復獲得指標“長出突起的神經節的比例”的6個觀測值，因而2個實驗因素都是重復測量因素。該資料所對應的實驗設計類型為具有2個重復測量的兩因素設計。
4.11.3 分析目的與統計分析方法的選擇
分析目的都是比較不同實驗條件下定量觀測指標的平均值之間的差別是否具有統計學意義。
例4-12，應選用具有一個重復測量的單因素設計一元定量資料方差分析；例4-13，應選用具有一個重復測量的三因素設計一元定量資料方差分析；例4-14，應選用具有2個重復測量的兩因素設計一元定量資料方差分析。
4.11.4 SAS程序
對例4-12進行具有一個重復測量的單因素設計一元定量資料方差分析，SAS程序名為SASTJFX4_12.SAS。
SAS程序中第1步為建立數據集，person代表“患者編號”，month代表“觀測時間”，y代表觀測指標“血紅蛋白含量”。第2、3、4、5、6步分別調用MIXED過程，采用VC、CS、UN、AR(1)、SP(POW)5種協方差結構模型對資料進行方差分析。其中，SP(POW)模型在使用時要求將“repeated/type=sp(pow) (c-list)”語句中的“c-list”替換成重復測量因素（可以是多個）的名稱，但這個或這些重復測量因素應為定量變量且其水平數不能太少。在本例中，重復測量因素為觀測時間，其包含4個水平（0、1、2、3），可近似視為定量變量，故可選用SP(POW)模型。若某重復測量因素為定性變量（如部位等）或雖可視為定量變量但水平數較少（如小于等于3）時，是不適宜選用此模型的。此外，定量的重復測量因素在賦值時需保證賦值數值的間隔與實際觀測時間間隔成比例，較好且最簡單的方法是以其真實水平代入。如本例，可以0、1、2、3代表治療前、服藥1個月、服藥2個月、服藥3個月，而不宜隨意賦值，如賦為1、3、5、9則不妥。
對例4-13進行具有一個重復測量的三因素設計一元定量資料方差分析，SAS程序名為SASTJFX4_13.SAS。
SAS程序第1步為建立數據集，style代表“攝食方式”，forage代表“飼料種類”，time代表“時間”，subject代表“小鼠個體”，y代表觀測指標“體重”。第2、3、4、5、6步分別調用MIXED過程，采用VC、CS、UN、AR(1)、SP(POW)5種協方差結構模型對資料進行方差分析，其中model語句后的“style|forage|time”相當于“style forage time style*forage style*time forage*time style*forage*time”，即這三個因素的主效應及各級交互作用的效應。第7步為建立宏shuju，以實現對數據集中已有變量value的更名。第8、9步均用來實現對不同數據集的橫向合并。第10、11步均用來將數據集中的內容輸出到output窗口中。第2、3、4、5、6步分別調用MIXED過程，幾個過程所用語句基本相同，僅在“type=”后的選項不同，5個過程分別指定了5種協方差結構模型。MIXED過程中的repeated語句用來規定個體的重復測量的協方差結構，“/”后的sub（也可寫為subject）用來指定數據集中的個體，若不含有分組因素，直接在“sub=”后面給出受試對象個體變量名稱即可，如程序SASTJFX4_12；若含有分組因素，則在“sub=”后面給出受試對象個體變量名稱的同時，還需在后面加注“（）”，在括號內填入分組變量名稱，如程序SASTJFX4_13。在調用MIXED過程進行方差分析時，使用了兩個ods語句，分別用來將模型擬合的有關信息（fitstatistics）和模型維度有關參數（dimensions）輸出。
對例4-14進行具有兩個重復測量的兩因素設計一元定量資料方差分析，SAS程序名為SASTJFX4_14.SAS。
SAS程序第1步為建立數據集，patient代表“病例號”，treat代表“處理”，type代表“神經節類型”，y代表觀測指標“長出凸起的神經節的比例”。第2、3、4、5步分別調用MIXED過程，采用VC、CS、UN、AR(1)4種協方差結構模型對資料進行方差分析。第6、7步將不同數據集的內容進行合并查詢，這與程序SASTJFX4_13中的宏shuju功能相同。
4.11.5 主要分析結果及解釋
下面給出例4-12的分析結果，即程序SASTJFX4_12.SAS的部分輸出結果。
以下是采用VC（方差分量型）協方差結構模型進行混合效應模型分析的輸出結果。
以下是采用CS（復合對稱型）協方差結構模型進行混合效應模型分析的輸出結果。
以下是采用UN（無結構型）協方差結構模型進行混合效應模型分析的輸出結果。
以下是采用AR(1)（一階自回歸型）協方差結構模型進行混合效應模型分析的輸出結果。
這是采用SP（POW）（空間冪相關型）協方差結構模型進行混合效應模型分析的輸出結果。
以上是采用MIXED過程進行混合效應模型分析的結果。這里給出了模型的擬合信息和固定效應假設檢驗的結果。通常關注的是固定效應假設檢驗的結果。但從上面的結果中可以看出，采用不同的協方差結構模型得出的固定效應假設檢驗的結果是不完全相同的，其中GLM過程計算出來的結果與CS協方差結構模型（采用REML即約束最大似然估計法）得出的結果是相同的，它是最簡單的模型。在采用這5種協方差結構模型計算所得的結果中，應以哪個結果為準呢？這就是模型選擇的問題了。
通常情況下，可以Akaike的信息準則AIC值或Schwarz的信息準則BIC（或稱SBC）值來選擇協方差結構模型。AIC值和BIC值越小，協方差結構模型擬合越好；若兩個協方差結構模型擬合的AIC值和BIC值接近，還可參考?2 LogL（?2 Res Log Likelihood）的數值，小者為優。若兩個協方差結構模型分別包含q+v和q個參數時，可用這兩個模型的?2倍的對數似然函數值構造出似然比統計量，采用χ2檢驗進行推斷，見式（4-1）。若似然比統計量對應的P值大于設定的臨界值，則兩個協方差結構模型對資料的擬合效果之間的差異無統計學意義，此時可選擇參數個數較少的那個協方差結構模型。若似然比統計量對應的P值小于設定的臨界值，則兩個協方差結構模型對資料的擬合效果之間的差異有統計學意義，此時可選擇參數個數較多的那個協方差結構模型。式中，χ2v服從自由度為v的χ2分布，?2 logLq和?2 logLq+v分別為含q和q+v個協方差參數模型的?2倍的對數似然函數值。
現將本資料輸出結果的模型擬合信息部分的主要內容匯總在表4-17中，以利于比較。
表4-17 用5種類型的協方差結構模型擬合本資料的擬合效果比較
注：各模型中待估計的協方差結構中參數的個數依次為1、2、10、2、2。
由上面的結果可知：采用UN（無結構型）協方差結構模型擬合本資料，其迭代過程并沒有收斂，擬合效果不好。對比其他4種協方差結構模型的擬合情況，可以看出采用CS協方差結構模型擬合本資料效果最好，因為在評價擬合效果的三個準則中，CS模型對應的值均為最小且其?2 Log L值較參數個數比其少的模型VC的值差異也比較大，所以此時應以CS（復合對稱型）協方差結構模型進行混合模型分析的輸出結果為準，其中F=1267.75，P<0.0001。說明不同時間點上測得的HbH患者血紅蛋白含量均值之間的差異存在統計學意義。欲得出在各時間點上患者血紅蛋白含量均值之間兩兩比較的結果，可將原程序中所有過程步刪除，在data步后換上以下過程步：
ods html;
proc mixed;
class person month;
model y=month/ddfm=sat;
repeated/type=CS sub=person;
lsmeans month/pdiff;
run;
ods html close;
其中，“ddfm=sat”說明采用Satterthwaite近似的方法來計算分母的自由度，從而得到一個更加準確的近似程度更高的F值。但是，此法計算的結果有時會與GLM過程計算的結果有所不同。
程序運行后，得結果如下：
這是各時間點上患者血紅蛋白含量均值之間兩兩比較的結果。首先，給出了各時間點上患者血紅蛋白含量均值與0比較的檢驗結果，沒有太大實際意義。然后，給出了兩兩比較的結果，month和_month列中的0、1、2、3分別代表時間的4個水平，即治療前、服藥1個月、服藥2個月、服藥3個月。讀者可查看最后兩列給出的兩兩比較的結果，可發現這4個時間點上患者血紅蛋白含量均值之間的差異均有統計學意義，說明4個時間點上患者血紅蛋白含量各不相同，根據上面給出的均值大小可判斷，各時間點患者血紅蛋白含量的大小關系為：服藥1個月>服藥3個月>服藥2個月>治療前。
以下是對例4-13的分析結果，即程序SASTJFX4_13.SAS的輸出結果。
這是上述程序中ODS（Output Delivery System）輸出的結果。首先，給出了5種協方差結構模型擬合本資料的有關情況，然后給出了協方差結構的有關信息（Covariance Parameters表示模型中待估計的協方差結構中的參數個數）。比較5種模型擬合資料情況的AIC、BIC數值，可發現UN與AR(1)兩種協方差結構模型擬合資料較好。現比較這兩種協方差結構模型擬合資料的效果之間的差異是否有統計學意義。
χv2=?2 logLq?(?2 logLq+v)=235.5?219.1=16.4
由ODS輸出結果的第2部分可知，q=2，q+v=10，所以v=8。因χ0.05(8)2=15.51<16.4，故P<0.05。因為可認為不適合用AR(1)模型取代UN模型，所以最后的結論應按UN協方差結構模型計算出來的結果來下。其假設檢驗結果為：
由上述結果可知：攝食方式（style）、觀測時間（time）、攝食方式與觀測時間的交互作用（style*time）、飼料種類與觀測時間的交互作用（forage*time）均有統計學意義。
欲得出各均數之間兩兩比較的結果，可將原程序中所有過程步刪除，換上以下過程步：
ods html;
proc mixed data=sastjfx4_13;
class style forage subject time;
model y=style|forage|time/ddfm=sat;
repeated/type=UN sub=subject(style forage);
lsmeans style*forage*time/pdiff;
run;
ods html close;
此步運行結果大部分與上述TYPE=UN計算結果相同，僅多了多個均數兩兩比較的結果，因輸出結果所占篇幅較多，故此處從略。
其實，對于含分組因素的重復測量設計定量資料，如果觀測了各受試對象某指標的“基礎”狀態（即實驗前狀態），對資料進行分析時，最合適的做法是以“基礎值”為定量的影響因素，即協變量，采用相應設計類型定量資料的協方差分析。例如本例，若第一個時間點觀測的數據不是實驗后第1周小鼠的體重，而是實驗前小鼠的體重，則本資料最合適的分析方法是以實驗前小鼠的體重為基礎值（協變量），采用具有一個重復測量的三因素設計定量資料一元協方差分析。有關做法可參閱4.12節。
以下是對例4-14的分析結果，即程序SASTJFX4_14.SAS的輸出結果。
這是上述程序中ODS輸出的結果。首先給出了4種協方差結構模型擬合本資料的有關情況，然后給出了協方差結構的有關信息。比較4種模型擬合資料情況的AIC、BIC數值（越小越好），可發現VC協方差結構模型擬合資料較好且協方差結構中參數的個數最少。因此，最后的結論應按VC協方差結構模型計算出來的結果來下。其假設檢驗結果為：
由上述結果可知：處理（treat）和神經節類型（type）均有統計學意義。
欲得出各均數之間兩兩比較的結果，可將原程序中所有過程步刪除，換上以下過程步：
ods html;
proc mixed data=sastjfx4_16;
class patient treat type;
model y=treat|type/ddfm=sat;
repeated/type=VC sub=patient;
lsmeans treat*type/pdiff;
run;
ods html close;
此步運行結果大部分與上述type=VC計算結果基本相同，僅多了多個均數兩兩比較的結果。因輸出結果所占篇幅較多，故此處從略。