1. 基础知识
方差分析适用于多组均数的比较（在完全随机设计的实验中，两组均数的t检验和方差分析是完全等价的。但t检验只能用于两组的均数比较，对于三组和三组以上的均数比较，就需要用到方差分析。）
-
方差分析的目的：比较多组均数的差异-->分析一个数值型变量与一个(单因素)或多个(多因素)分类变量是否相关 -
方差分析的结果：不同组别的均数差异显著-->该数值型变量与分类变量是相关的 -
方差分析需要满足的条件：
1. 在控制变量的不同水平下，观测变量服从正态分布
2. 在控制变量的不同水平下，观测变量方差齐性
3. 各组数据相互独立 -
方差分析基本原理：
1. 实验条件，即不同的处理造成的差异，称为组间差异。用变量在各组的均值与总均值之偏差平方和的总和表示，记做SSb，组间自由度dfb。
2. 随机误差，如测量误差造成的差异或个体间的差异，称为组内差异，用变量在各组的均值与该组内变量值之偏差平方和的总和表示，记做SSw，组间自由度dfw。
总偏差平方和SSt=SSb+SSw
整个方差分析的基本步骤如下：
1、建立检验假设；
H0：多个样本总体均值相等；
H1：多个样本总体均值不相等或不全等。
检验水准为0.05。
2、计算检验统计量F值；
3、确定P值并作出推断结果。
2. 正态性检验和方差齐性检验
# 载入演示数据集
# 以multcomp包中的cholesterol数据集为例
library('multcomp')
data('cholesterol')
head(cholesterol)
# trt response
# 1 1time 3.8612
# 2 1time 10.3868
# 3 1time 5.9059
# 4 1time 3.0609
# 5 1time 7.7204
# 6 1time 2.7139
str(cholesterol)
# 'data.frame': 50 obs. of 2 variables:
# $ trt : Factor w/ 5 levels "1time","2times",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ response: num 3.86 10.39 5.91 3.06 7.72 ...
可以看到这个数据集只有两个变量，其中治疗是分类变量（因子型），有5个水平。response是数值型变量。要对每种治疗所对应的response的均值进行比较，就只能用方差分析而不能用t检验。
2.1 进行正态性检验
shapiro.test(cholesterol$response)
# Shapiro-Wilk normality test
# data: cholesterol$response
# W = 0.97722, p-value = 0.4417
符合正态分布
2.2 进行方差齐性检验（如下四个函数都可进行多组数据的方差齐性检验，t检验中提到的var.test函数只用于两组数据的方差齐性检验）
- a.
bartlett.test函数
bartlett.test(response~trt,data = cholesterol)
# Bartlett test of homogeneity of variances
# data: response by trt
# Bartlett's K-squared = 0.57975, df = 4, p-value =
# 0.9653
要比较均值的数据写～左边，分组变量写右边。p=0.9653，方差齐。
- b.
flinger.test函数
fligner.test(response~trt,data = cholesterol)
# Fligner-Killeen test of homogeneity of variances
# data: response by trt
# Fligner-Killeen:med chi-squared = 0.74277, df = 4,
# p-value = 0.946
写法同上，方差齐。
- c. car包中的
ncvTest函数
library(car)
ncvTest(lm(response~trt,data = cholesterol)) #传入的是线性模型
# Non-constant Variance Score Test
# Variance formula: ~ fitted.values
# Chisquare = 0.1338923, Df = 1, p = 0.71443
- d. car包中的
leveneTest函数
leveneTest(response~trt,data = cholesterol)
# Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
# group 4 0.0755 0.9893
# 45
需要注意的是，如果检验出方差不齐，我们第一步不是立马选择进行非参数检验，而是首先要判断有无异常值存在，因为异常值对方差的影响很大。当然，到这一步才来检验有无异常值是不符合数据分析的流程的，异常值在进行数据初步处理的时候就因该被发现和处理掉。
3 进行方差分析
方差分析包括单因素方差分析，多因素方差分析，协方差分析，多元方差分析，重复测量数据方差分析。
3.1单因素方差分析
-
aov函数
library(multcomp)
data("cholesterol")
head(cholesterol)
fit <- aov(response~trt,data = cholesterol)
summary(fit)
# Df Sum Sq Mean Sq F value Pr(>F)
# trt 4 1351.4 337.8 32.43 9.82e-13 ***
# Residuals 45 468.8 10.4
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##Df是自由度，Sum Sq是总的离差平方和，Mean Sq是平均的离差平方和，F value是检验统计量F值
##p值=0.01364，小于0.05，拒绝原假设，接受备择假设，response均值有差异
gplots包的plotmeans函数 对上述结果进行可视化
plotmeans(response~trt,data = cholesterol)
-
oneway.test函数
oneway.test(response~trt,data = cholesterol)
# One-way analysis of means (not assuming equal variances)
# data: response and trt
# F = 30.709, num df = 4.00, denom df = 22.46, p-value = 8.227e-09
oneway.test(response~trt,data = cholesterol,var.equal = T) #⚠️var.equal参数是对方差是否齐进行设置
# One-way analysis of means
# data: response and trt
# F = 32.433, num df = 4, denom df = 45, p-value = 9.819e-13
fit <- aov(response~trt,data = cholesterol,var.equal = T) #和aov函数得到的结果是一样的
summary(fit)
# Df Sum Sq Mean Sq F value Pr(>F)
# trt 4 1351.4 337.8 32.43 9.82e-13 ***
# Residuals 45 468.8 10.4
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
3.2 两因素方差分析
使用ToothGrowth数据集进行演示
data('ToothGrowth')
head(ToothGrowth)
# len supp dose
# 1 4.2 VC 0.5
# 2 11.5 VC 0.5
# 3 7.3 VC 0.5
# 4 5.8 VC 0.5
# 5 6.4 VC 0.5
# 6 10.0 VC 0.5
levels(ToothGrowth$supp)
# [1] "OJ" "VC"
levels(as.factor(ToothGrowth$dose))
# [1] "0.5" "1" "2"
table(as.factor(ToothGrowth$dose))
# 0.5 1 2
# 20 20 20
aov函数
fangcha <- aov(len~supp+dose,data = ToothGrowth)
summary(fangcha)
# Df Sum Sq Mean Sq F value Pr(>F)
# supp 1 205.4 205.4 11.45 0.0013 **
# dose 1 2224.3 2224.3 123.99 6.31e-16 ***
# Residuals 57 1022.6 17.9
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
不考虑supp和dose之间的交互作用的情况。结果显示两个因素都对小鼠牙齿生长影响显著。
fangcha <- aov(len~supp*dose,data = ToothGrowth)
summary(fangcha)
# Df Sum Sq Mean Sq F value Pr(>F)
# supp 1 205.4 205.4 12.317 0.000894 ***
# dose 1 2224.3 2224.3 133.415 < 2e-16 ***
# supp:dose 1 88.9 88.9 5.333 0.024631 *
# Residuals 56 933.6 16.7
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
考虑两个因素之间的交互作用：将上面的+换成*。结果显示两个因素都对小鼠牙齿生长影响显著而且两者间的互相影响也不容忽视。
可视化
interaction.plot(ToothGrowth$dose,ToothGrowth$supp,ToothGrowth$len,type = 'b',pch = c(1,10),col=c('red','green'))
3.3 两两比较（多组比较）的方差分析
上述结果已经知道了再五组数据中的均值不全相等，下一步想知道哪些相等哪些不等，就要对这五组进行两两比较。
思考🤔：这里的两两比较为什么不能用t检验？
背景假设：共有4组，进行两两比较的话，则一共需要进行6次，设α=0.05，6次比较都不发生假阳性错误的概率是(1-α)6=0.735
-
TukeyHSD函数
TukeyHSD(fit) #直接使用前面得到的fit
# Tukey multiple comparisons of means
# 95% family-wise confidence level
# Fit: aov(formula = response ~ trt, data = cholesterol, var.equal = T)
# $trt
# diff lwr upr p adj
# 2times-1time 3.44300 -0.6582817 7.544282 0.1380949
# 4times-1time 6.59281 2.4915283 10.694092 0.0003542
# drugD-1time 9.57920 5.4779183 13.680482 0.0000003
# drugE-1time 15.16555 11.0642683 19.266832 0.0000000
# 4times-2times 3.14981 -0.9514717 7.251092 0.2050382
# drugD-2times 6.13620 2.0349183 10.237482 0.0009611
# drugE-2times 11.72255 7.6212683 15.823832 0.0000000
# drugD-4times 2.98639 -1.1148917 7.087672 0.2512446
# drugE-4times 8.57274 4.4714583 12.674022 0.0000037
# drugE-drugD 5.58635 1.4850683 9.687632 0.0030633
输出的结果：从左往右依次是：两两比较、两两间的差值、lwr是95%可信线的下限，upr是上限。最后是p值。
将结果可视化：
plot(TukeyHSD(fit))
线段中点是均值，两端是95%置信区间，跨过0说明没有显著差异。
3.4 单因素协方差分析
在进行方差分析时，所有混杂因素统称为协变量。
协方差分析的前提条件：
1. 协变量是连续型变量
2. 协变量与Y存在线性关系。且在X的不同水平下，斜率大致相同，仅有截距不同。
# 使用litter数据集进行演示
head(litter)
# dose weight gesttime number
# 1 0 28.05 22.5 15
# 2 0 33.33 22.5 14
# 3 0 36.37 22.0 14
# 4 0 35.52 22.0 13
# 5 0 36.77 21.5 15
# 6 0 29.60 23.0 5
检验dose对weight的影响。出生时间gesttime是协变量。
- Step 1: 协变量与因变量之间的线性关系是否斜率相同？
# 使用aov函数查看
facn <- aov(weight~gesttime + dose + gesttime:dose,data = litter)
summary(facn)
# Df Sum Sq Mean Sq F value Pr(>F)
# gesttime 1 134.3 134.30 8.289 0.00537 **
# dose 3 137.1 45.71 2.821 0.04556 *
# gesttime:dose 3 81.9 27.29 1.684 0.17889
# Residuals 66 1069.4 16.20
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
aov后面小括号里写的顺序：结果变量～协变量+自变量。如果要看协变量和自变量之间是否存在交互，在后面写+协变量：自变量。最后是data=数据集。
结果显示两个变量之间不存在交互效应(p=0.17889, >0.05)，可以认为它们的斜率是相同的。
- Step2: 因变量的正态性和方差齐性
library(car)
qqPlot(lm(weight~gesttime+dose, data = litter), simulate=TRUE,
main='QQ Plot',labels=FALSE)
- Step3: 方差分析
fangcha=aov(weight~gesttime+dose,data = litter)
summary(fangcha)
# Df Sum Sq Mean Sq F value Pr(>F)
# gesttime 1 134.3 134.30 8.049 0.00597 **
# dose 3 137.1 45.71 2.739 0.04988 *
# Residuals 69 1151.3 16.69
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- Step4: 使用effects包里的effect函数去除协变量的影响
library(effects)
effect('dose',facn)
# NOTE: dose is not a high-order term in the model
# dose effect
# dose
# 0 5 50 500
# 32.30722 28.50625 30.61078 29.29005
3.5 多元方差分析MANOVA
因变量不止一个，但是需要将它们作为一个整体同时进行分析。例如：某种药物对患者血红蛋白浓度，红细胞计数，外周血细胞因子水平等多种因素的影响。
一元方差分析不能满足，必须使用多元方差分析的情况：
1. 因变量之间存在暧昧关系（性质类似或存在相关性）
2. 有时候单个变量难以反应出组间真实差异
3. 自变量为分类变量
使用manovs()函数进行性多元方差分析
library(MASS)
attach(UScereal)
y <- cbind(calories,fat,sugars)
fit_m <- manova(y~shelf)
summary(fit_m)
# Df Pillai approx F num Df den Df Pr(>F)
# shelf 1 0.19594 4.955 3 61 0.00383
# Residuals 63
# shelf **
# Residuals
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
3.6 重复测量方差分析 REPEAT ANOVA
基本特点:
1. 同一受试者在不同时间点上进行多次测量
2. 所测指标大多为连续型变量
3. 处理因素：g个水平（g>=1），每个水平下有n个受试者，共计g*n个实验对象
people <- rep(1:20,each=3)
time <- rep(LETTERS[1:3],20)
value <- sample(60:180,60)
dat <- data.frame(peo=people,time=time,value=value)
DT::datatable(dat)
fit1 <- aov(value~time,data=dat)
summary(fit1)
# Df Sum Sq Mean Sq F value Pr(>F)
# time 2 3810 1905 1.865 0.164
# Residuals 57 58213 1021
#其本质就是随机区组方差分析，只不过时间是控制变量
4. 拓展：
参考：https://blog.csdn.net/dingming001/article/details/72822270
