R語言方差分析ANOVA

自己整理編寫的R語言常用數(shù)據(jù)分析模型的模板，原文件為Rmd格式，直接復制粘貼過來，作為個人學習筆記保存和分享。
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)
#各組的標準差
aggregate(medicine$Response,by=list(medicine$Treatment),FUN=sd)
#調(diào)用aov函數(shù)進行方差分析（檢驗組間差異）
medicine.aov <- aov(Response ~ Treatment,data=medicine)
#summary提取方差分析的結(jié)果
summary(medicine.aov)

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

但是我們注意到，這個結(jié)果并不完整。直接用summary()函數(shù)時候，只有因素A和誤差兩行，沒有總和，這里編個小程序(anova.tab.R)作改進，計算方法為：將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檢驗表明，4種藥品用于緩解術后疼痛的療效不同，但是并不能得出哪種藥品療法與其他不同。多重比較可以解決這個問題.e.g. TukeyHSD()函數(shù)提供了對各組均值差異的成對檢驗；multcomp包中的glht()函數(shù)提供了多重均值比較更為全面的方法，既適用于線性模型，也適用于廣義線性模型；多重t檢驗方法針對每組數(shù)據(jù)進行t檢驗。代碼如下：

TukeyHSD(medicine.aov)
#par()函數(shù)旋轉(zhuǎn)軸標簽，增大左邊界面積，使標簽擺放更美觀。
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選項為設置的顯著性水平（這里的0.05對應95%置信區(qū)間）
plot(cld(tuk, level = 0.05), col = "lightgrey")

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

多次重復使用t檢驗會增大犯第一類錯誤的概率，為了克服這一缺點，需要調(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檢驗。這里用到的pairwise.t.test()函數(shù)用來得到多重比較的p值
pairwise.t.test(Response, Treatment, p.adjust.method = "none")

#觀察兩個作調(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.評估檢驗的假設條件

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

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

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

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

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

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

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

(2)方差的其次性檢驗

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

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

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

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

library(car)
outlierTest(medicine.aov)

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

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

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

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

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

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ù)還有另外兩種寫法如下：
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ù)作正太檢驗和方差齊次檢驗，如果全部通過檢驗，則該數(shù)據(jù)也可以作方差分析

(4)Friedman秩和檢驗

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

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

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

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個數(shù)據(jù)，總共16個
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é)方差分析（顯著因素下的水平間差異檢驗）

單因素協(xié)方差分析（ANCOVA）擴展了單因素方差分析（ANOVA），包含一個或多個定量的協(xié)變量。下面的例子來自于multcomp包中的litter數(shù)據(jù)集，懷孕小鼠被分為四個小組，每個小組接受不同劑量（0、5、50、500）的藥物處理，產(chǎn)下幼崽的體重均值為因變量，懷孕時間為協(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é)變量效應后的組均值，可使用effects包中的effects()函數(shù)來計算調(diào)整的均值

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

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

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

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

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

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

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

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

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

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

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

II.雙因素方差分析

1.不考慮交互作用

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

class <- read.csv("class.csv",header=T)
#預處理表明該設計不是均衡設計（各設計單元中樣本大小不一致）
table(class$Sex,class$Age)
#獲得各單元的均值和標準差
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()，顯示計算結(jié)果
source("anova.tab.R");anova.tab(class.aov)

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

2.考慮交互作用

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

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

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")

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

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

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é)果，圖形對任意順序的因子設計的主效應和交互效應都會進行展示
library(HH)
interaction2wt(class$Weight ~ class$Sex*class$Age)

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

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

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

對于存在交互作用的兩因素，我們應當固定一個因素的水平，對另一個因素的水平進行水平間差異檢驗?
library(effects)
effect("humidity",fruit.aov)

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