如何用 SAS 軟件正確分析生物醫(yī)學(xué)科研資料XIV.用 SAS 軟件實現(xiàn)具有一個重復(fù)測量的單因素設(shè)計定量資料的統(tǒng)計分析

王琪，胡良平

1 重復(fù)測量設(shè)計的概念及應(yīng)用場合

重復(fù)測量設(shè)計是在不同條件下，從同一受試對象上重復(fù)獲得某指標(biāo)觀測值的一種試驗設(shè)計類型。這里所說的“在不同條件下”，通常是指“時間因素”取不同水平，有時指受試者身上的幾個“對稱部位”或具有可比性的不同部位，有時也指“時間因素”和“對稱部位”的各種水平組合。如果“不同條件”僅與一個因素有關(guān)，就叫做“具有一個重復(fù)測量的M因素設(shè)計”；如果“不同條件”與兩個因素的水平組合有關(guān)，就叫做“具有兩個重復(fù)測量的M因素設(shè)計”。這里的M因素指試驗中涉及到的全部試驗因素的個數(shù)，包括與重復(fù)測量有關(guān)的試驗因素。與重復(fù)測量無關(guān)的因素被稱為試驗分組因素，即受試對象被完全隨機(jī)地分配到這些因素的水平組合所形成的各小組中去。

在某項研究中，若專業(yè)上需要了解隨時間推移或部位改變時，定量觀測指標(biāo)的動態(tài)變化情況時，就需要運(yùn)用此設(shè)計。

2 具有一個重復(fù)測量的單因素設(shè)計定量資料的SAS 實現(xiàn)

【例 1】某課題組研究益髓生血顆粒治療_SEA/ααCS基因型患者血紅蛋白 H（hemoglobin H，HbH）病的臨床療效，選取 13 名_SEA/ααCS基因型 HbH 患者，給以益髓生血顆粒治療，療程為3 個月。分別記錄患者治療前、服藥 1 個月、服藥 2 個月及服藥 3 個月時血紅蛋白的含量，結(jié)果見表1。請對資料進(jìn)行合適的統(tǒng)計分析。

表1 13 例 HbH 患者服藥前后血紅蛋白含量變化

對數(shù)據(jù)結(jié)構(gòu)的分析與 SAS 實現(xiàn)：對每一位患者來說，在 4 個時間點上分別被測量血紅蛋白含量，說明“觀測時間”是一個重復(fù)測量的因素，且所有患者均接受同一種治療方法——服用益髓生血顆粒，因而這是具有一個重復(fù)測量的單因素設(shè)計定量資料。直接采用隨機(jī)區(qū)組設(shè)計定量資料方差分析處理此數(shù)據(jù)是不妥的，因為用隨機(jī)區(qū)組設(shè)計定量資料方差分析處理具有一個重復(fù)測量的單因素設(shè)計定量資料有特殊的前提條件——球?qū)ΨQ條件。當(dāng)球?qū)ΨQ條件滿足時，采用隨機(jī)區(qū)組設(shè)計定量資料方差分析處理與采用具有一個重復(fù)測量的單因素設(shè)計定量資料方差分析處理的結(jié)果一致；當(dāng)球?qū)ΨQ條件不滿足時，若采用隨機(jī)區(qū)組設(shè)計定量資料方差分析來處理，與采用具有一個重復(fù)測量的單因素設(shè)計定量資料方差分析時未校正的結(jié)果一致，當(dāng)然會增大犯 I 類錯誤的概率。應(yīng)先檢查資料是否滿足球?qū)ΨQ條件，若滿足，可采用隨機(jī)區(qū)組設(shè)計定量資料方差分析；若不滿足，可采用與此設(shè)計對應(yīng)的混合效應(yīng)模型處理。

94.98 102.33 98.86 102.31 91.41 98.79 95.46 98.66 78.75 86.22 83.38 85.70 83.13 90.58 87.56 90.19 71.20 78.73 76.18 77.98 87.28 94.69 91.51 94.43 65.01 72.58 70.27 71.64 66.93 74.49 72.11 73.61 103.31 110.61 106.81 110.84 80.87 88.32 85.40 87.87 71.97 79.49 76.91 78.77 65.20 72.77 70.46 71.84 79.16 86.63 83.77 86.12

程序說明：SAS 程序中第1 步為建立數(shù)據(jù)集，patient代表“患者編號”，m0、m1、m2、m3 分別代表治療前、服藥 1 個月、服藥 2 個月、服藥 3 個月4 個時間點獲得的觀測指標(biāo)“血紅蛋白含量”的值。第2 步是調(diào)用 GLM 過程進(jìn)行具有一個重復(fù)測量的單因素設(shè)計定量資料一元方差分析和多元方差分析，其中 model 語句“/”后的 nouni 用來指明不要將 m0、m1、m2、m3 看成 4 個相互獨立的結(jié)果變量進(jìn)行方差分析，repeated 語句用來指定與重復(fù)測量有關(guān)的因素及其水平數(shù)，并對資料進(jìn)行一元方差分析和多元方差分析。在進(jìn)行一元方差分析時，必須對協(xié)方差陣進(jìn)行球?qū)ΨQ性檢驗，這只需在 repeated 語句“/”后加上 printe 即可，選項 summary 可輸出方差分析表。第3 步為調(diào)整數(shù)據(jù)集，將原本單獨成列的 m0、m1、m2、m3 轉(zhuǎn)換成變量 month 的4 個水平，并將其取值全部賦給變量 y，后面幾步均針對此新產(chǎn)生的數(shù)據(jù)集 a 進(jìn)行分析。第4 步對新數(shù)據(jù)集 a 進(jìn)行隨機(jī)區(qū)組設(shè)計定量資料的方差分析。第5、6、7、8、9 步分別調(diào)用 MIXED 過程，采用 VC、CS、UN、AR(1)、SP（POW）五種協(xié)方差結(jié)構(gòu)模型對資料進(jìn)行方差分析。其中，SP（POW）模型在使用時要求將“repeated/type = sp(pow)(c-list)”語句中“c-list”替換成重復(fù)測量因素的名稱，但此重復(fù)測量因素應(yīng)為數(shù)值型變量。本例中，重復(fù)測量因素為觀測時間，其包含 4 個水平（0、1、2、3），屬數(shù)值型變量，故可選用 SP（POW）模型。若某重復(fù)測量因素為定性變量（如部位等）或雖可視作定量變量但水平數(shù)較少（如 ≤ 3）時，是不適宜選用此模型的。此外，在使用 SP（POW）模型時，定量的重復(fù)測量因素在賦值時應(yīng)以其真實水平代入。

SAS 輸出結(jié)果與結(jié)果解釋：

Sphericity tests

The GLM procedure repeated measures analysis of variance univariate tests of hypotheses for within subject effects

Greenhouse-Geisser epsilon 0.3334 Huynh-Feldt epsilon 0.3334

這是輸出結(jié)果的第一部分，是第一個 GLM 過程產(chǎn)生的具有一個重復(fù)測量的單因素設(shè)計定量資料一元方差分析的結(jié)果。首先給出的是對協(xié)方差陣進(jìn)行球?qū)ΨQ性檢驗的結(jié)果，由P< 0.0001 可知，此資料不滿足球?qū)ΨQ性條件，因而在查看方差分析結(jié)果時，應(yīng)查看校正后的P值。本資料采用 Greenhouse-Geisser 法和 Huynh-Feldt 法校正后的P值均小于 0.0001，說明各時間點上測得的 HbH 患者血紅蛋白含量均值之間的差異存在統(tǒng)計學(xué)意義。當(dāng)協(xié)方差矩陣滿足球?qū)ΨQ性要求時，（Epsilon）（對P值進(jìn)行校正的系數(shù)）等于 1，越大，說明協(xié)方差矩陣越接近球?qū)ΨQ性。Huynh-Feldt 法校正的ε? 取值可能大于 1，當(dāng)> 1 時，取= 1。

MANOVA test criteria and exact F statistics for the hypothesis of no month effect

這是輸出結(jié)果的第二部分，是第一個 GLM 過程產(chǎn)生的具有一個重復(fù)測量的單因素設(shè)計定量資料多元方差分析的結(jié)果。分別給出了 Wilks' Lambda、Pillai's Trace、Hotelling-Lawley Trace、Roy's Greatest Root 四個檢驗統(tǒng)計量，它們都是從離均差平方和矩陣的特征根推導(dǎo)出來的，都是特征根的函數(shù)，彼此密切相關(guān)。當(dāng)結(jié)果變量的第一個即最大特征根完全解釋了所有變異時，四個統(tǒng)計量對應(yīng)的F值相等，并精確服從F分布。否則，這四個統(tǒng)計量所對應(yīng)的F值不相等且只是近似服從F分布。在多元方差分析中，通常用 Wilks' Lambda 進(jìn)行統(tǒng)計學(xué)推斷。由上述結(jié)果可知：P< 0.0001，所以各時間點上測得的 HbH 患者血紅蛋白含量均值之間的差異具有統(tǒng)計學(xué)意義。

The GLM procedure repeated measures analysis of variance analysis of variance of contrast variables month_N represents the contrast between the nth level of month and the last

這是輸出結(jié)果的第三部分，是第一個 GLM 過程產(chǎn)生的各時間點與最后一個時間點患者血紅蛋白含量兩兩比較的結(jié)果?？梢?，各時間點與第4 個時間點上患者血紅蛋白平均含量之間的差異均有統(tǒng)計學(xué)意義。當(dāng)然，repeated 語句默認(rèn)的是各時間點與最后一個時間點進(jìn)行比較。其實，在此語句中，可設(shè)置以某個組作為對照組，其他各組與之進(jìn)行比較。如本資料，欲比較各時間點與第1 個時間點上血紅蛋白平均含量之間的差異有無統(tǒng)計學(xué)意義，將“repeated month 4（0 1 2 3）/summary printe;”改為“repeated month 4(0 1 2 3)contrast(1)/summary printe;”即可，即添加 contrast 選項，并在其后的括號內(nèi)規(guī)定何為對照組。

?

R-Square Coeff Var Root MSE yMean 0.999406 0.408029 0.345574 84.69365

Source DF Type III SS Mean square F value PR> F Patient 12 6777.885931 564.823828 4729.66 < .0001 Month 3 454.188898 151.396299 1267.75 < .0001

這是輸出結(jié)果的第四部分，是第二個 GLM 過程產(chǎn)生的隨機(jī)區(qū)組設(shè)計定量資料一元方差分析的結(jié)果。由方差分析結(jié)果可以看出：F= 1267.75，P< 0.0001，說明不同時間點上測得的 HbH 患者血紅蛋白平均含量之間的差異具有統(tǒng)計學(xué)意義。這個方差分析的結(jié)果與具有一個重復(fù)測量的單因素設(shè)計定量資料一元方差分析結(jié)果一致，這是因為資料雖然不滿足球?qū)ΨQ性條件但F值 >>F臨界值，即采用具有一個重復(fù)測量的單因素設(shè)計定量資料一元方差分析和隨機(jī)區(qū)組設(shè)計定量資料一元方差分析處理單因素重復(fù)測量設(shè)計定量資料結(jié)果一致的兩種情況之一。

Fit statistics

Type 3 tests of fixed effects

這是采用 VC（方差分量型）協(xié)方差結(jié)構(gòu)模型進(jìn)行混合效應(yīng)模型分析的輸出結(jié)果。

Fit statistics

Type 3 tests of fixed effects

這是采用 CS（復(fù)合對稱型）協(xié)方差結(jié)構(gòu)模型進(jìn)行混合效應(yīng)模型分析的輸出結(jié)果。

Iteration history

Iteration Evaluations –2 Res Log like Criterion 6 0 –210.88040371 0.00000214 7 0 –210.88040371 0.00000214 8 0 –210.88040371 0.00000214 9 0 –210.88040371 0.00000214 10 0 –210.88040371 0.00000214 11 0 –210.88040371 0.00000214

WARNING: Did not converge.

這是采用 UN（無結(jié)構(gòu)型模型）協(xié)方差結(jié)構(gòu)模型進(jìn)行混合效應(yīng)模型分析的輸出結(jié)果。

Fit statistics

Type 3 tests of fixed effects

這是采用 AR(1)（一階自回歸型）協(xié)方差結(jié)構(gòu)模型進(jìn)行混合效應(yīng)模型分析的輸出結(jié)果。

Fit statistics

Type 3 tests of fixed effects

這是采用 SP（POW）（空間冪相關(guān)型）協(xié)方差結(jié)構(gòu)模型進(jìn)行混合效應(yīng)模型分析的輸出結(jié)果。

以上是采用 MIXED 過程進(jìn)行混合效應(yīng)模型分析的結(jié)果。這里給出了模型的擬合信息和固定效應(yīng)假設(shè)檢驗的結(jié)果。通常關(guān)注的是固定效應(yīng)假設(shè)檢驗的結(jié)果，但從上面的結(jié)果中可以看出，采用不同的協(xié)方差結(jié)構(gòu)模型得出的固定效應(yīng)假設(shè)檢驗的結(jié)果是不完全相同的，其中 GLM 過程計算出來的結(jié)果與 CS 協(xié)方差結(jié)構(gòu)模型（采用 REML 即約束最大似然估計法）得出的結(jié)果是相同的，它是次簡單的模型。那么，采用這 5 種協(xié)方差結(jié)構(gòu)模型計算所得的結(jié)果中應(yīng)以哪個結(jié)果為準(zhǔn)呢？這就是模型選擇的問題了。

通常情況下，可以 Akaike 的信息準(zhǔn)則 AIC 值或Schwarz 的信息準(zhǔn)則 BIC（或稱 SBC）值來選擇協(xié)方差結(jié)構(gòu)模型。AIC 值和 BIC 值越小，協(xié)方差結(jié)構(gòu)模型擬合給定資料越好；若兩個協(xié)方差結(jié)構(gòu)模型擬合的 AIC 值和 BIC值接近，還可參考 –2LogL（–2 Res Log Likelihood）的數(shù)值，小者為優(yōu)。若兩個協(xié)方差結(jié)構(gòu)模型分別包含q+v和q個參數(shù)時，可用這兩個 –2 倍的對數(shù)似然函數(shù)值構(gòu)造出似然比統(tǒng)計量，采用χ2檢驗進(jìn)行推斷，見下面的公式。若似然比統(tǒng)計量對應(yīng)的P值大于設(shè)定的顯著性水平，則兩個協(xié)方差結(jié)構(gòu)模型對資料的擬合效果之間的差異無統(tǒng)計學(xué)意義，此時可選擇參數(shù)個數(shù)較少的那個協(xié)方差結(jié)構(gòu)模型。若似然比統(tǒng)計量對應(yīng)的P值小于設(shè)定的顯著性水平，則兩個協(xié)方差結(jié)構(gòu)模型對資料的擬合效果之間的差異有統(tǒng)計學(xué)意義，此時應(yīng)選擇參數(shù)較多的那個協(xié)方差結(jié)構(gòu)模型。

現(xiàn)將本資料輸出結(jié)果的模型擬合信息部分的主要內(nèi)容匯總在表2 中，以利于比較。

表2 用 5 種類型的協(xié)方差結(jié)構(gòu)模型擬合本資料的擬合效果比較

由上面的結(jié)果可知：采用 UN （無結(jié)構(gòu)型模型）協(xié)方差結(jié)構(gòu)模型擬合本資料，其迭代并沒有收斂，擬合效果不好。對比其他四種協(xié)方差結(jié)構(gòu)模型的擬合情況，可以看出采用CS協(xié)方差結(jié)構(gòu)模型擬合本資料效果較好，因為評價擬合效果的三個準(zhǔn)則中這個模型對應(yīng)的值較小。所以此時應(yīng)以采用CS 協(xié)方差結(jié)構(gòu)模型進(jìn)行混合模型分析的輸出結(jié)果為準(zhǔn)，其中F= 1267.75，P< 0.0001。說明不同時間點上測得的 HbH患者血紅蛋白平均含量之間的差異存在統(tǒng)計學(xué)意義。欲得出各時間點上 HbH 患者血紅蛋白含量均值之間兩兩比較的結(jié)果，可將原程序中所有過程步刪除，在第二個 DATA 步后換上以下過程步：

其中，“ddfm = sat”說明采用 Satterthwaite 近似的方法來計算分母的自由度，從而得到一個更加準(zhǔn)確的近似F值。但是，此法計算的結(jié)果有時會與 glm 過程計算的結(jié)果有所不同。

程序運(yùn)行后，得結(jié)果如下：

Least squares means

Differences of least squares means

這是各時間點上 HbH 患者血紅蛋白含量均值之間兩兩比較的結(jié)果，首先給出了各時間點 HbH 患者血紅蛋白含量均值與 0 比較的檢驗結(jié)果，沒有太大實際意義。后面給出了兩兩比較的結(jié)果，month 和 _month 列中的 0、1、2、3 分別代表時間的 4 個水平，即治療前、服藥 1 個月、服藥 2 個月、服藥 3 個月。讀者可查看最后兩列給出的兩兩比較的結(jié)果，可發(fā)現(xiàn)這 4 個時間點上 HbH 患者血紅蛋白含量均值之間的差異均有統(tǒng)計學(xué)意義，說明 4 個時間點上HbH 患者血紅蛋白平均含量各不相同。

[1]Hu LP.Application of statistical triple-type theory in the experiment design.Beijing: People’s Military Medical Press, 2006:107-120.(in Chinese)胡良平.統(tǒng)計學(xué)三型理論在實驗設(shè)計中的應(yīng)用.北京: 人民軍醫(yī)出版社, 2006:107-120.

[2]Hu LP.Scientific research design and statistical analysis of cardiovascular disease.Beijing: People’s Military Medical Press,2010:93-111.(in Chinese)胡良平.心血管病科研設(shè)計與統(tǒng)計分析.北京: 人民軍醫(yī)出版社,2010:93-111.