R語(yǔ)言實(shí)戰(zhàn)topic7方差分析

1、Topic7 方差分析一、相關(guān)術(shù)語(yǔ)以焦慮癥治療為例，現(xiàn)有兩種治療方案：認(rèn)知行為療法（CBT）和眼動(dòng)脫敏再加工法（EMDR）。我們招募10位焦慮癥患者作為志愿者，隨機(jī)分配一半的人接受為期五周的CBT，另外一半接受為期五周的EMDR，設(shè)計(jì)方案如表9-1所示。在治療結(jié)束時(shí)，要求每位患者都填寫(xiě)狀態(tài)特質(zhì)焦慮問(wèn)卷（STAI），也就是一份焦慮度測(cè)量的自我評(píng)測(cè)報(bào)告。在這個(gè)實(shí)驗(yàn)設(shè)計(jì)中，治療方案是兩水平（CBT、 EMDR）的組間因子。之所以稱(chēng)其為組間因子，是因?yàn)槊课换颊叨純H被分配到一個(gè)組別中，沒(méi)有患者同時(shí)接受CBT和EMDR。表中字母s代表受試者（患者）。 STAI是因變量，治療方案是自變量。由于在每種治療方案

2、下觀(guān)測(cè)數(shù)相等，因此這種設(shè)計(jì)也稱(chēng)為均衡設(shè)計(jì)（balanced design）；若觀(guān)測(cè)數(shù)不同，則稱(chēng)作非均衡設(shè)計(jì)（unbalanceddesign）。因?yàn)閮H有一個(gè)類(lèi)別型變量，表9-1的統(tǒng)計(jì)設(shè)計(jì)又稱(chēng)為單因素方差分析（one-way ANOVA），或進(jìn)一步稱(chēng)為單因素組間方差分析。方差分析主要通過(guò)F檢驗(yàn)來(lái)進(jìn)行效果評(píng)測(cè)，若治療方案的F檢驗(yàn)顯著，則說(shuō)明五周后兩種療法的STAI得分均值不同。假設(shè)你只對(duì)CBT的效果感興趣，則需將10個(gè)患者都放在CBT組中，然后在治療五周和六個(gè)月后分別評(píng)價(jià)療效，設(shè)計(jì)方案如表9-2所示。療法（therapy）和時(shí)間（time）都作為因子時(shí)，我們既可分析療法的影響（時(shí)間跨度上的平均）和

3、時(shí)間的影響（療法類(lèi)型跨度上的平均），又可分析療法和時(shí)間的交互影響。前兩個(gè)稱(chēng)作主效應(yīng)，交互部分稱(chēng)作交互效應(yīng)。當(dāng)設(shè)計(jì)包含兩個(gè)甚至更多的因子時(shí)，便是多因素方差分析設(shè)計(jì)，比如兩因子時(shí)稱(chēng)作雙因素方差分析，三因子時(shí)稱(chēng)作三因素方差分析，以此類(lèi)推。若因子設(shè)計(jì)包括組內(nèi)和組間因子，又稱(chēng)作混合模型方差分析，當(dāng)前的例子就是典型的雙因素混合模型方差分析。本例中，你將做三次F檢驗(yàn)：療法因素一次，時(shí)間因素一次，兩者交互因素一次。若療法結(jié)果顯著，說(shuō)明CBT和EMDR對(duì)焦慮癥的治療效果不同；若時(shí)間結(jié)果顯著，說(shuō)明焦慮度從五周到六個(gè)月發(fā)生了變化；若兩者交互效應(yīng)顯著，說(shuō)明兩種療法隨著時(shí)間變化對(duì)焦慮癥治療影響不同（也就是說(shuō)，焦慮度從五

4、周到六個(gè)月的改變程度在兩種療法間是不同的）。現(xiàn)在，我們對(duì)上面的實(shí)驗(yàn)設(shè)計(jì)稍微做些擴(kuò)展。眾所周知，抑郁癥對(duì)病癥治療有影響，而且抑郁癥和焦慮癥常常同時(shí)出現(xiàn)。即使受試者被隨機(jī)分配到不同的治療方案中，在研究開(kāi)始時(shí)，兩組療法中的患者抑郁水平就可能不同，任何治療后的差異都有可能是最初的抑郁水平不同導(dǎo)致的，而不是由于實(shí)驗(yàn)的操作問(wèn)題。抑郁癥也可以解釋因變量的組間差異，因此它常稱(chēng)為混淆因素（confounding factor）。由于你對(duì)抑郁癥不感興趣，它也被稱(chēng)作干擾變數(shù)（nuisance variable）。假設(shè)招募患者時(shí)使用抑郁癥的自我評(píng)測(cè)報(bào)告，比如白氏抑郁癥量表（BDI），記錄了他們的抑郁水平，那么你可以在

5、評(píng)測(cè)療法類(lèi)型的影響前，對(duì)任何抑郁水平的組間差異進(jìn)行統(tǒng)計(jì)性調(diào)整。本案例中， BDI為協(xié)變量，該設(shè)計(jì)為協(xié)方差分析（ANCOVA）。以上設(shè)計(jì)只記錄了單個(gè)因變量情況（STAI），為增強(qiáng)研究的有效性，可以對(duì)焦慮癥進(jìn)行其他的測(cè)量（比如家庭評(píng)分、醫(yī)師評(píng)分，以及焦慮癥對(duì)日常行為的影響評(píng)價(jià)）。當(dāng)因變量不止一個(gè)時(shí)，設(shè) 計(jì) 被 稱(chēng) 作 多 元 方 差 分 析 （MANOVA ）， 若 協(xié) 變 量 也 存 在 ， 那 么 就 叫 多 元 協(xié) 方 差 分析（MANCOVA）。二、單因素方差分析（1）單因素方差分析以multcomp包中的cholesterol數(shù)據(jù)集為例（取自Westfall、 Tobia、 Rom、 H

6、ochberg， 1999）， 50個(gè)患者均接受降低膽固醇藥物治療（trt）五種療法中的一種療法。其中三種治療條件使用藥物相同，分別是20mg一天一次（1time）、 10mg一天兩次（2times）和5mg一天四次（4times）。剩下的兩種方式（drugD和drugE）代表候選藥物。哪種藥物療法降低膽固醇（響應(yīng)變量）最多呢？#單因素方差分析（感興趣地是比較分類(lèi)因子定義的兩個(gè)或多個(gè)組別中的因變量均值） install.packages("multcomp") library(multcomp) attach(cholesterol) str(cholesterol) ch

7、olesterol table(trt) aggregate(response,by=list(trt),FUN=mean) aggregate(response,by=list(trt),FUN=sd) fit<-aov(responsetrt) summary(fit) library(gplots) plotmeans(responsetrt,xlab="Treatment",ylab="Response",main="Mean Plotn with 95%CI")（2）多重比較雖然ANOVA對(duì)各療法的F檢驗(yàn)表明五種藥物療

8、法效果不同，但是并沒(méi)有告訴你哪種療法與其他療法不同。多重比較可以解決這個(gè)問(wèn)題。TukeyHSD(fit) par(las=2) par(mar=c(5,4,6,2) plot(TukeyHSD(fit) （3）評(píng)估檢驗(yàn)的假設(shè)條件當(dāng)因變量服從正態(tài)分布，各組方差相等時(shí)，可用Q-Q圖來(lái)檢驗(yàn)正態(tài)性假設(shè)qqPlot（）要求用lm（）擬合，若數(shù)據(jù)落在95%的置信區(qū)間范圍內(nèi)，說(shuō)明滿(mǎn)足正態(tài)性假設(shè)。library(car)qqPlot(lm(response trt, data=cholesterol), simulate=TRUE, main="Q-Q Plot", labels=FALS

9、E)R提供了一些可用來(lái)做方差齊性檢驗(yàn)的函數(shù)bartlett.test(response trt, data=cholesterol)離群點(diǎn)檢驗(yàn)library(car)outlierTest(fit)三、單因素協(xié)方差分析單因素協(xié)方差分析（ANCOVA）擴(kuò)展了單因素方差分析（ANOVA），包含一個(gè)或多個(gè)定量的協(xié)變量。下面的例子來(lái)自于multcomp包中的litter數(shù)據(jù)集（見(jiàn)Westfall et al.， 1999）。懷孕小鼠被分為四個(gè)小組，每個(gè)小組接受不同劑量（0、 5、 50或500）的藥物處理。產(chǎn)下幼崽的體重均值為因變量，懷孕時(shí)間為協(xié)變量。data(litter, package=&quo

10、t;multcomp")attach(litter)table(dose) aggregate(weight, by=list(dose), FUN=mean)fit <- aov(weight gesttime + dose) summary(fit)（1）評(píng)估檢驗(yàn)的假設(shè)條件ANCOVA還假定回歸斜率相同。本例中，假定四個(gè)處理組通過(guò)懷孕時(shí)間來(lái)預(yù)測(cè)出生體重的回歸斜率都相同。library(multcomp) fit2<-aov(weightgesttime*dose,data=litter) summary(fit2)（2）結(jié)果可視化HH包中的ancova()函數(shù)可以繪制

11、因變量、協(xié)變量和因子之間的關(guān)系圖。library(HH)ancova(weight gesttime + dose, data=litter)四、雙因素方差分析在 雙 因 素 方 差 分 析 中 ， 受 試 者 被 分 配 到 兩 因 子 的 交 叉 類(lèi) 別 組 中 。 以 基 礎(chǔ) 安 裝 中 的ToothGrowth數(shù)據(jù)集為例，隨機(jī)分配60只豚鼠，分別采用兩種喂食方法（橙汁或維生素C），各喂食方法中抗壞血酸含量有三種水平（0.5mg、 1mg或2mg），每種處理方式組合都被分配10只豚鼠。牙齒長(zhǎng)度為因變量attach(ToothGrowth)table(supp,dose)aggregate

12、(len, by=list(supp,dose), FUN=mean)aggregate(len, by=list(supp,dose), FUN=sd)dose <- factor(dose)fit <- aov(len supp*dose)summary(fit)五、重復(fù)測(cè)量方差分析所謂重復(fù)測(cè)量方差分析，即受試者被測(cè)量不止一次。本節(jié)重點(diǎn)關(guān)注含一個(gè)組內(nèi)和一個(gè)組間因子的重復(fù)測(cè)量方差分析（這是一個(gè)常見(jiàn)的設(shè)計(jì)）。示例來(lái)源于生理生態(tài)學(xué)領(lǐng)域，研究方向是生命系統(tǒng)的生理和生化過(guò)程如何響應(yīng)環(huán)境因素的變異（此為應(yīng)對(duì)全球變暖的一個(gè)非常重要的研究領(lǐng)域）?；A(chǔ)安裝包中的CO2數(shù)據(jù)集包含了北方和南方牧草類(lèi)

13、植物Echinochloa crus-galli （Potvin、Lechowicz、 Tardif， 1990）的寒冷容忍度研究結(jié)果，在某濃度二氧化碳的環(huán)境中，對(duì)寒帶植物與非寒帶植物的光合作用率進(jìn)行了比較。研究所用植物一半來(lái)自于加拿大的魁北克?。≦uebec），另一半來(lái)自美國(guó)的密西西比州（Mississippi）。首先，我們關(guān)注寒帶植物。因變量是二氧化碳吸收量（uptake），單位為ml/L，自變量是植物類(lèi)型Type（魁北克VS.密西西比）和七種水平（951000 umol/m2 sec）的二氧化碳濃度（conc）。另外， Type是組間因子， conc是組內(nèi)因子。 Type已經(jīng)被存儲(chǔ)為一個(gè)

14、因子變量，但你還需要先將conc轉(zhuǎn)換為因子變量含一個(gè)組間因子和一個(gè)組內(nèi)因子的重復(fù)測(cè)量方差分析CO2$conc <- factor(CO2$conc)w1b1 <- subset(CO2, Treatment='chilled')fit <- aov(uptake (conc*Type) + Error(Plant/(conc), w1b1)summary(fit)par(las=2)par(mar=c(10,4,4,2)with(w1b1, interaction.plot(conc,Type,uptake, type="b", col=c

15、("red","blue"), pch=c(16,18), main="Interaction Plot for Plant Type and Concentration")boxplot(uptake Type*conc, data=w1b1, col=(c("gold","green"), main="Chilled Quebec and Mississippi Plants", ylab="Carbon dioxide uptake rate (umol/m2

16、 sec)")par(opar)六、多元方差分析當(dāng)因變量（結(jié)果變量）不止一個(gè)時(shí)，可用多元方差分析（MANOVA）對(duì)它們同時(shí)進(jìn)行分析。以MASS包中的UScereal數(shù)據(jù)集為例（Venables， Ripley（1999），我們將研究美國(guó)谷物中的卡路里、脂肪和糖含量是否會(huì)因?yàn)閮?chǔ)存架位置的不同而發(fā)生變化；其中1代表底層貨架， 2代表中層貨架， 3代表頂層貨架?？防?、脂肪和糖含量是因變量，貨架是三水平（1、 2、 3）的自變量。library(MASS)attach(UScereal)shelf <- factor(shelf)y <- cbind(calories, fat, sugars)aggregate(y, by=list(shelf), FUN=mean)cov(y)fit <- manova(y shelf)summary(fit)summary.aov(fit)（1）評(píng)估假設(shè)檢驗(yàn)單因素多元方差分析有兩個(gè)前提假設(shè)，一個(gè)是多元正態(tài)性，一個(gè)是方差-協(xié)方差矩陣同質(zhì)性。第一個(gè)假設(shè)即指因變量組合成的向量服從一個(gè)多元正態(tài)分布?？梢杂肣-Q圖來(lái)檢驗(yàn)該假設(shè)條件center <- colMeans(y)n <- nrow(y)p <