R語言方差分析ANOVA

自己整理編寫的R語言常用數(shù)據(jù)分析模型的模板，原文件為Rmd格式，直接復(fù)制粘貼過來，作為個(gè)人學(xué)習(xí)筆記保存和分享。
I. 單因素方差分析
#用data frame的格式輸入數(shù)據(jù)
medicine <- data.frame(
   Response=c(7,5,3,1,6,5,3,3,7,9,9,9,4,3,4,3),
   Treatment=factor(c(rep(1,4),rep(2,4),rep(3,4),rep(4,4)))
   )       

#各組樣本大小
table(medicine$Treatment)
#各組的均值
aggregate(medicine$Response,by=list(medicine$Treatment),FUN=mean)
#各組的標(biāo)準(zhǔn)差
aggregate(medicine$Response,by=list(medicine$Treatment),FUN=sd)
#調(diào)用aov函數(shù)進(jìn)行方差分析（檢驗(yàn)組間差異）
medicine.aov <- aov(Response ~ Treatment,data=medicine)
#summary提取方差分析的結(jié)果
summary(medicine.aov)

分析上述計(jì)算結(jié)果，Df表示自由度，Sum Sq 表示平方和，Mean Sq 表示均方，F(xiàn) value 是F值，Pr(>F)是p值，A即為因子A，Residuals 是殘差。

但是我們注意到，這個(gè)結(jié)果并不完整。直接用summary()函數(shù)時(shí)候，只有因素A和誤差兩行，沒有總和，這里編個(gè)小程序(anova.tab.R)作改進(jìn)，計(jì)算方法為：將summary函數(shù)得到表中的第一行與第二行求和，得到總和行的值。

#anova.tab.R程序
anova.tab <- function(fm){
  tab <- summary(fm)
  k <- length(tab[[1]]-2)
  temp <- c(sum(tab[[1]][,1]),sum(tab[[1]][,2]),rep(NA,k))
  tab[[1]]["Total",] <- temp
}

將anova.tab.R函數(shù)保存在工作目錄中。

getwd()
#利用anova.tab.R函數(shù)，得到完整的方差分析表
source("anova.tab.R");anova.tab(medicine.aov)
#畫圖
plot(medicine$Response~medicine$Treatment)
#繪制各組均值及其置信區(qū)間的圖形
library(gplots)
plotmeans(medicine$Response~medicine$Treatment,xlab = "Treatment",ylab = "Response",main = "Mean Plot\nwith 95% CI")

1.多重比較

ANOVA對各療法的F檢驗(yàn)表明，4種藥品用于緩解術(shù)后疼痛的療效不同，但是并不能得出哪種藥品療法與其他不同。多重比較可以解決這個(gè)問題.e.g. TukeyHSD()函數(shù)提供了對各組均值差異的成對檢驗(yàn)；multcomp包中的glht()函數(shù)提供了多重均值比較更為全面的方法，既適用于線性模型，也適用于廣義線性模型；多重t檢驗(yàn)方法針對每組數(shù)據(jù)進(jìn)行t檢驗(yàn)。代碼如下：

TukeyHSD(medicine.aov)
#par()函數(shù)旋轉(zhuǎn)軸標(biāo)簽，增大左邊界面積，使標(biāo)簽擺放更美觀。
par(las = 2)
par(mar = c(5, 8, 4, 2))
plot(TukeyHSD(medicine.aov))

圖形中置信區(qū)間包含0的藥品對比，說明差異不顯著。
library(multcomp)
#為適合字母陣列擺放，par語句用來增大頂部邊界面積
par(mar = c(5, 4, 6, 2))
tuk <- glht(medicine.aov, linfct = mcp(Treatment = "Tukey"))
#cld()函數(shù)中l(wèi)evel選項(xiàng)為設(shè)置的顯著性水平（這里的0.05對應(yīng)95%置信區(qū)間）
plot(cld(tuk, level = 0.05), col = "lightgrey")

有相同字母的組（用箱線圖表示）說明均值差異不顯著。

多次重復(fù)使用t檢驗(yàn)會(huì)增大犯第一類錯(cuò)誤的概率，為了克服這一缺點(diǎn)，需要調(diào)整p-值。R軟件調(diào)整p-值用的是p.adjust()函數(shù)，函數(shù)使用的不同參數(shù)代表不同的調(diào)整方法。
attach(medicine)
#求數(shù)據(jù)在各水平下的均值
mu<-c(mean(Response[Treatment==1]), mean(Response[Treatment==2]), mean(Response[Treatment==3]),mean(Response[Treatment==4])); mu
#作多重t檢驗(yàn)。這里用到的pairwise.t.test()函數(shù)用來得到多重比較的p值
pairwise.t.test(Response, Treatment, p.adjust.method = "none")

#觀察兩個(gè)作調(diào)整后的p值的情況。p.adjust.method()函數(shù)的參數(shù)也可換為"hochberg","hommel","bonferroni","BH","BY","fdr"等。
pairwise.t.test(Response, Treatment, p.adjust.method = "holm")
#繪制箱線圖
plot(medicine$Response~medicine$Treatment)

從上述結(jié)果可見，124無顯著差異，3與124均有顯著差異，即緩解疼痛的4種藥品，3與124有顯著差異，124間差異不顯著

2.評估檢驗(yàn)的假設(shè)條件

擬合結(jié)果的可信度來源于，做統(tǒng)計(jì)檢驗(yàn)時(shí)數(shù)據(jù)滿足假設(shè)條件的程度

(1)誤差的正態(tài)性檢驗(yàn)

單因素方差分析中，我們假設(shè)因變量服從正態(tài)分布，各組方差相等。可用Q-Q圖來檢驗(yàn)正態(tài)性假設(shè)。擬合診斷如下：
library(car)
qqPlot(lm(Response ~ Treatment, data = medicine), simulate = TRUE,
    main = "QQ Plot", labels = FALSE)

數(shù)據(jù)幾乎都落在95%的置信區(qū)間范圍內(nèi)，說明滿足正態(tài)性假設(shè)

也可用W檢驗(yàn)（shapiro.test()函數(shù)）方法對數(shù)據(jù)作正態(tài)性檢驗(yàn)

attach(medicine)
shapiro.test(Response[Treatment==1])
shapiro.test(Response[Treatment==2])
shapiro.test(Response[Treatment==3])
shapiro.test(Response[Treatment==4])

計(jì)算結(jié)果表明，數(shù)據(jù)在四種水平下的均是正態(tài)的

(2)方差的其次性檢驗(yàn)

方差的其次性檢驗(yàn)就是檢驗(yàn)數(shù)據(jù)在不同的水平下方差是否相同，常用方法是Bartlett檢驗(yàn)

#R里用bartlett.test()函數(shù)來提供Bartlett檢驗(yàn)。另外還有Fligner-Killeen檢驗(yàn)（fligner.test()函數(shù))和Brown-Forsythe檢驗(yàn)（HH包中的hov()函數(shù)）
bartlett.test(Response~Treatment,data=medicine)

p值0.1285>0.05，接受原假設(shè)，認(rèn)為各組的數(shù)據(jù)是等方差的

方差其次性分析對離群點(diǎn)非常敏感，可用car包的outlierTest()函數(shù)來檢測離群點(diǎn)

library(car)
outlierTest(medicine.aov)

從p值結(jié)果看，并沒有證據(jù)說明該數(shù)據(jù)中含有離群點(diǎn)

根據(jù)Q-Q圖，Bartlett檢驗(yàn)和離群點(diǎn)檢驗(yàn)，該數(shù)據(jù)似乎可以用ANOVA模型擬合得很好，這些方法反過來增強(qiáng)了我們對于所得結(jié)果的信心

數(shù)據(jù)的總體分布類型未知；或數(shù)據(jù)的總體分布類型已知，但不符合正態(tài)分布；或某些變量可能無法精確測量時(shí)，可以使用非參數(shù)統(tǒng)計(jì)方法.非參數(shù)統(tǒng)計(jì)是拋開總體分布類型不考慮，對總體參數(shù)不做比較，比較的是總體分布的位置是否相同的統(tǒng)計(jì)方法．秩和檢驗(yàn)是非參數(shù)統(tǒng)計(jì)中一種經(jīng)常使用的檢驗(yàn)方法．這里的“秩”又可被稱為等級，即按照數(shù)據(jù)大小排定的次序號(hào)．此次序號(hào)的總和被稱為“秩和”．

方差分析過程需要滿足若干條件，F(xiàn)檢驗(yàn)才能奏效，可惜有時(shí)候采集到的數(shù)據(jù)并不能滿足這樣的要求。像兩樣本比較時(shí)一樣，嘗試將數(shù)據(jù)轉(zhuǎn)換為秩統(tǒng)計(jì)量，因?yàn)橹冉y(tǒng)計(jì)量的分布與總體分布無關(guān)，這樣就可以避開總體分布的要求．上述問題就可以通過數(shù)據(jù)的秩統(tǒng)計(jì)量就解決了。在比較兩個(gè)以上的總體時(shí)，廣泛使用的是Kruskal-Wallis秩和檢驗(yàn)，它是對兩個(gè)以上樣本進(jìn)行比較的非參數(shù)檢驗(yàn)方法。實(shí)質(zhì)上，它是兩樣本的Wilcoxon方法在多于兩個(gè)樣本時(shí)的推廣。

R軟件提供了Kruskal-Wallis秩和檢驗(yàn)，函數(shù)為kruskal.test()
(3)Kruskal-Wallis秩和檢驗(yàn)

medicine <- data.frame(
   Response=c(7,5,3,1,6,5,3,3,7,9,9,9,4,3,4,3),
   Treatment=factor(c(rep(1,4),rep(2,4),rep(3,4),rep(4,4)))
   )  
kruskal.test(Response~Treatment,data=medicine)
p值=0.0344<0.05，拒絕原假設(shè)，認(rèn)為四種藥物緩解疼痛效果有顯著差異
該函數(shù)還有另外兩種寫法如下：
kruskal.test(medicine$Response,medicine$Treatment)
A <- c(7,5,3,1)
B <- c(6,5,3,3)
C <- c(7,9,9,9)
D <- c(4,3,4,3)
kruskal.test(list(A,B,C,D))

之后再對上述數(shù)據(jù)作正太檢驗(yàn)和方差齊次檢驗(yàn)，如果全部通過檢驗(yàn)，則該數(shù)據(jù)也可以作方差分析

(4)Friedman秩和檢驗(yàn)

在配伍組設(shè)計(jì)中，多個(gè)樣本的比較，如果它們的總體不能滿足正態(tài)性和方差齊性的要求，可采用Friedman秩和檢驗(yàn)

Friedman秩和檢驗(yàn)的基本思想與前面介紹的方法類似，但是配伍組設(shè)計(jì)的隨機(jī)化是在配伍組內(nèi)進(jìn)行的，而配伍組間沒有進(jìn)行隨機(jī)化。因此在進(jìn)行Friedman秩和檢驗(yàn)時(shí)，是分別在每個(gè)配伍組里將數(shù)據(jù)從小到大編秩，如果相同的數(shù)據(jù)取平均秩次。

r軟件中，函數(shù)friedman.test()提供了Friedman秩和檢驗(yàn)

medicine.matrix <- matrix(
   c(7,5,3,1,6,5,3,3,7,9,9,9,4,3,4,3),
   ncol = 4,dimnames = list(1:4,c("A","B","C","D"))
   )  
friedman.test(medicine.matrix)

該函數(shù)還有另外兩種寫法如下：

x <- c(7,5,3,1,6,5,3,3,7,9,9,9,4,3,4,3)
#4行4列，每行4個(gè)數(shù)據(jù)，總共16個(gè)
g <- gl(4,4)
b <- gl(4,1,16)
friedman.test(x,g,b)
medicine <- data.frame(
   x=c(7,5,3,1,6,5,3,3,7,9,9,9,4,3,4,3),
   g = gl(4,4),b = gl(4,1,16)
   )
friedman.test(x~g|b,data = medicine)

3.單因素協(xié)方差分析（顯著因素下的水平間差異檢驗(yàn)）

單因素協(xié)方差分析（ANCOVA）擴(kuò)展了單因素方差分析（ANOVA），包含一個(gè)或多個(gè)定量的協(xié)變量。下面的例子來自于multcomp包中的litter數(shù)據(jù)集，懷孕小鼠被分為四個(gè)小組，每個(gè)小組接受不同劑量（0、5、50、500）的藥物處理，產(chǎn)下幼崽的體重均值為因變量，懷孕時(shí)間為協(xié)變量。

（1）單因素ANCOVA

data(litter, package = "multcomp")
attach(litter)
#table()函數(shù)，看到每種劑量下所產(chǎn)幼崽數(shù)量并不同
table(dose)
#aggrgate()函數(shù)獲得各組均值，可以發(fā)現(xiàn)未用藥組幼崽體重均值最高
aggregate(weight, by = list(dose), FUN = mean)
fit <- aov(weight ~ gesttime + dose)
summary(fit)

由于使用了協(xié)變量，如果想要獲取調(diào)整的組均值–即去除協(xié)變量效應(yīng)后的組均值，可使用effects包中的effects()函數(shù)來計(jì)算調(diào)整的均值

library(effects)
effect("dose",fit)

（2）對用戶定義的對照的多重比較

想得知具體哪種處理方式與其他不同，使用multcomp包來對所有均值進(jìn)行成對比較（多重比較）

library(multcomp)
contrast <- rbind(`no drug vs. drug` = c(3, -1, -1, -1))
summary(glht(fit, linfct = mcp(dose = contrast)))

對照c(3, -1, -1, -1)設(shè)定第一組與其他三組飛均值進(jìn)行比較。其他對照可用rbind()函數(shù)添加。從結(jié)果來看，假設(shè)檢驗(yàn)的t統(tǒng)計(jì)量在p<0.05水平下顯著，可以得出未用藥組比其他用藥條件下的出生體重高的結(jié)論
（3）評估檢驗(yàn)的假設(shè)條件–檢驗(yàn)同歸斜率的同質(zhì)性

ANCOVA與ANOVA相同，都需要正態(tài)性和方差齊次性假設(shè)，可用上述ANOVA的假設(shè)檢驗(yàn)的相同步驟來檢驗(yàn)。另外ANCOVA還假定回歸斜率相同。ANCOVA模型包含懷孕時(shí)間*劑量的交互項(xiàng)時(shí)，可對回歸斜率的同質(zhì)性進(jìn)行檢驗(yàn)。交互效應(yīng)若顯著，則意味著時(shí)間和幼崽出生體重間的關(guān)系依賴于藥物劑量的水平

library(multcomp)
fit2 <- aov(weight ~ gesttime * dose)
summary(fit2)

結(jié)果可以看到交互效應(yīng)不顯著，支持了斜率相等的假設(shè)。若假設(shè)不成立，可以嘗試變換協(xié)變量或因變量，或使用能對每個(gè)斜率獨(dú)立解釋的模型，或使用不需要假設(shè)回歸斜率同質(zhì)性的非參數(shù)ANCOVA方法。（如sm包中的sm.ancova()函數(shù)）

（4）結(jié)果可視化

HH包中的ancova()函數(shù)可以繪制因變量、協(xié)變量和因子之間的關(guān)系圖，代碼如下：

library(HH)
ancova(weight ~ gesttime + dose, data = litter)

從圖中可看出，用懷孕時(shí)間來預(yù)測出生體重的回歸線相互平行，只是截距項(xiàng)不同。隨著懷孕時(shí)間增加，幼崽出生體重也會(huì)增加。另外，還可以看到0劑量組截距項(xiàng)最大，5劑量組截距項(xiàng)最小。由于之前的設(shè)置，直線會(huì)保持平行，若用anvova(weight~gesttime*dose),生成的圖形將允許斜率和截距項(xiàng)依據(jù)組別而發(fā)生變化，這對可視化那些違背回歸斜率同質(zhì)性的實(shí)例非常有用

II.雙因素方差分析

1.不考慮交互作用

SAS自帶數(shù)據(jù)集sashelp.class中包含了學(xué)生的姓名、性別與身高。導(dǎo)出數(shù)據(jù)存為csv格式，現(xiàn)在分析年齡與性別是否是影響體重的顯著因素。該問題屬于不均衡數(shù)據(jù)集的方差分析

class <- read.csv("class.csv",header=T)
#預(yù)處理表明該設(shè)計(jì)不是均衡設(shè)計(jì)（各設(shè)計(jì)單元中樣本大小不一致）
table(class$Sex,class$Age)
#獲得各單元的均值和標(biāo)準(zhǔn)差
aggregate(class$Weight,by=list(class$Sex,class$Age),FUN=mean)
aggregate(class$Weight,by=list(class$Sex,class$Age),FUN=sd)
#作雙因素方差分析
class.aov <- aov(Weight ~ Sex+Age,data=class)
#調(diào)用自變函數(shù)anova.tab()，顯示計(jì)算結(jié)果
source("anova.tab.R");anova.tab(class.aov)

根據(jù)p值不同說明年齡和性別對體重有顯著影響

2.考慮交互作用

（1）3種方式對結(jié)果進(jìn)行可視化處理

用interaction.plot()函數(shù)來展示雙因素方差分析的交互效應(yīng)

interaction.plot(class$Sex,class$Age,class$Weight, type = "b", col = c("red",  "blue"), pch = c(16, 18), main = "Interaction between Dose and Supplement Type")

圖形展示了各年齡下，學(xué)生體重的均值

或者用gplots包中的plotmeans()函數(shù)來展示交互效應(yīng)

library(gplots)
plotmeans(class$Weight ~ interaction(class$Sex,class$Age, sep = ","),
    connect = list(c(1, 3, 5), c(2, 4, 6)),
    col = c("red", "darkgreen"),
    main = "Interaction Plot with 95% CIs",
    xlab = "Sex and Age Combination")

圖形展示了均值、誤差棒（95%CI）和樣本大小

用HH包中的interaction2wt()函數(shù)來可視化結(jié)果，圖形對任意順序的因子設(shè)計(jì)的主效應(yīng)和交互效應(yīng)都會(huì)進(jìn)行展示
library(HH)
interaction2wt(class$Weight ~ class$Sex*class$Age)

（2）有交互作用的方差分析

數(shù)據(jù)集fruit記錄了在不同濕度和溫度下某種植物的查處。這是一個(gè)雙因素方差分析的情形。假設(shè)方差分析的假設(shè)條件滿足，在顯著性水平0.05的前提下，欲分析不同溫度、不同濕度下產(chǎn)出是否有顯著差異，以及溫度和濕度的交互是否顯著差異，如果交互有差異，分析在濕度一定的情況下，溫度對產(chǎn)出的影響。

fruit <- read.csv("fruit.csv",header=T)
#output分別對于A、B、A&B的方差檢驗(yàn)
fruit.aov <- aov(output_lbs ~ humidity+temperature+humidity:temperature, data=fruit)
source("anova.tab.R"); anova.tab(fruit.aov)
output對于A&B高度顯著，說明交互效應(yīng)顯著

對于存在交互作用的兩因素，我們應(yīng)當(dāng)固定一個(gè)因素的水平，對另一個(gè)因素的水平進(jìn)行水平間差異檢驗(yàn)?
library(effects)
effect("humidity",fruit.aov)

SUMMARY:方差分析是一種常見的統(tǒng)計(jì)模型，用于檢驗(yàn)樣本間均值是否相等。方差分析適用于處理因素類型為分類變量、響應(yīng)變量類型為連續(xù)的情形。根據(jù)因素個(gè)數(shù)，方差分析可以分為單因素方差分析與多因素方差分析。在多因素方差分析中，要特別注意判斷因素間是否存在交互作用。此外，在實(shí)際應(yīng)用中，可以通過設(shè)計(jì)合理的試驗(yàn)，在盡可能排除外部因素的干擾后，再對試驗(yàn)數(shù)據(jù)進(jìn)行方差分析，這樣結(jié)果會(huì)更準(zhǔn)確。
write.csv(medcine,"test_medcine.csv")
write.csv(class,"test_class.csv")