方差分析二重復(fù)測量資料方差分析沉毅

1、方差分析（二）：重復(fù)測量資料的方差分析沈 毅2005.31 重復(fù)測量（repeated measure）是指對同一觀察對象的同一觀察指標(biāo)在不同時間點上進(jìn)行的多次測量，用于分析該觀察指標(biāo)在不同時間上的變化特點。這類測量資料在臨床和流行病學(xué)研究中比較常見。例如，為研究某種藥物對高血壓病人的治療效果，需要定時多次測定受試者的血壓，以分析其血壓的變動情況。一些傳統(tǒng)的統(tǒng)計方法，如t檢驗、方差分析、線性回歸模型等，都要求各次觀察是相互獨立的。而重復(fù)測量資料由于是對同一受試者的某項觀察指標(biāo)進(jìn)行的多次測量，在同一受試者的多次測量之間可能存在某種相關(guān)性，用通常的統(tǒng)計方法就不能充分揭示出其內(nèi)在的特點，有時甚至?xí)?/p>

2、出錯誤的結(jié)論。因此，有關(guān)重復(fù)測量資料的分析方法是近代統(tǒng)計學(xué)研究的熱點之一。2重復(fù)測量數(shù)據(jù)結(jié)構(gòu)與獨立數(shù)據(jù)結(jié)構(gòu)的區(qū)別及其優(yōu)缺點 為了解重復(fù)測量數(shù)據(jù)與獨立觀察數(shù)據(jù)之間的區(qū)別，下面用一個完全隨機(jī)設(shè)計的獨立數(shù) 據(jù)結(jié)構(gòu)與具有一個受試者內(nèi)因素(時間)的重復(fù)測量數(shù)據(jù)結(jié)構(gòu)進(jìn)行比較。 完全隨機(jī)設(shè)計的獨立數(shù)據(jù)結(jié)構(gòu)：從正常人、可疑硅沉著病者及一期硅沉著病病人中各隨機(jī)抽取5人，測量他們的血清黏蛋白含量(mgL)，結(jié)果列于下表中第一部分。 重復(fù)測量的數(shù)據(jù)結(jié)構(gòu)：對5名粉塵作業(yè)工人的血清黏蛋白含量(mg幾)連續(xù)3年的測量 結(jié)果列于下表中第二部分。3 從上表的例子可以看出，獨立數(shù)據(jù)結(jié)構(gòu)的各個觀察值是彼此獨立的，它適宜于用通常的

3、方差分析方法做統(tǒng)計分析。重復(fù)測量數(shù)據(jù)結(jié)構(gòu)是對每一受試者的同一觀察指標(biāo)(血清黏蛋白含量)進(jìn)行的多次測量。由于這種多次測量之間可能存在相關(guān)性，就需要用特殊的統(tǒng)計方法進(jìn)行分析。 重復(fù)測量設(shè)計的主要優(yōu)點是可以減少樣本含量，其次是能夠控制個體變異，即個體差異。例如在單因素實驗中，可以用隨機(jī)區(qū)組(或稱配伍組)設(shè)計方法來縮小隨機(jī)誤差。而重復(fù)測量設(shè)計是以同一個受試者作為一個區(qū)組，故可以把它看成為是隨機(jī)區(qū)組設(shè)計的一種極端形式。但在隨機(jī)區(qū)組設(shè)計下的每一測量都是在不同受試者身上進(jìn)行的，它們對某種處理因素的反應(yīng)是獨立的，符合獨立性的假定。而在重復(fù)測量設(shè)計下的測量是在同一受試者身上進(jìn)行的，它們對同一處理因素在不同時間上

4、的反應(yīng)可能是不獨立的，后一次的測量結(jié)果可能受到前一次測量結(jié)果的影響。因此，對同一個體在不同時間上的測量值之間就可能存在相關(guān)關(guān)系。這給分析工作帶來了一定的復(fù)雜性。 在實際工作中，重復(fù)測量資料比獨立觀察資料往往更為多見。如在臨床研究中，需要觀察病人在不同時間的某些生理、生化或病理指標(biāo)的變化趨勢，不同時間或療程的治療效果。在流行病學(xué)研究中觀察隊列人群在不同時間上的發(fā)病情況。在衛(wèi)生學(xué)研究中，縱向觀察兒童的生長發(fā)育規(guī)律等。重復(fù)測量資料在自然科學(xué)和社會科學(xué)的很多領(lǐng)域內(nèi)都可見到。4 第一節(jié) 重復(fù)測量資料方差分析對協(xié)方差陣的要求 在對重復(fù)測量資料進(jìn)行方差分析時，除要求樣本是隨機(jī)的、在處理的同一個水平上的觀察是

5、獨立的以及每一水平上的測定值都來自正態(tài)總體外，特別強(qiáng)調(diào)協(xié)方差陣（covariance matrix）的球形性（sphericity）或為符合對稱性（compound symmetry）。Box（1954）指出，若球形性質(zhì)得不到滿足，則方差分析的F值是有偏的，這會造成過多的拒絕本來是真的無效假設(shè)(增加型錯誤)。 5一、協(xié)方差陣的球形性檢驗 方差是指在某一時點上測定值變異性的大小，而協(xié)方差是指在兩個不同時點上測定值相互變異性的大小。如果在某個時點上的取值不影響其他時點上的取值，則協(xié)方差為0，反之，則不為0。由方差協(xié)方差構(gòu)成的矩陣稱協(xié)方差陣。 6 設(shè) k、l為兩個測定時點， 代表協(xié)方差陣中的元素。當(dāng)

6、kl時為方差，kl時為協(xié)方差。共有a個測定時點，將這a個方差和（a-1）/2個協(xié)方差排成協(xié)方差陣V為：7 協(xié)方差陣的球形性質(zhì)是指該矩陣主對角線元素（方差）相等、非主對角線元素（協(xié)方差）為零。用Mauchly氏法檢驗協(xié)方差陣的球形性質(zhì)。Mauchly氏檢驗的 P值若大于研究者所選擇的顯著性水準(zhǔn)時，說明協(xié)方差陣的球形性質(zhì)得到滿足。否則，必須對與時間有關(guān)的F統(tǒng)計量的分子、分母自由度進(jìn)行調(diào)整，以便減少犯類錯誤的概率。調(diào)整系數(shù)為（epsilon）。 8二、自由度調(diào)整方法 1. 調(diào)整系數(shù)的計算 有兩種調(diào)整系數(shù)。 （1）Greenhouse-Geisser調(diào)整系數(shù) 為： 式（10-2）中的 是矩陣（10-1

7、）中第k行第l列元素， 是所有元素的總平均值， 是主對角線元素的平均值， 是第k行的平均值。 的取值在1.0與1/（a-1）之間。 9 （2）HuynhFeldt調(diào)整系數(shù) 研究表明，當(dāng)真值在0.7以上時，用 進(jìn)行自由度凋整后的統(tǒng)計學(xué)結(jié)論偏于保守，故Huynh和Feldt提出用平均調(diào)整值 值進(jìn)行調(diào)整。 值的計算公式為 式（10-3）中的g是對受試對象的某種特征（如性別或年齡）進(jìn)行分組的組數(shù)，n是每組的觀察例數(shù)。 10 2調(diào)整規(guī)則 只對具有重復(fù)測定性質(zhì)的時間效應(yīng)的F值的自由度，和處理時間交互作用的F值的自由度進(jìn)行調(diào)整。由于F值有兩個自由度v1和v2，調(diào)整的分子自由度 ，分母自由度 具體計算時可用

8、或 代替。用調(diào)整所得的 及 的值查 F界值表，得 。由于1，所以調(diào)整后的F臨界值要大于調(diào)整前的F臨界值。 11第二節(jié) 單因素重復(fù)測量資料的方差分析 一、單因素重復(fù)測量資料的例子12 例10-1 一項關(guān)于不同藥物治療心律失常效果的對比研究。對9例經(jīng)常出現(xiàn)心室早搏的病人于用藥前測定其心率后進(jìn)行隨機(jī)化給藥。一部分病人按A藥安慰劑（C藥）B藥的順序給藥，另一部分病人按B藥安慰劑（C藥）A藥的順序給藥。安慰劑（C藥）持續(xù)一周作為藥物后效的清除期。比較用藥前與各種藥物及A藥與B藥之間的心律差別。表10-1列出9名受試病人在用藥前、安慰劑（C藥）期及用藥（A與B）期的心率。 上表的下部及右側(cè)所列的數(shù)據(jù)均從表

9、體部分的原始測量值派生而來，是為方差分析而準(zhǔn)備的。本例藥物水平數(shù)a4，每組觀察例數(shù)n9，觀察值總個數(shù)Nan36。 13二、方差分析的步驟 1提出檢驗假設(shè) 檢驗假設(shè)為： H0：1234； H1：ih，至少有一個不等式成立。14 2計算離均差平方和、自由度及均方 有總離均差平方和、處理因素離均差平方和、受試對象間離均差平方和及受試對象內(nèi)離均差平方和等。計算公式為： （1）總離均差平方和SS總及總自由度v總的計算 本例，SS總201647-（2665）2/（49）4362.97，v總36-135。 （2）處理因素的離均差平方和SS處理及自由度v處理的計算15 （3）受試對象間離均差平方和SS對象間及

10、自由度v對象間的計算 （4）受試對象內(nèi)離均差平方和SS對象內(nèi)及自由度v對象內(nèi)的計算16 （5）誤差的離均差平方和SS誤差與自由度v誤差的計算 SS誤差SS總-SS處理-SS對象間，v誤差（n-1）（a-1） （10-8） 本例，SS誤差4362.97-1185.42-2023.721153.83，v誤差（9-1）（4-1）24 根據(jù)以上4種離均差平方和與自由度計算所得的均方見表 10-2。 3. 計算F值 由于主要是處理因素的統(tǒng)計學(xué)檢驗，故只計算處理因素的F值。 F處理MS處理/MS誤差，F(xiàn)處理服從v1v處理與v2v誤差的F分布 本例，F(xiàn)處理395.14/48.088.22；v13，v224。

11、所有計算結(jié)果均列于表10-2中。 17 SAS程序給出本例的協(xié)方差陣Mauchly球形性檢驗的結(jié)果為P0.1628，故不必進(jìn)行自由度的調(diào)整。查F界值表得：F0.01（3，24）3.01，F(xiàn)0.01（3，24）4.72。本例處理因素的F值為8.22，大于F0.01（3，24），故拒絕無效假設(shè)，說明處理因素間的差別具有統(tǒng)計學(xué)意義。 由表10-2看出，受試對象內(nèi)離均差平方和等于處理因素的離均差平方和與誤差的離均差平方和兩項之和。 SAS程序給出 。用 調(diào)整的處理因素的分子自由度為 ；分母自由度為 。查F界值表得調(diào)整自由度后的F臨界值F0.05（2，19）3.52，比未調(diào)整的F臨界值大。未調(diào)整的概率

12、P0.0006， 調(diào)整概率P0.0020。 18三、平均值之間的多重比較 以上用單因素重復(fù)測量方差分析方法對心率資料進(jìn)行分析之后所得到的統(tǒng)計學(xué)結(jié)論是：拒絕無效假設(shè)，即在治療藥物的四個水平中，至少有一個水平的總體平均值不同于其他水平的總體平均值。為了確定這個特殊總體，必須進(jìn)行平均值之間的多重比較。但此處不能采用上一講中介紹的多重比較方法，因為那些方法都是建立在獨立樣本基礎(chǔ)之上的。這里可采用配對樣本的差值t檢驗，因為配對樣本就是重復(fù)測量試驗中一種最簡單的對比研究設(shè)計。其檢驗步驟如下： 19 1. 計算每一個病人在不同給藥情況下的差值：di（j-h）Yij-Yih，i為病人號，j，h為藥物水平號，各

13、di（j-h）列于表10-3的右邊各列中。由于設(shè)計時只考慮用藥前與各種藥物及A藥與B藥之間的差縱情況，表中只列出了di（1-2）、di（1-3）、di（1-4）及di（2-4）四種組合，而不是所有可能的六種組合。 2. 進(jìn)行差值t檢驗（參閱成對測量值之差的比較），可分別得到t值為：t（1：2）4.41，t（1：3）0.03，t（1：4）3.19，t（2：4）-0.96 3計算校正臨界值t 由于是對同一份資料進(jìn)行多重比較，為克服累積類錯誤對結(jié)果判斷所造成的影響，根據(jù)Bonferroni不等式的原理對臨界t值進(jìn)行調(diào)整。 20 首先確定比較的次數(shù) c。因該研究已事先確定只作 4次比較，故c4。若在方

14、差分析之后再作多重比較，則只能取所有可能的比較次數(shù)。例如本例在方差分析之后再進(jìn)行比較時，則比較的次數(shù)應(yīng)為 c4（41）/26。 其次是選擇累積 類錯誤的概率0.10。采用雙側(cè)檢驗，每次檢驗所用的I類錯誤概率水準(zhǔn)為0.10/40.0125，自由度vn-19-18，查附表2：t界值表得t0.0125/2（8）3.21。與前面計算出的t值相比較，可見用藥前心率與服用A藥后心率之差具有統(tǒng)計學(xué)意義。用藥后心率平均降低12.44次/分，而用藥前心率與服安慰劑后心率之間以及A藥與B藥之間心率之差無統(tǒng)計學(xué)意義。用藥前心率與用B藥后心率之差接近顯著性水平。 21第三節(jié) 兩因素重復(fù)測定資料的方差分析 兩因素重復(fù)測

15、定資料中的兩因素是指一個組間因素（處理因素）和一個組內(nèi)因素（時間因素）。組間因素是指分組或分類變量，它把所有受試對象按分類變量的水平分為幾個組。組內(nèi)因素是指重復(fù)測定的時間變量。例10-1只有組內(nèi)因素，沒有組間因素。 22一、例子 例10-2 一項藥物代謝動力學(xué)研究，目的是對比某種藥物的不同劑型在體內(nèi)的代謝速度。劑型分膠囊型和片劑型。將16名受試對象隨機(jī)分為兩組，每組8名。一組給予膠囊，另一組給予片劑，分別在服藥后1、2、4、6及8小時測定血中的藥物濃度。測定結(jié)果見表10-4。 S（9.73）2+（54.61）2+（30.30）2+（70.20）2238012.38 本例的組間因素是藥物劑型，組

16、內(nèi)因素是測定時間。各下標(biāo)的意義是：i（i1，2，g）為組間因素的分組號，j（j1，2，p）為測定時間點的序號，k（k1，2，ni）為組間因素第i水平的受試對象號，受試對象總數(shù)為n1n2ng。當(dāng)各ni相等時，則用 n代替ni。測量值總個數(shù)Ngnp本例g2；各組受試對象數(shù)n8，p5，受試對象總數(shù)為2816例，測量值總個數(shù)N80。23 二、方差分析模型 一個組間因素，一個組內(nèi)因素的方差分析模型為： Yijk+i+i+（）ij+（i）k+ijk 模型中各參數(shù)的意義是：為總體平均值；i為處理組i的效應(yīng)；i為第j個測定時間點的效應(yīng)；（）ij為第i組在第j個測定時點上的效應(yīng)，屬交互作用，為固定效應(yīng)；（i）k

17、為第i組第k個觀察對象的效應(yīng)，屬隨機(jī)效應(yīng)；ijk為誤差項。給定限制條件為： 模型中的參數(shù)估計值與平均值之間的關(guān)系見 表 10-5。24三、方差分析的步驟 1.離均差平方和、自由度及均方的計算 令 為觀察值總和， 為觀察值平方總和， 為第i組觀察值之和， 第j時點觀察值之和， 為第i組第k個受試對象的觀察值之和， 為在（ij）水平上的觀察值之和。觀察值總個數(shù)為Ngpn。表10-4中同時列出了Tij，Tik，Ti，Tj，T和S的值。將各離均差平方和、自由度及均方的計算公式列于表 10-6中。 用表10-6中所列的計算公式計算出各離均差平方和、自由度及均方見表10-7。 2計算F值本例的各F值見表1

18、0-7。25 3確定P值并作出統(tǒng)計推斷 在查閱F臨界值表前，要對組內(nèi)變量（時間）及組間 組內(nèi)（劑型時間）交互作用的自由度作調(diào)整。調(diào)整時要應(yīng)用 的估計值。G-G調(diào)整系數(shù) ； H-F調(diào)整系數(shù) 。如用 計算的F組內(nèi)自由度調(diào)整值為： F組間*組內(nèi)的分子自由度的調(diào)整值為： 對組間（劑型）效應(yīng)仍用原來的自由度 查F臨界值表得：F0.05（1，14）4.60。對組內(nèi)（時間）效應(yīng)及組間組內(nèi)交互效應(yīng)的F臨界值要用調(diào)整自由度得F0.05（2，29）3.33。 統(tǒng)計學(xué)推斷為：不同劑型藥物的血液濃度間的差別無統(tǒng)計學(xué)意義（ P0.0645），不同時間的血藥濃度間的差別具有統(tǒng)計學(xué)意義（ P0.0001），劑型與時間之間無

19、明顯的交互作用（ P0.3413）。26第四節(jié) 趨勢分析 趨勢分析（trend analysis）用于分析某項觀察指標(biāo)在整個時間區(qū)間內(nèi)的變化趨勢，或者對比不同處理組在趨勢上的差異。例如分析不同類型糖尿病病人在一天中的血糖變化規(guī)律，分析不同性別兒童的生長發(fā)育趨勢等。本節(jié)所介紹的趨勢分析方法是把應(yīng)變量作為時間的函數(shù)來配合多項式曲線，不要求資料的協(xié)方差陣滿足球形性條件。 27 一、多項式曲線模型 一般的多項式回歸模型可表示為： 本節(jié)所介紹的多項式曲線是采用正交多項式系數(shù)來配合模型。其模型結(jié)構(gòu)為： 模型中的 是在（r）階下對第j時點平均值 的預(yù)報值，Cj(1)，Cj(r)是從第（1）階起直到（r）階的

20、第j個正交多項式系數(shù)。0，1，r為模型參數(shù)，從實際資料中求得它們的估計值，b0，b1，br。設(shè)測定的時間點有p個，從理論上講，冪的階次可以達(dá)到rp。在這種情況下，理論模型與觀察資料完全一致。但在實際工作中，一般取冪的階次 r（p-1）。當(dāng) r1時為線形模型。利用多項式曲線作趨勢分析是將兩個組（如治療組與對照組或不同的治療組）之間的比較轉(zhuǎn)換為兩個組參數(shù)之間的比較。如果任何一對參數(shù)之間的差別具有統(tǒng)計學(xué)意義，就可認(rèn)為這兩條曲線的趨勢不一致。 28二、例 子 例10-3 一項藥物代謝動力學(xué)研究，目的是比較某種藥物的兩種劑型（膠囊型和片劑型）在體內(nèi)的代謝過程。每組8名受試對象，給藥后等時距地每隔1小時測

21、定一次血藥濃度（g/ml，用曲線下面積表示），共測定5次?？倻y定值的個數(shù)N40。資料見表10-8。 為了解用正交多項式作趨勢分析的過程，先對膠囊型組進(jìn)行單組資料的趨勢分析，然后再把片劑型資料包括進(jìn)來，作兩組趨勢的對比分析。 29三、單組資料的趨勢分析 用表10-8中膠囊型資料作單組資料趨勢分析。計算所需的數(shù)據(jù)經(jīng)整理后列于表10-9。由于在模型配合過程中需要反復(fù)使用誤差均方進(jìn)行假設(shè)檢驗，故首先進(jìn)行單向方差分析，得到誤差均方的估計值為：MS誤差899.98，v誤差35。 趨勢分析的過程可概括為六個步驟，各步驟的計算結(jié)果見表 10-11。 1.從正交多項式系數(shù)表10-10中選擇一組系數(shù)C（r）（c1

22、（r） c2（r）cp（r）。系數(shù)C（r）的下標(biāo)（r）是多項式曲線的階數(shù)。本例 p5。從表中水平數(shù)為 5的四種類型的多項式中選擇一種，得到一組多項式系 數(shù)C（r）。配合過程通常先由低階向高階逐步過渡，故首先選擇一組線性多項式系數(shù)為C（1）（-2 -1 0 1 2）。第二步配合二次曲線時選擇一組系數(shù)為C（2）（2 -1 -2 -1 2），余類推。 30 2計算r次冪下的（r）值： 3計算離均差平方和SS（r） 式（10-11）的分子中的n為觀察例數(shù)，分母為系數(shù)的平方和（見表 10-10 的最后一列）。 4計算F值，確定概率P并作出統(tǒng)計推斷 F值的計算公式為 統(tǒng)計量F服從v11及v2=v誤差的F分

23、布 5計算該次冪下的回歸系數(shù)r的估計值br 31 6確定回歸方程并求出預(yù)報值回歸方程為； 本例資料一直配合到四次冪的多項式曲線時仍然有統(tǒng)計學(xué)意義，說明該回歸曲線的彎曲度較大。最后得到的回歸方程如下： 用式（10-15）對本例資料的預(yù)報結(jié)果為： 同理得到 利用測量的平均值和預(yù)報的平均值制作的趨勢圖（ 圖 10-1）中，可見兩條曲線幾乎完全重疊。 經(jīng)配合多項式曲線后所得的總的回歸平方和 SS回歸為各次回歸平方和的合計，即： SS回歸113028.67+14366.44+10841.03+4005.44142241.53 即四次多項式曲線可以解釋總變異（ SS總173734.51）的約 82。32四、兩組重復(fù)測定資料的趨勢對比 對例10-3的實