方差分析(ANOVA)
什么是方差分析
模型可以歸納為在解決預(yù)測(cè)問題
Q：如何評(píng)價(jià)預(yù)測(cè)結(jié)果的好壞漱贱？
A：預(yù)測(cè)和模型越接近越好
Q：如何刻畫”接近“蔗包？
A：需要定義”距離“蛛壳，一般使用殘差的平方和(sum of squares)
總的差異=組內(nèi)差異+組間差異
SST=SSE + SSM
SSE越小寞缝，SSM越大派敷，說明模型預(yù)測(cè)得越好
將SSM所占SST的比例理解為方差被解釋的程度蛹批，記為R方撰洗。
但R方無法直接回答模型是否統(tǒng)計(jì)顯著這個(gè)問題；因?yàn)槎哑鰠?shù)會(huì)不斷地增加R方(過擬合)
因此我們引入平方和的平均值腐芍。
R語言中的實(shí)現(xiàn)
data("mtcars")
mod1 <- lm(mpg~wt,data = mtcars)
anova(mod1)
兩個(gè)模型的比較
之前的ANOVA可以看作是比較當(dāng)前模型和只保留截距項(xiàng)的最簡(jiǎn)單模型差导，類似地可推廣至比較全模型(full model)和簡(jiǎn)化模型(Reduced model)
mod2 <- lm(mpg ~ wt + am,data = mtcars)
anova(mod1,mod2)
單因素方差分析(one-way anova)
如果x(一維)是以類別變量的形式出現(xiàn)
則之前的模型檢驗(yàn)可以看作為檢驗(yàn)各組之間是否存在顯著性差異
如果差異存在，我們應(yīng)用各組均值分別進(jìn)行預(yù)測(cè)猪勇；
如果差異不存在设褐，我們應(yīng)用全體均值進(jìn)行統(tǒng)一預(yù)測(cè)。
那么此時(shí)問題又回到了方差分析泣刹。
注意此時(shí)系數(shù)的含義：與參數(shù)組的差值
mtcars$cyl <- factor(mtcars$cyl)
mod4 <- lm(mpg ~ cyl,data = mtcars)
summary(mod4)
也可以直接使用aov函數(shù)
mod4_anova <- aov(mpg~cyl,data = mtcars)
summary(mod4_anova)
均值的多重比較
anova僅能告訴我們各組之間是否存在差異
即便是只有一組與其他組不同助析，也是為存在差異，也就是說各組之間不一定都存在差異
我們自然關(guān)系究竟哪些組相同/不同
多重t檢驗(yàn)
pairwise.t.test(mtcars$mpg,mtcars$cyl)
#對(duì)修正方式進(jìn)行指定项玛，詳看幫助文檔
pairwise.t.test(mtcars$mpg,mtcars$cyl,p.adjust.method = "none")
p值的修正
單純地重復(fù)使用t檢驗(yàn)會(huì)增加第一類錯(cuò)誤地概率：H0正確貌笨，卻拒絕了H0
需要對(duì)p value進(jìn)行修正
一種標(biāo)準(zhǔn)為控制family-wise error rate
常用方法包括Bonferroni；Holm襟沮；scheffe锥惋，Tukey
如果我們需要對(duì)所有地兩兩組均進(jìn)行檢驗(yàn)，推薦使用Tukey
TukeyHSD只接受aov生成的對(duì)象
TukeyHSD(mod4_anova)
另一種標(biāo)準(zhǔn)為控制 False Discovery Rate(FDR)开伏，一般針對(duì)組數(shù)較大的情況膀跌，只要求一定比例的檢驗(yàn)結(jié)果是可靠的。
常見方法為BH固灵。
方差齊性檢驗(yàn)
對(duì)于類別型數(shù)據(jù)捅伤，同分布假設(shè)意味著各組方差相同，可以使用bartlett檢驗(yàn)
但該檢驗(yàn)對(duì)正態(tài)性假設(shè)比較敏感【如果不是正態(tài)假設(shè)巫玻，結(jié)果可能不靠譜】
bartlett.test(mpg~am,data = mtcars)
如果方差齊性假設(shè)被拒絕丛忆，則不應(yīng)該使用簡(jiǎn)單的單因素方差分析
bartlett.test(mpg~cyl,data = mtcars)
不要求方差齊性的檢驗(yàn)
oneway.test(mpg~cyl,data = mtcars)
不過前面的方法(比如多重比較)，也需要修正仍秤。
雙因素方差分析(Two-way ANOVA)
考慮兩個(gè)影響因素
mod5 <- aov(mpg~am+cyl,data = mtcars)
summary(mod5)
TukeyHSD(mod5)
上述模型包含可加性假設(shè)
即兩個(gè)因素之間沒有相互影響
如果認(rèn)為一個(gè)因素的效應(yīng)會(huì)受到另一因素的影響熄诡，則應(yīng)該加入交互作用。
data(iris)
aov1 <- aov(Sepal.Length~Species, iris)
aov1
summary(aov1)
tukey <- TukeyHSD(aov1)
tukey = as.data.frame(tukey$Species)
tukey$pair = rownames(tukey)
一個(gè)例子
????科研中常見的單因素方差分析的圖還是不同分組的差異用字母abcd來表示诗力，下面我以示例數(shù)據(jù)mtcars為例凰浮，結(jié)合ggplot，展示一下如何做一個(gè)方差分析的分組統(tǒng)計(jì)圖苇本。
rm(list = ls())
df <- data.frame(mtcars)
#進(jìn)行方差分析
fit <- aov(cyl~mpg,data = df)
#查看一下是否有差異
summary(fit)
p value遠(yuǎn)小于0.01袜茧，所以認(rèn)為，不同的分組之間有差異瓣窄；
但我們想知道笛厦，具體那一組和那一組比較有差異，下面我將用到LSD.test康栈。
#install.packages("agricolae")
library(agricolae)
res<-LSD.test(fit,"cyl",p.adj="none")
res
sta <- data.frame(res$groups)
plot(res)
????從結(jié)果中我們看到递递，不同分組之間均有差異喷橙；如果要簡(jiǎn)易畫圖看趨勢(shì)的話，直接用plot()函數(shù)登舞，但我們想做的好看一點(diǎn)贰逾，就選用了ggplot2。
library(ggplot2)
ggplot(df,aes(factor(cyl),mpg,group=cyl))+
geom_boxplot()+
geom_jitter(aes(fill=mpg),width=0.2) +
annotate("text",label="c",x=1.2,y=32) +
annotate("text",label="b",x=2.2,y=23) +
annotate("text",label="a",x=3.2,y=20) +
theme_bw()+
theme(panel.grid = element_blank(),
legend.position = "none")+
xlab("")+
ylab("")
效果圖：
參考鏈接：
1.https://blog.csdn.net/yijiaobani/article/details/115494887
2.https://zhuanlan.zhihu.com/p/38412409
3.https://zhuanlan.zhihu.com/p/296726829
4.http://www.reibang.com/p/aa80b6f65399
