R语言方差分析ANOVA

自己整理编写的R语言常用数据分析模型的模板，原文件为Rmd格式，直接复制粘贴过来，作为个人学习笔记保存和分享。
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)
#各组的均值
aggregate(medicine$Response,by=list(medicine$Treatment),FUN=mean)
#各组的标准差
aggregate(medicine$Response,by=list(medicine$Treatment),FUN=sd)
#调用aov函数进行方差分析（检验组间差异）
medicine.aov <- aov(Response ~ Treatment,data=medicine)
#summary提取方差分析的结果
summary(medicine.aov)

分析上述计算结果，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")