各类统计方法R语言实现（三）
各位小伙伴们大家好，今天是我们的系列推文“各类统计方法R语言实现”第三篇，今天介绍的主要内容是方差分析。
方差分析适用用于两个及两个以上样本均数差别的显著性检验，其中两组之间的方差分析等价于t检验，因此常用于两组以上样本均数比较。
方差分析也必须满足正态性与方差齐性，代码与上次推文类似。
以下主要介绍单因素方差分析，单因素协变量方差分析，重复测量方差分析。
data(mtcars)
head(mtcars)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
#计算描述统计分析
mtcars$cyl<-as.factor(mtcars$cyl)
mtcars$vs<-as.factor(mtcars$vs)
mtcars$am<-as.factor(mtcars$am)
mtcars$gear<-as.factor(mtcars$gear)
mtcars$carb<-as.factor(mtcars$carb)
#计算描述统计分析
summary(mtcars)
## mpg cyl disp hp drat
## Min. :10.40 4:11 Min. : 71.1 Min. : 52.0 Min. :2.760
## 1st Qu.:15.43 6: 7 1st Qu.:120.8 1st Qu.: 96.5 1st Qu.:3.080
## Median :19.20 8:14 Median :196.3 Median :123.0 Median :3.695
## Mean :20.09 Mean :230.7 Mean :146.7 Mean :3.597
## 3rd Qu.:22.80 3rd Qu.:326.0 3rd Qu.:180.0 3rd Qu.:3.920
## Max. :33.90 Max. :472.0 Max. :335.0 Max. :4.930
## wt qsec vs am gear carb
## Min. :1.513 Min. :14.50 0:18 0:19 3:15 1: 7
## 1st Qu.:2.581 1st Qu.:16.89 1:14 1:13 4:12 2:10
## Median :3.325 Median :17.71 5: 5 3: 3
## Mean :3.217 Mean :17.85 4:10
## 3rd Qu.:3.610 3rd Qu.:18.90 6: 1
## Max. :5.424 Max. :22.90 8: 1
正态性检验
使用Shapiro-Wilk法
shapiro.test(mtcars$mpg[mtcars$cyl == 4])
##
## Shapiro-Wilk normality test
##
## data: mtcars$mpg[mtcars$cyl == 4]
## W = 0.91244, p-value = 0.2606
shapiro.test(mtcars$mpg[mtcars$cyl == 6])
##
## Shapiro-Wilk normality test
##
## data: mtcars$mpg[mtcars$cyl == 6]
## W = 0.89903, p-value = 0.3252
shapiro.test(mtcars$mpg[mtcars$cyl == 8])
##
## Shapiro-Wilk normality test
##
## data: mtcars$mpg[mtcars$cyl == 8]
## W = 0.93175, p-value = 0.3229
可以看到结果中各组p均大于0.05，服从正态分布。
方差齐性检验
方差齐性检验的方法有很多，本节主要依旧介绍以下三种：
Bartlett检验：需要数据服从正态分布
Levene检验：不依赖总体分布具体形式，更为稳健。
Fligner-Killeen检验：不依赖总体分布具体形式。
##绘制箱线图，
plot(mpg~cyl,data = mtcars)
##Bartlett检验
bartlett.test(mpg~cyl,data = mtcars)
##
## Bartlett test of homogeneity of variances
##
## data: mpg by cyl
## Bartlett's K-squared = 8.3934, df = 2, p-value = 0.01505
可以看到结果中p<0.05，方差不齐。
library(car)
## Loading required package: carData
##Levene检验
leveneTest(mpg~cyl,data = mtcars)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 5.5071 0.00939 **
## 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
可以看到结果中p<0.05，方差不齐。
##Fligner-Killeen检验
fligner.test(mpg~cyl,data = mtcars)
##
## Fligner-Killeen test of homogeneity of variances
##
## data: mpg by cyl
## Fligner-Killeen:med chi-squared = 6.8113, df = 2, p-value = 0.03319
可以看到结果中p<0.05，方差不齐。
对于多组数据，正态分布且方差齐的资料，可以使用单因素方差分析；非正态分布或方差不齐，可以进行变量转换或者适用K-W检验。
单因素方差分析
从以上结果可知，该数据服从正态分布，但方差不齐，说明不适用单因素方差分析，以下代码仅为演示。
##单因素方差分析
fit<-aov(mpg~cyl,data = mtcars)
summary(fit)
## Df Sum Sq Mean Sq F value Pr(>F)
## cyl 2 824.8 412.4 39.7 4.98e-09 ***
## Residuals 29 301.3 10.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p<0.05，表示各组均数不全相等，若要比较任意两组均数，则要进行多重比较。
##多重比较，Tukey方法
TukeyHSD(fit, conf.level = 0.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mpg ~ cyl, data = mtcars)
##
## $cyl
## diff lwr upr p adj
## 6-4 -6.920779 -10.769350 -3.0722086 0.0003424
## 8-4 -11.563636 -14.770779 -8.3564942 0.0000000
## 8-6 -4.642857 -8.327583 -0.9581313 0.0112287
plot(TukeyHSD(fit, conf.level = 0.95),xlim=c(-16,2))
p<0.05，结果均有统计学差异。
根据图片判断，未越过虚线则表示有统计学差异。
图片与计算结果等价。
##a使用multcomp多重比较，Tukey方法
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.6.3
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
tuk<-glht(fit,linfct=mcp(cyl="Tukey"))
summary(tuk)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = mpg ~ cyl, data = mtcars)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 6 - 4 == 0 -6.921 1.558 -4.441 <0.001 ***
## 8 - 4 == 0 -11.564 1.299 -8.905 <0.001 ***
## 8 - 6 == 0 -4.643 1.492 -3.112 0.011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
与以上结果一致
单因素协方差分析
有时我们的研究涉及协变量的干扰，导致研究的因变量出现改变，此时可使用单因素协方差分析矫正协变量影响,此处将hp当作协变量，分析代码如下：
##单因素协方差分析
fit<-aov(mpg~hp+cyl,data = mtcars)
summary(fit)
## Df Sum Sq Mean Sq F value Pr(>F)
## hp 1 678.4 678.4 68.530 5.22e-09 ***
## cyl 2 170.5 85.3 8.612 0.00122 **
## Residuals 28 277.2 9.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
得到的结果表明：
（1）hp与mpg相关。
（2）控制hp，cyl与mpg相关
##计算协变量调整后的均值
library(effects)
## Warning: package 'effects' was built under R version 3.6.3
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
effect('cyl', fit)
##
## cyl effect
## cyl
## 4 6 8
## 25.12392 19.15627 16.60307
单因素协方差分析的假设条件
除了服从正态分布与方差齐之外，单因素协方差分析还假定回归斜率相等。
##检验斜率是否相等
fit2<-aov(mpg~hp*cyl,data = mtcars)
summary(fit2)
## Df Sum Sq Mean Sq F value Pr(>F)
## hp 1 678.4 678.4 73.964 4.44e-09 ***
## cyl 2 170.5 85.3 9.295 0.000901 ***
## hp:cyl 2 38.7 19.4 2.110 0.141533
## Residuals 26 238.5 9.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
可以看到交互效应不显著（即hp:cyl p>0.05），支持斜率相等的假设。
若假设不成立，可以转化数据或使用sm包中的sm.ancova()函数
重复测量方差分析
重复测量方差分析用于受试者被测量不止一次的数据，比如在细胞培养中第1，3，5，7天分别测试细胞的增殖情况，或者在治疗1，3，5，7天抽取患者外周血用于检测血中的某些指标。
此处使用R中自带的CO2数据集。此处因变量是二氧化碳吸收量（uptake），自变量是植物类型Type（魁北克VS.密西西比）和七种水平的二氧化碳浓度（conc）。Type是组间因子，conc是组内因子。（此案例选自《R语言实战》，略作修改）
data(CO2)
w1b1 <- subset(CO2, Treatment == 'chilled') #只选择其中的寒带植物
#将分组转变为因子类型
w1b1$Type <- factor(w1b1$Type)
w1b1$conc <- factor(w1b1$conc)
#重复测量方差分析
fit <- aov(uptake ~ conc*Type + Error(Plant/conc), data = w1b1)
summary(fit)
##
## Error: Plant
## Df Sum Sq Mean Sq F value Pr(>F)
## Type 1 2667.2 2667.2 60.41 0.00148 **
## Residuals 4 176.6 44.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Plant:conc
## Df Sum Sq Mean Sq F value Pr(>F)
## conc 6 1472.4 245.40 52.52 1.26e-12 ***
## conc:Type 6 428.8 71.47 15.30 3.75e-07 ***
## Residuals 24 112.1 4.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
结果显示，类型和浓度以及交叉效应“类型×浓度”都非常显著。
par(las=2)
par(mar=c(10,4,4,2))
#法1：点线图
with(w1b1,
interaction.plot(conc,Type,uptake,
type="b", col=c("red","blue"), pch=c(16,18),
main="Interaction Plot for Plant Type and Concentration"))
#法2：箱线图
boxplot(uptake ~ Type*conc, data=w1b1, col=(c("gold","green")),
main="Chilled Quebec and Mississippi Plants",
ylab="Carbon dioxide uptake rate (umol/m^2 sec)")
结合上图及分析结果可得结论：魁北克省（箱线图黄色）植物比密西西比（箱线图绿色）植物的二氧化碳吸收率高，且随着CO2浓度升高，差异越来越明显。
好了，今天的R语言实现统计方法系列推文暂时告一段落，我们下次再见吧！小伙伴们如果有什么统计上的问题，或者如果想要学习什么方面的生物信息内容，可以在微信群或者知识星球提问，没准哪天的推文就是专门解答你的问题哦！