第9講協(xié)方差分析及混合線性模型

1、第第9 9講講 協(xié)方差分析與混合線協(xié)方差分析與混合線性模型性模型理學(xué)院理學(xué)院 汪曉銀汪曉銀 教授教授華中農(nóng)業(yè)大學(xué)數(shù)學(xué)建模創(chuàng)新實(shí)踐基地課件華中農(nóng)業(yè)大學(xué)數(shù)學(xué)建模創(chuàng)新實(shí)踐基地課件 如果在單因素、雙因素或多因素試驗(yàn)如果在單因素、雙因素或多因素試驗(yàn)中有無法控制的因素中有無法控制的因素x x影響試驗(yàn)的結(jié)果影響試驗(yàn)的結(jié)果Y Y，且且x x可以測量、可以測量、x x與與Y Y之間又有顯著的線性回之間又有顯著的線性回歸時(shí)，常常利用線性回歸來矯正歸時(shí)，常常利用線性回歸來矯正Y Y的觀測值、的觀測值、消去消去x x的差異對的差異對Y Y的影響。的影響。 例如，研究施肥對蘋果樹產(chǎn)量的影響，例如，研究施肥對蘋果樹產(chǎn)量的

2、影響，由于蘋果樹的長勢不齊，必須消去長勢對由于蘋果樹的長勢不齊，必須消去長勢對產(chǎn)量的影響。又如，研究飼料對動(dòng)物增重產(chǎn)量的影響。又如，研究飼料對動(dòng)物增重的影響，由于動(dòng)物的初重不同，必須消去的影響，由于動(dòng)物的初重不同，必須消去初重對增重的影響。初重對增重的影響。協(xié)方差分析協(xié)方差分析 這種不是在試驗(yàn)中控制某個(gè)因素，這種不是在試驗(yàn)中控制某個(gè)因素，而是在試驗(yàn)后對該因素的影響進(jìn)行估計(jì)，而是在試驗(yàn)后對該因素的影響進(jìn)行估計(jì)，并對試驗(yàn)指標(biāo)的值作出調(diào)整的方法稱為并對試驗(yàn)指標(biāo)的值作出調(diào)整的方法稱為統(tǒng)計(jì)控制，可以作為試驗(yàn)控制的輔助手統(tǒng)計(jì)控制，可以作為試驗(yàn)控制的輔助手段。以統(tǒng)計(jì)控制為目的，綜合線性回歸段。以統(tǒng)計(jì)控制為目

3、的，綜合線性回歸分析與方差分析所得到的統(tǒng)計(jì)分析方法，分析與方差分析所得到的統(tǒng)計(jì)分析方法，稱為協(xié)方差分析，所需要統(tǒng)計(jì)控制的一稱為協(xié)方差分析，所需要統(tǒng)計(jì)控制的一個(gè)或多個(gè)因素，例如蘋果樹的長勢，又個(gè)或多個(gè)因素，例如蘋果樹的長勢，又如動(dòng)物的初重等等稱為協(xié)變量。如動(dòng)物的初重等等稱為協(xié)變量。1 1 協(xié)方差分析思想原理協(xié)方差分析思想原理2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -理論理論2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -理論理論2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -理論理論2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -理論理論2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -計(jì)算計(jì)算dat

4、a ex; do a= ex; do a=1 to to 3;do i=;do i=1 to to 8; ;input x y ;output ;end;end;input x y ;output ;end;end;cards;cards;47 54 58 66 53 63 46 51 49 56 56 66 54 61 44 50 47 54 58 66 53 63 46 51 49 56 56 66 54 61 44 50 52 54 53 53 64 67 58 62 59 62 61 63 63 64 66 69 52 54 53 53 64 67 58 62 59 62 61 63

5、63 64 66 69 44 52 48 58 46 54 50 61 59 70 57 64 58 69 53 6644 52 48 58 46 54 50 61 59 70 57 64 58 69 53 66 ; ; proc glm;class a;model y=x a/solution;class a;model y=x a/solution;lsmeans a/stderr pdiff;lsmeans a/stderr pdiff;run; ;2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -計(jì)算計(jì)算2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -計(jì)算計(jì)算 施用三種肥料的產(chǎn)量矯正后有極顯

6、施用三種肥料的產(chǎn)量矯正后有極顯著的差異著的差異 2 2 單因素協(xié)方差分析單因素協(xié)方差分析- -計(jì)算計(jì)算3 3雙因素協(xié)方差分析雙因素協(xié)方差分析- -不考慮交互作用不考慮交互作用 方差來源方差來源平方和平方和自由度自由度 均方和均方和 F值值 顯著性顯著性 A QA r-1 MQA FA B QB s-1 MQB FB 誤差誤差 QErs-r-s MQE 總和總和 QT rs-23 3雙因素協(xié)方差分析雙因素協(xié)方差分析- -不考慮交互作用不考慮交互作用 data ex;do a=1 to 3 ;do b=1 to 5 ; data ex;do a=1 to 3 ;do b=1 to 5 ;input

7、 x y ;output; end; end;input x y ;output; end; end;cards;cards;8 2.85 10 4.24 12 3.00 11 4.94 10 2.888 2.85 10 4.24 12 3.00 11 4.94 10 2.8810 3.14 12 4.50 7 2.75 12 5.84 10 4.0610 3.14 12 4.50 7 2.75 12 5.84 10 4.0612 3.88 10 3.86 9 2.82 10 4.94 9 2.8912 3.88 10 3.86 9 2.82 10 4.94 9 2.89; ;proc glm

8、;class a b ;model y=x a b/solution;proc glm;class a b ;model y=x a b/solution;lsmeans a b/stderr pdiff;run;lsmeans a b/stderr pdiff;run;3 3雙因素協(xié)方差分析雙因素協(xié)方差分析- -不考慮交互作用不考慮交互作用3 3雙因素協(xié)方差分析雙因素協(xié)方差分析- -不考慮交互作用不考慮交互作用 方差來方差來源源平方和平方和 自由自由度度 均方均方和和 F值值 顯著顯著性性 A0.6046 2 0.3023 2.49 N B7.1245 4 1.7811 14.66 * 誤差

9、誤差0.8502 7 0.1215 總和總和8.5793 1各小區(qū)的產(chǎn)量矯正后沒有顯著的差異，各品各小區(qū)的產(chǎn)量矯正后沒有顯著的差異，各品種的產(chǎn)量矯正后有極顯著的差異。種的產(chǎn)量矯正后有極顯著的差異。3 3雙因素協(xié)方差分析雙因素協(xié)方差分析- -不考慮交互作用不考慮交互作用 方差來方差來源源平方平方和和自由度自由度 均方和均方和 F值值 顯著顯著性性 A QA r-1 MQA FA B QB s-1 MQB FB AB QAB(r-1)(s-1) MQAB FAB 誤差誤差 QErs(m-1)-1 MQE 總和總和 QT rsm-24 4雙因素協(xié)方差分析雙因素協(xié)方差分析- -考慮交互作用考慮交互作用

10、4 4雙因素協(xié)方差分析雙因素協(xié)方差分析- -考慮交互作用考慮交互作用data ex; do a=1 to 4; do b=1 to 2;data ex; do a=1 to 4; do b=1 to 2;do i=1 to 2;input x y;output;end;end;end;do i=1 to 2;input x y;output;end;end;end;cards;cards;14.6 97.8 12.1 94.2 19.5 11.2 18.814.6 97.8 12.1 94.2 19.5 11.2 18.8110.1 1.6 100.3 12.9 98.5 18.5 119.4

11、 110.1 1.6 100.3 12.9 98.5 18.5 119.4 12.2 114.7 12.8 99.2 10.7 89.612.2 114.7 12.8 99.2 10.7 89.612.212.2122.2 16.9 105.3 12.0 102.1 12.4 103.8 122.2 16.9 105.3 12.0 102.1 12.4 103.8 16.4 117.2 17.2 117.916.4 117.2 17.2 117.9proc glm; class a b;model y=x a bproc glm; class a b;model y=x a ba a* *b/

12、solution;lsmeans a b/stderr pdiff;b/solution;lsmeans a b/stderr pdiff;run;run;4 4雙因素協(xié)方差分析雙因素協(xié)方差分析- -考慮交互作用考慮交互作用4 4雙因素協(xié)方差分析雙因素協(xié)方差分析- -考慮交互作用考慮交互作用 方差來源方差來源 平方和平方和自由度自由度 均方和均方和 F值值顯著性顯著性 A277.43485 392.4782866.51 * B 2.845259 3 2.845259 0.20 N AB 12.848100 1 4.282700 0.30 N 誤差誤差 99.441171 714.205882

13、A與與B的交互作用矯正后不顯著，促生長劑的交互作用矯正后不顯著，促生長劑之間的差異極顯著，試驗(yàn)批次間的差異不顯著之間的差異極顯著，試驗(yàn)批次間的差異不顯著 4 4雙因素協(xié)方差分析雙因素協(xié)方差分析- -考慮交互作用考慮交互作用3.混合線性模型混合線性模型通過一個(gè)例子講述混合線性模型的使用通過一個(gè)例子講述混合線性模型的使用艾滋病療法的評價(jià)艾滋病療法的評價(jià)艾滋病是當(dāng)前人類社會最嚴(yán)重的瘟疫之一，從艾滋病是當(dāng)前人類社會最嚴(yán)重的瘟疫之一，從1981年發(fā)現(xiàn)以年發(fā)現(xiàn)以來的來的20多年間，它已經(jīng)吞噬了近多年間，它已經(jīng)吞噬了近3000萬人的生命。萬人的生命。艾滋病的醫(yī)學(xué)全名為艾滋病的醫(yī)學(xué)全名為“獲得性免疫缺損綜合癥

14、獲得性免疫缺損綜合癥”，英文簡稱，英文簡稱AIDS，它是由艾滋病毒（醫(yī)學(xué)全名為，它是由艾滋病毒（醫(yī)學(xué)全名為“人體免疫缺損病毒人體免疫缺損病毒”, 英文簡稱英文簡稱HIV）引起的。這種病毒破壞人的免疫系統(tǒng)，使人體）引起的。這種病毒破壞人的免疫系統(tǒng)，使人體喪失抵抗各種疾病的能力，從而嚴(yán)重危害人的生命。人類免喪失抵抗各種疾病的能力，從而嚴(yán)重危害人的生命。人類免疫系統(tǒng)的疫系統(tǒng)的CD4細(xì)胞在抵御細(xì)胞在抵御HIV的入侵中起著重要作用，當(dāng)?shù)娜肭种衅鹬匾饔茫?dāng)CD4被被HIV感染而裂解時(shí)，其數(shù)量會急劇減少，感染而裂解時(shí)，其數(shù)量會急劇減少，HIV將迅速增加，將迅速增加，導(dǎo)致導(dǎo)致AIDS發(fā)作。發(fā)作。3.混合線

15、性模型混合線性模型艾滋病治療的目的，是盡量減少人體內(nèi)艾滋病治療的目的，是盡量減少人體內(nèi)HIV的的數(shù)量，同時(shí)產(chǎn)生更多的數(shù)量，同時(shí)產(chǎn)生更多的CD4，至少要有效地降，至少要有效地降低低CD4減少的速度，以提高人體免疫能力。減少的速度，以提高人體免疫能力。迄今為止人類還沒有找到能根治迄今為止人類還沒有找到能根治AIDS的療法，的療法，目前的一些目前的一些AIDS療法不僅對人體有副作用，療法不僅對人體有副作用，而且成本也很高。許多國家和醫(yī)療組織都在積而且成本也很高。許多國家和醫(yī)療組織都在積極試驗(yàn)、尋找更好的極試驗(yàn)、尋找更好的AIDS療法。療法。3.混合線性模型混合線性模型請你完成以下問題：請你完成以下問

16、題：（1）利用附件）利用附件1的數(shù)據(jù)，預(yù)測繼續(xù)治療的效果，或者確定的數(shù)據(jù)，預(yù)測繼續(xù)治療的效果，或者確定最佳治療終止時(shí)間（繼續(xù)治療指在測試終止后繼續(xù)服藥，最佳治療終止時(shí)間（繼續(xù)治療指在測試終止后繼續(xù)服藥，如果認(rèn)為繼續(xù)服藥效果不好，則可選擇提前終止治療）。如果認(rèn)為繼續(xù)服藥效果不好，則可選擇提前終止治療）。（2）利用附件利用附件2的數(shù)據(jù)，評價(jià)的數(shù)據(jù)，評價(jià)4種療法的優(yōu)劣（僅以種療法的優(yōu)劣（僅以CD4為標(biāo)準(zhǔn)），為標(biāo)準(zhǔn)），并對較優(yōu)的療法預(yù)測繼續(xù)治療的效果，或者確并對較優(yōu)的療法預(yù)測繼續(xù)治療的效果，或者確定最佳治療終止時(shí)間。定最佳治療終止時(shí)間。(3) 艾滋病藥品的主要供給商對不發(fā)達(dá)國家提供的藥品價(jià)艾滋病藥品的

17、主要供給商對不發(fā)達(dá)國家提供的藥品價(jià)格如下：格如下：600mg zidovudine 1.60美元，美元，400mg didanosine 0.85美元，美元，2.25 mg zalcitabine 1.85美元，美元，400 mg nevirapine 1.20美元。如果病人需要考慮美元。如果病人需要考慮4種療種療法的費(fèi)用，對（法的費(fèi)用，對（2）中的評價(jià)和預(yù)測（或者提前終止）有）中的評價(jià)和預(yù)測（或者提前終止）有什么改變。什么改變。3.混合線性模型混合線性模型ID 療法療法 年齡年齡 時(shí)間時(shí)間 Log(CD4 count+1) 1236.427103.1355 1236.42717.57143.

18、0445 1236.427115.57142.7726 1236.427123.57142.8332 1236.427132.57143.2189 1236.4271403.0445 2447.846703.0681 2447.846783.8918 2447.8467163.9703 2447.8467233.6109 2447.846730.71433.3322 2447.8467393.0910 3160.287503.7377 4336.596904.1190 4336.59697.14294.1109 4336.596916.14294.7095 3.混合線性模型混合線性模型1對對4

19、種療法的療效評價(jià)的分析種療法的療效評價(jià)的分析對題目所給的附件對題目所給的附件2的數(shù)據(jù)進(jìn)行分析可知，決定病人的數(shù)據(jù)進(jìn)行分析可知，決定病人的的CD4的濃度的因素有年齡，檢查的時(shí)刻，治療方的濃度的因素有年齡，檢查的時(shí)刻，治療方案這三個(gè)因素。因此我們將年齡分成案這三個(gè)因素。因此我們將年齡分成5類，檢查的時(shí)類，檢查的時(shí)刻分為刻分為4個(gè)時(shí)間段，治療方案有個(gè)時(shí)間段，治療方案有4種。而問題是以種。而問題是以CD4的濃度為標(biāo)準(zhǔn)來評價(jià)療效的優(yōu)劣，即的濃度為標(biāo)準(zhǔn)來評價(jià)療效的優(yōu)劣，即CD4的濃的濃度越大，那么療效越好。由于考慮題中所給的樣本度越大，那么療效越好。由于考慮題中所給的樣本有有6000多個(gè)（病人的個(gè)數(shù)多個(gè)（

20、病人的個(gè)數(shù)*各個(gè)病人檢查的次數(shù)），各個(gè)病人檢查的次數(shù)），因此我們考慮用因此我們考慮用MIXED（混合線性模型）。下面我（混合線性模型）。下面我們就混合線性模型的原理進(jìn)行說明。們就混合線性模型的原理進(jìn)行說明。3.混合線性模型混合線性模型混合線性模型過程是擬合許多不同數(shù)據(jù)的混合線性混合線性模型過程是擬合許多不同數(shù)據(jù)的混合線性模型，并利用所擬合的模型對數(shù)據(jù)進(jìn)行統(tǒng)計(jì)推斷。模型，并利用所擬合的模型對數(shù)據(jù)進(jìn)行統(tǒng)計(jì)推斷。首先混合線性模型的主要假設(shè)是數(shù)據(jù)服從正態(tài)分布，首先混合線性模型的主要假設(shè)是數(shù)據(jù)服從正態(tài)分布，由于本題所給的數(shù)據(jù)的樣本容量有由于本題所給的數(shù)據(jù)的樣本容量有6000多個(gè)，因此，多個(gè)，因此，我們可

21、認(rèn)為它服從正態(tài)分布。又由于正態(tài)分布的數(shù)我們可認(rèn)為它服從正態(tài)分布。又由于正態(tài)分布的數(shù)據(jù)可完全有均值和方差確定，因此一個(gè)混合線性模據(jù)可完全有均值和方差確定，因此一個(gè)混合線性模型是由兩個(gè)模型決定的。分別是均值模型和方差模型是由兩個(gè)模型決定的。分別是均值模型和方差模型。型。MIXED使用約束最大似然的方法來擬合數(shù)據(jù)的。使用約束最大似然的方法來擬合數(shù)據(jù)的。一旦數(shù)據(jù)的模型已經(jīng)建立，我們可以使用該模型通一旦數(shù)據(jù)的模型已經(jīng)建立，我們可以使用該模型通過固定效應(yīng)參數(shù)和協(xié)方差參數(shù)進(jìn)行統(tǒng)計(jì)推斷。用這過固定效應(yīng)參數(shù)和協(xié)方差參數(shù)進(jìn)行統(tǒng)計(jì)推斷。用這些統(tǒng)計(jì)量可以對模型進(jìn)行評價(jià)。些統(tǒng)計(jì)量可以對模型進(jìn)行評價(jià)。3.混合線性模型混合

22、線性模型再者，該分析的重要假設(shè)是數(shù)據(jù)是正態(tài)分布的，由再者，該分析的重要假設(shè)是數(shù)據(jù)是正態(tài)分布的，由于我們將附錄于我們將附錄2的數(shù)據(jù)進(jìn)行了分類。由于數(shù)據(jù)出現(xiàn)在的數(shù)據(jù)進(jìn)行了分類。由于數(shù)據(jù)出現(xiàn)在類（如可能是同一年齡段中），那可能的情況是來類（如可能是同一年齡段中），那可能的情況是來自同一年齡段的這些自同一年齡段的這些CD4的值是相關(guān)的，不是獨(dú)立的值是相關(guān)的，不是獨(dú)立的。鑒于此，由于附錄的。鑒于此，由于附錄2所給的數(shù)據(jù)是高度數(shù)據(jù)（即所給的數(shù)據(jù)是高度數(shù)據(jù)（即樣本容量很大），那么考慮這種因素是相關(guān)。因此樣本容量很大），那么考慮這種因素是相關(guān)。因此我們對這些因素進(jìn)行相關(guān)性建模，我們使用隨機(jī)效我們對這些因素進(jìn)行

23、相關(guān)性建模，我們使用隨機(jī)效用。在本題中，我們規(guī)定年齡、檢查時(shí)刻為隨機(jī)效用。在本題中，我們規(guī)定年齡、檢查時(shí)刻為隨機(jī)效應(yīng)，即使得具有相同年齡水平或相同檢查時(shí)刻的水應(yīng)，即使得具有相同年齡水平或相同檢查時(shí)刻的水平之間存在共同的相關(guān)性，那么，此模型才較為合平之間存在共同的相關(guān)性，那么，此模型才較為合理。理。3.混合線性模型混合線性模型2 對對4種療法的療效模型的建立種療法的療效模型的建立1）數(shù)據(jù)的處理）數(shù)據(jù)的處理a、所有病人的年齡是在、所有病人的年齡是在14.9021，74.193的區(qū)間內(nèi)的區(qū)間內(nèi),以（以（74.193-14.9021）/5=11.85818為區(qū)間長度。我為區(qū)間長度。我們將患者按照年齡階

24、段分為們將患者按照年齡階段分為5級，分別記作級，分別記作1至至5（如（如表表4所示）；所示）；b、我們將患者所接受的治療方法分為、我們將患者所接受的治療方法分為4種，分別記種，分別記作作1至至4；c、我們將病人的檢查的時(shí)刻、我們將病人的檢查的時(shí)刻0，40以以10為區(qū)間長為區(qū)間長度，分為度，分為4級，分別記作級，分別記作1至至4 3.混合線性模型混合線性模型3.混合線性模型混合線性模型2）確定固定效應(yīng)和隨機(jī)效應(yīng)）確定固定效應(yīng)和隨機(jī)效應(yīng)固定效用是設(shè)計(jì)者所研究的因素，在此題中為固定效用是設(shè)計(jì)者所研究的因素，在此題中為4種療法。因?yàn)椴∪说哪挲g是隨機(jī)的，并且病種療法。因?yàn)椴∪说哪挲g是隨機(jī)的，并且病人接受檢驗(yàn)的時(shí)刻是也是隨機(jī)進(jìn)行的，故病人人接受檢驗(yàn)的時(shí)刻是也是隨機(jī)進(jìn)行的，故病人的年齡和病人接受檢驗(yàn)的治療時(shí)刻應(yīng)該為隨機(jī)的年齡和病人接受檢驗(yàn)的治療時(shí)刻應(yīng)該為隨機(jī)效應(yīng)。效應(yīng)。3.混合線性模型混合線性模型3）協(xié)方差結(jié)構(gòu)的選擇）協(xié)方差結(jié)構(gòu)的選擇模型中具體選用哪種結(jié)構(gòu)矩陣：在相同模型結(jié)構(gòu)下，模型中具體選用哪種結(jié)構(gòu)矩陣：在相同模型結(jié)構(gòu)下，選擇幾個(gè)不同結(jié)構(gòu)的協(xié)方差矩陣，從中選取似然比選擇幾個(gè)不同結(jié)構(gòu)的協(xié)方差矩陣，從中選取似然比統(tǒng)計(jì)量（統(tǒng)計(jì)量（-2Log Likeli-hood）、）、Akaikes Information Schwartz Bayesian三個(gè)指標(biāo)均較小三個(gè)指標(biāo)