对于组数高于2组的,t检验不再适用,这时,方差分析适用
单因子方差分析(one-way ANOVA)We use the data,cholesterol,作为例子 完全随机设计的方差分析
library(multcomp)
library("cowplot")
data("cholesterol")
head(cholesterol)
建模:
1.正态性检验2.方差齐性检验3.进行方差分析步骤
aggregate(cholesteroltrt), FUN=mean)
plot(cholesteroltrt)
model=aov(cholesteroltrt)
summary(model)
多重比较 可用于two-way anova Turkey法用于各组数量一样 SNK-q法用于各组数量不等时的比较
方差分析只告诉我们这五组之间是不同的,但没有告诉我们哪两组之间有明显差别
此时需要使用TukeyHSD函数进行均值的多重比较分析,从结果中观察到有三个两两比较是不显著的。
result=TukeyHSD(model)
plot(result)
假设检验
方差分析需要一定的假设,即数据集应该符合正态和同方差,我们分别用下面的函数来进行检验,从P值观察到这两个假设是符合的。
对于不符合假设的情况,我们就要用到非参数方法,例如Kruskal-Wallis秩和检验 ##秩和检验,没听说过用
two-way anova 随机区间设计的方差分析
如果问两种喂食方法对豚鼠牙齿长度的影响是否有差异,则用t检验。原因是只有一个自变量即喂食方法,且该自变量包含2个水平;
如果问单一喂食方法中三种水平的抗坏血酸含量对豚鼠牙齿长度的影响是否有差异,则用单因素方差分析。原因是只有一个自变量(单因素)即抗坏血酸含量,且该自变量包含3个水平;
再如果我们对喂食方式不感兴趣,但它又会对实验结果造成影响,那么可以通过区组设计来控制,这便是一个完全随机区组设计Randomized Complete Block Design (RCBD)
回到正题,该实验设计包含了2个因素(factor),第一个因素包含2个水平(level),第二个因素包含3个水平,这是一个2*3设计。采用的数据分析方法为双因素方差分析。
我们仍然使用aov进行建模,然后使用HH包的绘图函数来展现双因子交互效果图
data("ToothGrowth")
head(ToothGrowth)
str(ToothGrowth)
library("dplyr")
ToothGrowth %>% group_by(supp, dose) %>% summarise(len_mean = mean(len))
ToothGrowthdose)
dose目前还是双精度浮点数(dbl)型变量,由于它是一个分类变量,所以需要将其改为因子型(factor)
attach(ToothGrowth) # 将ToothGrowth变成R默认处理数据集
fit <- aov(len ~ supp*dose)
那么为什么我们用*而不是+呢?
这是因为*表示所有可能的交互项。两因素之间可能有交互作用,如果用+,系统便不会计算supp与dose的交互作用。
summary(fit)
interaction.plot(dose, supp, len,
type = "b",
col = c("#FF6347", "#00CD66"),
pch = c(2,19),
xlab = "抗坏血酸含量",
ylab = "牙齿长度",
main = "抗坏血酸含量和喂食方式的交互作用")
interaction.plot(x.factor, trace.factor, response,
type = c("l", "p", "b", "o", "c"),
# col = 1,
# pch = c(1:9, 0, letters),
# main = " ",
# xlab = " ",
# ylab = "牙齿长度",
# ...)
x.factor:自变量,即因素,但记住这个变量将形成x轴;
trace.factor:两因素中的另一个因素;
response:因变量;
type:线(l)、点(p)、线和点(b)等;
col:添加颜色;
pch:绘制点符号类型;
main:标题;
xlab:x轴标题;
ylab:y轴标题。
ggplots包中的plotmeans()可以用来绘制交互效应,而且可以展示置信区间和样本大小
library(gplots)
plotmeans(len ~ interaction(supp, dose, sep = " "),
connect = list(c(1,3,5),c(2,4,6)),
col = c("#FF6347", "#00CD66"),
main = "抗坏血酸含量和喂食方式的交互作用",
xlab = "处理组",
ylab = "牙齿长度")
该图相对interaction.plot()的图来说,多了置信区间以及样本量的显示。
最后,HH包中的interaction2wt()可以同时展示主效应和交互效应
library(HH)
interaction2wt(len ~ supp*dose)