R語言方差分析ANOVA
自己整理編寫的R語言常用資料分析模型的模板，原檔案為Rmd格式，直接複製貼上過來，作為個人學習筆記儲存和分享。部分參考薛毅的《統計建模與R軟體》和《R語言實戰》
I. 單因素方差分析
#用data frame的格式輸入資料
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)
#各組的均值 分析上述計算結果，Df表示自由度，Sum Sq 表示平方和，Mean Sq 表示均方，F value 是F值，Pr(>F)是p值，A即為因子A，Residuals 是殘差。
但是我們注意到，這個結果並不完整。直接用summary()函式時候，只有因素A和誤差兩行，沒有總和，這裡編個小程式(anova.tab.R)作改進，計算方法為：將summary函式得到表中的第一行與第二行求和，得到總和行的值。
#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函式儲存在工作目錄中。
getwd()
#利用anova.tab.R函式，得到完整的方差分析表
source("anova.tab.R");anova.tab(medicine.aov)
#畫圖
plot(medicine$Response~medicine$Treatment)
#繪製各組均值及其置信區間的圖形
library(gplots)
plotmeans(medicine$Response~medicine$Treatment,xlab = "Treatment",ylab = "Response",main = "Mean Plot\nwith 95% CI")1.多重比較
ANOVA對各療法的F檢驗表明，4種藥品用於緩解術後疼痛的療效不同，但是並不能得出哪種藥品療法與其他不同。多重比較可以解決這個問題.e.g. TukeyHSD()函式提供了對各組均值差異的成對檢驗；multcomp包中的glht()函式提供了多重均值比較更為全面的方法，既適用於線性模型，也適用於廣義線性模型；多重t檢驗方法針對每組資料進行t檢驗。程式碼如下：
TukeyHSD(medicine.aov)
#par()函式旋轉軸標籤，增大左邊介面積，使標籤擺放更美觀。
par(las = 2)
par(mar = c(5, 8, 4, 2))
plot(TukeyHSD(medicine.aov))圖形中置信區間包含0的藥品對比，說明差異不顯著。
library(multcomp)
#為適合字母陣列擺放，par語句用來增大頂部邊介面積
par(mar = c(5, 4, 6, 2))
tuk <- glht(medicine.aov, linfct = mcp(Treatment = "Tukey"))
#cld()函式中level選項為設定的顯著性水平（這裡的0.05對應95%置信區間）
plot(cld(tuk, level = 0.05), col = "lightgrey")有相同字母的組（用箱線圖表示）說明均值差異不顯著。
多次重複使用t檢驗會增大犯第一類錯誤的概率，為了克服這一缺點，需要調整p-值。R軟體調整p-值用的是p.adjust()函式，函式使用的不同引數代表不同的調整方法。
attach(medicine)
#求資料在各水平下的均值
mu<-c(mean(Response[Treatment==1]), mean(Response[Treatment==2]), mean(Response[Treatment==3]),mean(Response[Treatment==4])); mu
#作多重t檢驗。這裡用到的pairwise.t.test()函式用來得到多重比較的p值
pairwise.t.test(Response, Treatment, p.adjust.method = "none")
#觀察兩個作調整後的p值的情況。p.adjust.method()函式的引數也可換為"hochberg","hommel","bonferroni","BH","BY","fdr"等。
pairwise.t.test(Response, Treatment, p.adjust.method = "holm")
#繪製箱線圖
plot(medicine$Response~medicine$Treatment)從上述結果可見，124無顯著差異，3與124均有顯著差異，即緩解疼痛的4種藥品，3與124有顯著差異，124間差異不顯著
2.評估檢驗的假設條件
擬合結果的可信度來源於，做統計檢驗時資料滿足假設條件的程度
(1)誤差的正態性檢驗
單因素方差分析中，我們假設因變數服從正態分佈，各組方差相等。可用Q-Q圖來檢驗正態性假設。擬合診斷如下：
library(car)
qqPlot(lm(Response ~ Treatment, data = medicine), simulate = TRUE,
main = "QQ Plot", labels = FALSE)資料幾乎都落在95%的置信區間範圍內，說明滿足正態性假設
也可用W檢驗（shapiro.test()函式）方法對資料作正態性檢驗
attach(medicine)
shapiro.test(Response[Treatment==1])
shapiro.test(Response[Treatment==2])
shapiro.test(Response[Treatment==3])
shapiro.test(Response[Treatment==4])計算結果表明，資料在四種水平下的均是正態的
(2)方差的其次性檢驗
方差的其次性檢驗就是檢驗資料在不同的水平下方差是否相同，常用方法是Bartlett檢驗
#R裡用bartlett.test()函式來提供Bartlett檢驗。另外還有Fligner-Killeen檢驗（fligner.test()函式)和Brown-Forsythe檢驗（HH包中的hov()函式）
bartlett.test(Response~Treatment,data=medicine)p值0.1285>0.05，接受原假設，認為各組的資料是等方差的
方差其次性分析對離群點非常敏感，可用car包的outlierTest()函式來檢測離群點
library(car)
outlierTest(medicine.aov)從p值結果看，並沒有證據說明該資料中含有離群點
根據Q-Q圖，Bartlett檢驗和離群點檢驗，該資料似乎可以用ANOVA模型擬合得很好，這些方法反過來增強了我們對於所得結果的信心
資料的總體分佈型別未知；或資料的總體分佈型別已知，但不符合正態分佈；或某些變數可能無法精確測量時，可以使用非引數統計方法.非引數統計是拋開總體分佈型別不考慮，對總體引數不做比較，比較的是總體分佈的位置是否相同的統計方法．秩和檢驗是非引數統計中一種經常使用的檢驗方法．這裡的“秩”又可被稱為等級，即按照資料大小排定的次序號．此次序號的總和被稱為“秩和”．
方差分析過程需要滿足若干條件，F檢驗才能奏效，可惜有時候採集到的資料並不能滿足這樣的要求。像兩樣本比較時一樣，嘗試將資料轉換為秩統計量，因為秩統計量的分佈與總體分佈無關，這樣就可以避開總體分佈的要求．上述問題就可以通過資料的秩統計量就解決了。在比較兩個以上的總體時，廣泛使用的是Kruskal-Wallis秩和檢驗，它是對兩個以上樣本進行比較的非引數檢驗方法。實質上，它是兩樣本的Wilcoxon方法在多於兩個樣本時的推廣。
R軟體提供了Kruskal-Wallis秩和檢驗，函式為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，拒絕原假設，認為四種藥物緩解疼痛效果有顯著差異
該函式還有另外兩種寫法如下：
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))之後再對上述資料作正太檢驗和方差齊次檢驗，如果全部通過檢驗，則該資料也可以作方差分析
(4)Friedman秩和檢驗
在配伍組設計中，多個樣本的比較，如果它們的總體不能滿足正態性和方差齊性的要求，可採用Friedman秩和檢驗
Friedman秩和檢驗的基本思想與前面介紹的方法類似，但是配伍組設計的隨機化是在配伍組內進行的，而配伍組間沒有進行隨機化。因此在進行Friedman秩和檢驗時，是分別在每個配伍組裡將資料從小到大編秩，如果相同的資料取平均秩次。
R軟體中，函式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)該函式還有另外兩種寫法如下：
x <- c(7,5,3,1,6,5,3,3,7,9,9,9,4,3,4,3)
#4行4列，每行4個數據，總共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.單因素協方差分析（顯著因素下的水平間差異檢驗）
單因素協方差分析（ANCOVA）擴充套件了單因素方差分析（ANOVA），包含一個或多個定量的協變數。下面的例子來自於multcomp包中的litter資料集，懷孕小鼠被分為四個小組，每個小組接受不同劑量（0、5、50、500）的藥物處理，產下幼崽的體重均值為因變數，懷孕時間為協變數。
（1）單因素ANCOVA
data(litter, package = "multcomp")
attach(litter)
#table()函式，看到每種劑量下所產幼崽數量並不同
table(dose)
#aggrgate()函式獲得各組均值，可以發現未用藥組幼崽體重均值最高
aggregate(weight, by = list(dose), FUN = mean)
fit <- aov(weight ~ gesttime + dose)
summary(fit)由於使用了協變數，如果想要獲取調整的組均值–即去除協變數效應後的組均值，可使用effects包中的effects()函式來計算調整的均值
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()函式新增。從結果來看，假設檢驗的t統計量在p<0.05水平下顯著，可以得出未用藥組比其他用藥條件下的出生體重高的結論
（3）評估檢驗的假設條件–檢驗同歸斜率的同質性
ANCOVA與ANOVA相同，都需要正態性和方差齊次性假設，可用上述ANOVA的假設檢驗的相同步驟來檢驗。另外ANCOVA還假定迴歸斜率相同。ANCOVA模型包含懷孕時間*劑量的互動項時，可對迴歸斜率的同質性進行檢驗。互動效應若顯著，則意味著時間和幼崽出生體重間的關係依賴於藥物劑量的水平
library(multcomp)
fit2 <- aov(weight ~ gesttime * dose)
summary(fit2)結果可以看到互動效應不顯著，支援了斜率相等的假設。若假設不成立，可以嘗試變換協變數或因變數，或使用能對每個斜率獨立解釋的模型，或使用不需要假設迴歸斜率同質性的非引數ANCOVA方法。（如sm包中的sm.ancova()函式）
（4）結果視覺化
HH包中的ancova()函式可以繪製因變數、協變數和因子之間的關係圖，程式碼如下：
library(HH)
ancova(weight ~ gesttime + dose, data = litter)從圖中可看出，用懷孕時間來預測出生體重的迴歸線相互平行，只是截距項不同。隨著懷孕時間增加，幼崽出生體重也會增加。另外，還可以看到0劑量組截距項最大，5劑量組截距項最小。由於之前的設定，直線會保持平行，若用anvova(weight~gesttime*dose),生成的圖形將允許斜率和截距項依據組別而發生變化，這對視覺化那些違揹回歸斜率同質性的例項非常有用
II.雙因素方差分析
1.不考慮互動作用
SAS自帶資料集sashelp.class中包含了學生的姓名、性別與身高。匯出資料存為csv格式，現在分析年齡與性別是否是影響體重的顯著因素。該問題屬於不均衡資料集的方差分析
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)
#呼叫自變函式anova.tab()，顯示計算結果
source("anova.tab.R");anova.tab(class.aov)根據p值不同說明年齡和性別對體重有顯著影響
2.考慮互動作用
（1）3種方式對結果進行視覺化處理
用interaction.plot()函式來展示雙因素方差分析的互動效應
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()函式來展示互動效應
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()函式來視覺化結果，圖形對任意順序的因子設計的主效應和互動效應都會進行展示
library(HH)
interaction2wt(class$Weight ~ class$Sex*class$Age)（2）有互動作用的方差分析
資料集fruit記錄了在不同溼度和溫度下某種植物的查處。這是一個雙因素方差分析的情形。假設方差分析的假設條件滿足，在顯著性水平0.05的前提下，欲分析不同溫度、不同溼度下產出是否有顯著差異，以及溫度和溼度的互動是否顯著差異，如果互動有差異，分析在溼度一定的情況下，溫度對產出的影響。
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:方差分析是一種常見的統計模型，用於檢驗樣本間均值是否相等。方差分析適用於處理因素型別為分類變數、響應變數型別為連續的情形。根據因素個數，方差分析可以分為單因素方差分析與多因素方差分析。在多因素方差分析中，要特別注意判斷因素間是否存在互動作用。此外，在實際應用中，可以通過設計合理的試驗，在儘可能排除外部因素的干擾後，再對試驗資料進行方差分析，這樣結果會更準確。
write.csv(medcine,"test_medcine.csv")
write.csv(class,"test_class.csv")