在第六節(jié)t檢驗(yàn)部分的多組間差異分析提到涉及的內(nèi)容比較多蒋纬，并未記錄，這次專門系統(tǒng)學(xué)習(xí)一下相關(guān)知識(shí)坚弱。Analysis of Variance颠锉，簡(jiǎn)稱ANOVA

0、相關(guān)術(shù)語

• 方差分析（ANOVA）：多組別差異的分析史汗；
• 單因素組間方差分析（one-way ANOVA）：實(shí)驗(yàn)設(shè)計(jì)為一個(gè)因素，兩/多水平的組間因子拒垃。例如10位病人分為兩組停撞，每組5人，各接受一種治療方案悼瓮，評(píng)價(jià)組間區(qū)別戈毒。又稱為單因素組間方差分析；

特點(diǎn)：一個(gè)病人只測(cè)量一次横堡；實(shí)驗(yàn)分為幾組就是幾水平的埋市。（較常見我認(rèn)為）

• 單因素組內(nèi)方差分析：實(shí)驗(yàn)設(shè)計(jì)為兩或多水平的組內(nèi)因子。例如評(píng)價(jià)一種治療方案隨時(shí)間的變化情況命贴，10位病人為一組在所有時(shí)間水平下都進(jìn)行測(cè)量道宅。又因?yàn)槊總€(gè)病人不止被測(cè)了一次，也稱為重復(fù)測(cè)量方差分析胸蛛。

特點(diǎn)：一大組病人一起測(cè)污茵；測(cè)量幾次就是幾水平的。

• 雙因素方差分析：上述的治療方案與時(shí)間均考慮時(shí)葬项，即分析不同治療方案的影響和治療時(shí)間的影響泞当，又分析治療方案與時(shí)間的交互作用。前兩個(gè)稱為主效應(yīng)民珍，最后一個(gè)稱為交互效應(yīng)襟士。將會(huì)做分別3次F檢驗(yàn)：若療法結(jié)果顯著盗飒，則說明兩種方案治療效果不同；若時(shí)間結(jié)果顯著陋桂，則說明時(shí)間長(zhǎng)短會(huì)影響療效逆趣；若交互效應(yīng)影響顯著，則說明兩種方案隨著時(shí)間變化治療效果變化不同章喉；
• 雙因素混合模型方差分析：兩因子一個(gè)是組內(nèi)因子汗贫，一個(gè)是組間因子。
• 混淆因素：實(shí)驗(yàn)設(shè)計(jì)對(duì)象本身其它因素造成的對(duì)象起始水平就存在差異秸脱，從而也能一定程度解釋因變量的組間差異的因素落包。因?yàn)閷?duì)其不感興趣，也可以稱為干擾常數(shù)摊唇，可以進(jìn)行調(diào)整咐蝇；
• 協(xié)方差分析（ANCOVA）：約定好混淆因素為協(xié)變量，在對(duì)目標(biāo)因素進(jìn)行結(jié)果評(píng)價(jià)前巷查，對(duì)任何混淆因素水平的組間差異進(jìn)行的統(tǒng)計(jì)性調(diào)整的設(shè)計(jì)有序；
• 多元方差分析（MANOVA）：上述的例子實(shí)驗(yàn)結(jié)果只記錄單個(gè)因變量y。當(dāng)因變量不止一個(gè)（結(jié)果的多元分析）時(shí)岛请，即為多元方差分析旭寿。若協(xié)變量也存在，即稱為多元協(xié)方差分析崇败。
• aov方差分析函數(shù)表達(dá)式中效應(yīng)的順序要嚴(yán)格注意：越基礎(chǔ)性的效應(yīng)越需要放在前面盅称；一般首先是協(xié)變量--主效應(yīng)--雙因素交互項(xiàng)--三因素交互項(xiàng)

1、單因素方差分析

• 一般是指單因素組間方差分析
• 實(shí)驗(yàn)示例來自multcomp包的cholesterol的數(shù)據(jù)后室，記載了50位病人分為5組各接受一種trt后的膽固醇含量缩膝。

library(multcomp)  #第一次使用需要安裝
cholesterol

cholesterol

1.1、數(shù)據(jù)合理性檢驗(yàn)

單因素方差分析中岸霹，因變量要服從正態(tài)分布疾层；并且各組方差相等。

（1）Q-Q圖檢驗(yàn)正態(tài)分布

library(car)
qqPlot(lm(response ~ trt, data=cholesterol),    #注意qqPlot() 函數(shù)要求用lm()擬合
       simulate=TRUE, main="Q-Q Plot", labels=FALSE)
#圖形結(jié)果表示數(shù)據(jù)均落在95%的置信區(qū)間范圍內(nèi)贡避，符合正態(tài)性假設(shè)

Q-Q Plot.png

（2）bartlett方差齊性檢驗(yàn)

bartlett.test(response ~ trt, data=cholesterol)
#p值越大痛黎，表明5組的方差沒有顯著不同

bartlett.test.png

• 值得注意的是方差齊性分析對(duì)離群點(diǎn)分析非常敏感，利用car包的outlierTest()函數(shù)檢測(cè)離群點(diǎn)

outlierTest(fit)
#結(jié)果表明 No Studentized residuals with Bonferroni

1.2刮吧、aov()單因素方差分析

attach(cholesterol)  
table(trt)     #查看水平分組以及各組樣本大芯艘荨（一般設(shè)計(jì)相同，為均衡設(shè)計(jì)）
aggregate(response, by=list(trt), FUN=mean)  #計(jì)算各組/水平平均值
aggregate(response, by=list(trt), FUN=sd)    #計(jì)算各組/水平方差
fit <- aov(response ~ trt)     #進(jìn)行單因素方差分析 皇筛，注意格式
summary(fit)

summary(fit)

結(jié)果表明ANOVA對(duì)trt（治療方式）的F檢驗(yàn)非常顯著（p<0.0001）,說明五種療法的效果不同

#將上述結(jié)果可視化
library(gplots)    #繪制帶有置信區(qū)間的組均值圖形
png("1.png")   #剛才在線顯示有問題琉历，存在本地顯示正常
plotmeans(response ~ trt, xlab="Treatment", ylab="Response", 
          main="Mean Plot\nwith 95% CI")
#繪制組均值+置信區(qū)間圖形
detach(cholesterol)
dev.off()

1.png

1.3、深入分析多重比較

• 在上述分析的基礎(chǔ)上，進(jìn)一步分析兩組的差異顯著與否旗笔，即具體哪種療法與其它療法不同彪置。有以下兩種方法--

#法1
TukeyHSD(fit)  # p adj 值越小，越顯著蝇恶；即該兩組均值差異顯著
#將上述數(shù)據(jù)可視化
opar <- par(no.readonly = TRUE)
par(las=2)      # 旋轉(zhuǎn)軸標(biāo)簽（主要是便于縱坐標(biāo)便簽展示）
par(mar=c(5,8,4,2))   #增大左邊界的面積（也是針對(duì)縱坐標(biāo)）
plot(TukeyHSD(fit)) 
#凡是置信區(qū)間包含0的（靠左邊的）說明差異不顯著（p>0.5）
par(opar)

TukeyHSD.png

#法2
library(multcomp)
opar <- par(no.readonly = TRUE)
par(mar=c(5,4,6,2))   #擴(kuò)大頂部邊界面積
tuk <- glht(fit, linfct=mcp(trt="Tukey"))   
plot(cld(tuk, level=.05),col="lightgrey") 
#cld()函數(shù)中的level選項(xiàng)設(shè)置了使用的顯著水平（0.05拳魁，即本例中的95%的置信區(qū)間）
par(opar)   

如下圖，圖形上方中撮弧，有相同字母的組說明均值差異不顯著潘懊。該圖此外也能反映各組的分布情況，因此信息更多贿衍。

箱圖.png

2授舟、單因素協(xié)方差分析（ANCOVA）

• litter為一實(shí)驗(yàn)數(shù)據(jù)框，列名為劑量贸辈、體重释树，懷孕時(shí)間（最后一列應(yīng)該是產(chǎn)崽數(shù)量，不關(guān)注）擎淤；目的為研究懷孕小鼠接受不同劑量藥物產(chǎn)崽體重的變化情況曹铃。懷孕時(shí)間為協(xié)變量街望。

  
  
  litter

2.1咙崎、數(shù)據(jù)合理性檢驗(yàn)

• ANCOVA也要進(jìn)行正態(tài)分布與齊方差檢驗(yàn)胳喷，步驟同前，不再贅述席吴；

同質(zhì)性檢驗(yàn)

• 此外還要進(jìn)行協(xié)變量的回歸率同質(zhì)性檢驗(yàn)正驻，協(xié)變量影響與劑量主因素不存在顯著的交互作用，即懷孕時(shí)間對(duì)體重的影響不會(huì)因劑量分組的不同而不同抢腐。

library(multcomp)
data(litter, package="multcomp") 
fit <- aov(weight ~ gesttime*dose, data=litter)
summary(fit)

注意，上述表達(dá)式 y ~ A*B 為雙因素方差分析表達(dá)式襟交。我認(rèn)為這里只是為了單純看下交互項(xiàng)是否顯著迈倍，而使用這個(gè)表達(dá)式的。

summary(fit)

• 返回結(jié)果表明交互效應(yīng)不顯著捣域，支持了斜率相等的假設(shè)啼染。此外也可繪制因變量、協(xié)變量與因子之間的關(guān)系圖焕梅。

library(HH)
ancova(weight ~ gesttime + dose, data=litter)
#可以清楚看出回歸線相互平行“同質(zhì)性”迹鹅，而截距項(xiàng)不同（劑量因素不同導(dǎo)致）

同質(zhì)性檢驗(yàn).png

2.2、ANCOVA分析

• 這里強(qiáng)調(diào)下表達(dá)式順序贞言，之前為單因素分析僅一個(gè)因子斜棚。當(dāng)存在協(xié)變量時(shí)，協(xié)變量要放在前面，其次是主效應(yīng)弟蚀，再是雙因素交互項(xiàng)

attach(litter)
table(dose)  #按實(shí)驗(yàn)劑量分為四組蚤霞，發(fā)現(xiàn)每組的數(shù)量都不同（為非均衡設(shè)計(jì)）
aggregate(weight, by=list(dose), FUN=mean)
fit2 <- aov(weight ~ gesttime + dose)   # 協(xié)變量+單因素  注意順序！
summary(fit2)
# F檢驗(yàn)結(jié)果表明懷孕時(shí)間域產(chǎn)崽體重相關(guān)义钉；控制懷孕時(shí)間（協(xié)變量）昧绣，藥物劑量與出生體重有關(guān)。

平均值1

summary(fit2)

2.3捶闸、衍生計(jì)算

（1）獲得去除協(xié)變量效應(yīng)的組均值

我的理解是將小鼠的懷孕時(shí)間都變換為相同值時(shí)后的劑量影響值

library(effects)
effect("dose", fit2)
#可以對(duì)比一下之間用aggregate()計(jì)算得到的平均值

平均值2

（2）化繁為簡(jiǎn)實(shí)驗(yàn)?zāi)康?/h3>

• 研究用藥與否對(duì)產(chǎn)崽體重的影響夜畴，即比較第一組（劑量為0）與第二~四組（劑量為5、50,500）的差異

library(multcomp)
contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))  
#設(shè)定為第一組與其它三組的均值比較
summary(glht(fit2, linfct=mcp(dose=contrast)))
#假設(shè)檢驗(yàn)的t統(tǒng)計(jì)量（2.581）在p<0.05水平下顯著删壮，可以得出未用藥比用藥組體重高的結(jié)論贪绘。

no drug vs. drug

3、雙因素方差分析

• 示例實(shí)驗(yàn)來自基礎(chǔ)數(shù)據(jù)ToothGrowth醉锅，為研究?jī)煞N試劑（supp）的各自三種劑量水平（dose）對(duì)豚鼠牙齒長(zhǎng)度的影響（每個(gè)豚鼠只測(cè)一次兔簇，為組間因子）
  
  ToothGrowth
  
  注意：雙因素方差分析時(shí)，要保證兩組樣本量相同硬耍；否則y ~ A:B與y ~ B:A兩個(gè)模型結(jié)果會(huì)有差異垄琐。

attach(ToothGrowth)
aggregate(len, by=list(supp,dose), FUN=mean)    #共3*2=6組
aggregate(len, by=list(supp,dose), FUN=sd)
dose <- factor(dose)  
#ToothGrowth數(shù)據(jù)框中supp已經(jīng)默認(rèn)為因子了，這樣aov() 函數(shù)就會(huì)將它當(dāng)做一個(gè)分組變量经柴，而不是一個(gè)數(shù)值型協(xié)變量.
fit <- aov(len ~ supp*dose)  #兩因素的函數(shù)表達(dá)格式
summary(fit)  
#結(jié)果表明主效應(yīng)與交互效應(yīng)都非常顯著

summary(fit)

• 可視化一下~

#法1
interaction.plot(dose, supp, len, type="b", 
                 col=c("red","blue"), pch=c(16, 18),
                 main = "Interaction between Dose and Supplement Type")
#法2
library(gplots)
plotmeans(len ~ interaction(supp, dose, sep=" "),
          connect=list(c(1, 3, 5),c(2, 4, 6)), 
          col=c("red","darkgreen"),
          main = "Interaction Plot with 95% CIs", 
          xlab="Treatment and Dose Combination")
#觀察與上圖的區(qū)別（顯示要問題的話配合 png("2.png")  dev.off()下載到本地）
#法3
library(HH)
interaction2wt(len~supp*dose)
#顯示為四幅圖的匯總大圖狸窘；箱線圖為控制一種水平下，另一種因素不同水平的評(píng)價(jià)坯认；折線圖為交互效應(yīng)互相影響變化圖翻擒。

1.png

2.png

3.png

4、重復(fù)測(cè)量方差分析（組內(nèi)因子）

• 常見的實(shí)驗(yàn)設(shè)計(jì)為含一個(gè)組內(nèi)因子牛哺，一個(gè)組間因子重復(fù)測(cè)量方差分析陋气。比如本節(jié)一開始的論述兩種治療方法隨時(shí)間長(zhǎng)短變化的影響。
• 這里示例實(shí)驗(yàn)來自基礎(chǔ)函數(shù)CO2引润，此處我們?nèi)パ芯績(jī)煞N植物（Type巩趁，為組間因子）在不同的CO2濃度范圍（conc，為組內(nèi)因子淳附，重復(fù)測(cè)量）下二氧化碳吸收量（uptake）的變化议慰。

CO2$conc <- factor(CO2$conc)
w1b1 <- subset(CO2, Treatment=='chilled')  #抽取目標(biāo)行數(shù)據(jù)
fit <- aov(uptake ~ (conc*Type) + Error(Plant/(conc)), w1b1) 
#注意含單個(gè)組間因子（Type）、單個(gè)組內(nèi)因子（conc）的重復(fù)測(cè)量方差分析的表達(dá)格式奴曙。
#(conc*Type)里兩者的順序可以顛倒别凹；
#Error(Plant/(conc))中Plant是相對(duì)于每組conc中的獨(dú)有的標(biāo)識(shí)變量。即重復(fù)測(cè)量的單位
summary(fit)
# 結(jié)果表明在0.01的水平下洽糟，主效應(yīng)type與conc以及交互效應(yīng)都非常顯著

summary(fit)

• 可視化一下~

par(las=2)
par(mar=c(10,4,4,2))
#法1：點(diǎn)線圖
with(w1b1, 
     interaction.plot(conc,Type,uptake, 
                      type="b", col=c("red","blue"), pch=c(16,18),
                      main="Interaction Plot for Plant Type and Concentration"))
#法2：箱線圖
boxplot(uptake ~ Type*conc, data=w1b1, col=(c("gold","green")),
        main="Chilled Quebec and Mississippi Plants", 
        ylab="Carbon dioxide uptake rate (umol/m^2 sec)")
par(opar)

1.png

2.png

結(jié)合上圖及分析結(jié)果可得結(jié)論：兩種植物吸收二氧化碳能力差異顯著炉菲，且隨二氧化碳濃度升高堕战，差異越明顯。

5颁督、多元方差分析（MANOVA）

• 因變量（結(jié)果變量）不止一個(gè)時(shí)（y1践啄，y2...），可用多元方差對(duì)它們同時(shí)進(jìn)行分析
• 示例實(shí)驗(yàn)設(shè)計(jì)為不同儲(chǔ)架層谷物（自變量）里沉御，卡路里calories屿讽、脂肪fat以及糖sugars含量（3個(gè)因變量）是否存在差異

5.1、多元正態(tài)性檢驗(yàn)

library(MASS)
attach(UScereal)    #數(shù)據(jù)包
shelf <- factor(shelf) #
y <- cbind(calories, fat, sugars)  #將3個(gè)因變量合并為一個(gè)目標(biāo)矩陣
center <- colMeans(y)
n <- nrow(y) #行數(shù)
p <- ncol(y) #列數(shù)
cov <- cov(y) 
d <- mahalanobis(y,center,cov)
coord <- qqplot(qchisq(ppoints(n),df=p),
                d, main="QQ Plot Assessing Multivariate Normality",
                ylab="Mahalanobis D2")
abline(a=0,b=1)  
#繪制y=x 的參考線吠裆；若數(shù)據(jù)服從多元正態(tài)分布伐谈，則點(diǎn)將落在直線上。
#比如本圖中Wheaties Honey Gold點(diǎn)與Wheaties點(diǎn)異常试疙，違反多元正態(tài)性诵棵。可以刪除這兩個(gè)點(diǎn)再重新分析祝旷。
identify(coord$x, coord$y, labels=row.names(UScereal))  #交互式標(biāo)注重點(diǎn)關(guān)注的點(diǎn)

Q-Q Plot

5.2履澳、MANOVA分析

aggregate(y, by=list(shelf), FUN=mean)  #獲取各個(gè)貨架的各個(gè)均值
cov(y)  #輸出個(gè)谷物間的方差和協(xié)方差（是什么？存疑）
fit <- manova(y ~ shelf)  #多元方差分析
summary(fit) 
#返回結(jié)果說明三個(gè)貨架的營(yíng)養(yǎng)成分測(cè)量值不同
summary.aov(fit)
#對(duì)每一個(gè)因變量做單因素方差分析怀跛，返回結(jié)果表明三個(gè)貨架中每種營(yíng)養(yǎng)成分的測(cè)量值都是不同的

summary(fit)

summary.aov(fit)

5.3距贷、穩(wěn)健單因素MANOVA

• 在存在不滿足多元正態(tài)性、方差-協(xié)方差均值假設(shè)吻谋，以及多元離群點(diǎn)忠蝗，可以考慮使用本方法--穩(wěn)健或非參數(shù)版本的MANOVA檢驗(yàn)

library(rrcov)
Wilks.test(y,shelf, method="mcd")   #計(jì)算過程需要等幾十秒時(shí)間（也是第一次遇到）~ 

穩(wěn)健MANOVA

• 返回結(jié)果表明穩(wěn)健檢驗(yàn)對(duì)離群點(diǎn)和違反MANOVA假設(shè)的情況不敏感，也驗(yàn)證了存儲(chǔ)在不同層貨架的營(yíng)養(yǎng)含量不同漓拾。

教材中還提到一種用回歸方法做ANOVA的介紹（P220）
參考書籍《R語言實(shí)戰(zhàn)（第2版）》