方差分析(ANOVA)
什么是方差分析
模型可以归纳为在解决预测问题
Q：如何评价预测结果的好坏？
A：预测和模型越接近越好
Q：如何刻画”接近“？
A：需要定义”距离“，一般使用残差的平方和(sum of squares)
总的差异=组内差异+组间差异
SST=SSE + SSM
SSE越小，SSM越大，说明模型预测得越好
将SSM所占SST的比例理解为方差被解释的程度，记为R方。
但R方无法直接回答模型是否统计显著这个问题；因为堆砌参数会不断地增加R方(过拟合)
因此我们引入平方和的平均值。
R语言中的实现
data("mtcars")
mod1 <- lm(mpg~wt,data = mtcars)
anova(mod1)
两个模型的比较
之前的ANOVA可以看作是比较当前模型和只保留截距项的最简单模型，类似地可推广至比较全模型(full model)和简化模型(Reduced model)
mod2 <- lm(mpg ~ wt + am,data = mtcars)
anova(mod1,mod2)
单因素方差分析(one-way anova)
如果x(一维)是以类别变量的形式出现
则之前的模型检验可以看作为检验各组之间是否存在显著性差异
如果差异存在，我们应用各组均值分别进行预测；
如果差异不存在，我们应用全体均值进行统一预测。
那么此时问题又回到了方差分析。
注意此时系数的含义：与参数组的差值
mtcars$cyl <- factor(mtcars$cyl)
mod4 <- lm(mpg ~ cyl,data = mtcars)
summary(mod4)
也可以直接使用aov函数
mod4_anova <- aov(mpg~cyl,data = mtcars)
summary(mod4_anova)
均值的多重比较
anova仅能告诉我们各组之间是否存在差异
即便是只有一组与其他组不同，也是为存在差异，也就是说各组之间不一定都存在差异
我们自然关系究竟哪些组相同/不同
多重t检验
pairwise.t.test(mtcars$mpg,mtcars$cyl)
#对修正方式进行指定，详看帮助文档
pairwise.t.test(mtcars$mpg,mtcars$cyl,p.adjust.method = "none")
p值的修正
单纯地重复使用t检验会增加第一类错误地概率：H0正确，却拒绝了H0
需要对p value进行修正
一种标准为控制family-wise error rate
常用方法包括Bonferroni；Holm；scheffe，Tukey
如果我们需要对所有地两两组均进行检验，推荐使用Tukey
TukeyHSD只接受aov生成的对象
TukeyHSD(mod4_anova)
另一种标准为控制 False Discovery Rate(FDR)，一般针对组数较大的情况，只要求一定比例的检验结果是可靠的。
常见方法为BH。
方差齐性检验
对于类别型数据，同分布假设意味着各组方差相同，可以使用bartlett检验
但该检验对正态性假设比较敏感【如果不是正态假设，结果可能不靠谱】
bartlett.test(mpg~am,data = mtcars)
如果方差齐性假设被拒绝，则不应该使用简单的单因素方差分析
bartlett.test(mpg~cyl,data = mtcars)
不要求方差齐性的检验
oneway.test(mpg~cyl,data = mtcars)
不过前面的方法(比如多重比较)，也需要修正。
双因素方差分析(Two-way ANOVA)
考虑两个影响因素
mod5 <- aov(mpg~am+cyl,data = mtcars)
summary(mod5)
TukeyHSD(mod5)
上述模型包含可加性假设
即两个因素之间没有相互影响
如果认为一个因素的效应会受到另一因素的影响，则应该加入交互作用。
data(iris)
aov1 <- aov(Sepal.Length~Species, iris)
aov1
summary(aov1)
tukey <- TukeyHSD(aov1)
tukey = as.data.frame(tukey$Species)
tukey$pair = rownames(tukey)
一个例子
科研中常见的单因素方差分析的图还是不同分组的差异用字母abcd来表示，下面我以示例数据mtcars为例，结合ggplot，展示一下如何做一个方差分析的分组统计图。
rm(list = ls())
df <- data.frame(mtcars)
#进行方差分析
fit <- aov(cyl~mpg,data = df)
#查看一下是否有差异
summary(fit)
p value远小于0.01，所以认为，不同的分组之间有差异；
但我们想知道，具体那一组和那一组比较有差异，下面我将用到LSD.test。
#install.packages("agricolae")
library(agricolae)
res<-LSD.test(fit,"cyl",p.adj="none")
res
sta <- data.frame(res$groups)
plot(res)
从结果中我们看到，不同分组之间均有差异；如果要简易画图看趋势的话，直接用plot()函数，但我们想做的好看一点，就选用了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.https://www.jianshu.com/p/aa80b6f65399
