R语言学习笔记（七）方差分析

写在前面

总结一下最近学的方差分析的R语言实现。

方差分析(ANOVA)

方差分析是统计检验的一种。由英国著名统计学家：R.A.Fisher推导，也叫F检验。

用于多个样本间均数的比较(分析类别变量(有序变量))。当包含的因子是解释变量时，关注的重点通常会从预测转向组别差异的分析。

一些术语

1. 组间因子和组内因子(同时含有两种因子->双因素混合模型方差分析)。
2. 主效应、交互效应。
3. 单因素方差分析、双因素方差分析。
4. 协方差分析(存在协变量，单元、多元)。

• 注意：(表达式中)变量输入的顺序会影响模型！

单因素方差分析

library(multcomp)

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

dt <- cholesterol
attach(dt)
table(trt)

## trt
##  1time 2times 4times  drugD  drugE 
##     10     10     10     10     10

aggregate(response, by=list(trt), FUN=mean)

##   Group.1        x
## 1   1time  5.78197
## 2  2times  9.22497
## 3  4times 12.37478
## 4   drugD 15.36117
## 5   drugE 20.94752

fit <- aov(response~trt)
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

多重比较

TukeyHSD(fit)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = response ~ trt)
## 
## $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

# 绘图
par(las=2) # 标签垂直于轴
par(mar=c(5,8,4,2)) # 设置边距便于显示标签
plot(TukeyHSD(fit)) # 绘图，触碰虚线则无明显差异

评估检验的假设条件

严谨性的评估(判断)，前提假设是否满足。

1. 因变量$y$服从正态分布；
2. 各组之间方差齐性(相等)。

正态性检验：

qqPlot(lm(response~trt, data=dt))

## [1] 19 38

散点基本位于通道之间即可证明正态性。

方差齐性检验：

bartlett.test(response~trt, data=dt)

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  response by trt
## Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653

# 对离群点进行检验
outlierTest(fit)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 19 2.251149           0.029422           NA

双因素方差分析

双因素方差分析的实现

dt <- ToothGrowth
attach(dt)
# View(dt)
table(supp, dose)

##     dose
## supp 0.5  1  2
##   OJ  10 10 10
##   VC  10 10 10

# 查看每一组别的均值
aggregate(len, by=list(supp, dose), FUN = mean)

##   Group.1 Group.2     x
## 1      OJ     0.5 13.23
## 2      VC     0.5  7.98
## 3      OJ     1.0 22.70
## 4      VC     1.0 16.77
## 5      OJ     2.0 26.06
## 6      VC     2.0 26.14

# 进行方差分析
dose <- factor(dose)
fit <- aov(len~supp*dose)
summary(fit)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## supp         1  205.4   205.4  15.572 0.000231 ***
## dose         2 2426.4  1213.2  92.000  < 2e-16 ***
## supp:dose    2  108.3    54.2   4.107 0.021860 *  
## Residuals   54  712.1    13.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

可视化方法实例

1. 使用外部包HH

# 导入包
library(HH)

## Loading required package: lattice

## Loading required package: grid

## Loading required package: latticeExtra

## Loading required package: multcomp

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

## Loading required package: gridExtra

interaction2wt(len~supp*dose)

2. 直接使用内置函数

interaction.plot(dose, supp, len, type="b",
                 col=c("red","blue"), pch=c(16, 18),
                 main = "Interaction between Dose and Supplement Type")

3. 使用外部包gplots

library(gplots)

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:HH':
## 
##     residplot

## The following object is masked from 'package:stats':
## 
##     lowess

plotmeans(len ~ interaction(supp, dose, sep=" "),
          connect=list(c(1,3,5),c(2,4,6)),# 用于声明连接点
          col=c("red", "darkgreen"),
          main = "Interaction Plot with 95% CIs",
          xlab="Treatment and Dose Combination")

正交试验设计与方差分析

• 多因素问题的方差分析，因素过多导致试验次数大大增加，全面试验不现实；此时一般使用正交表安排试验。

正交表

例：记号$L_9(3^4)$：

• $L$为正交表；
• $9$为正交表的行数，表示需要试验的次数；
• $4$为正交表的列数，表示最多可以安排的因素的个数；
• $3$为因素水平数，表示此表可以安排三水平的试验。

一个例子：

# 输入数据框，有规律的因子可以使用`gl()`输入
rate <- data.frame(
    a=gl(3, 3),
    b=gl(3, 1, 9),
    c=factor(c(1,2,3,2,3,1,3,1,2)),
    y=c(31, 54, 38, 53, 49, 42, 57, 62, 64)
)
k = matrix(0, nrow=3, ncol=3, 
           dimnames=list(1:3, c('a', 'b', 'c')))

for (i in 1:3)
  for (j in 1:3)
    k[j, i]=mean(rate$y[rate[i]==j])

k

##    a  b  c
## 1 41 47 45
## 2 48 55 57
## 3 61 48 48

# 绘图
plot(as.vector(k))
lines(1:3,k[,1])
lines(4:6,k[,2])
lines(7:9,k[,3])

# 进行方差分析
rate.aov <- aov(y~a+b+c, data = rate)
summary(rate.aov)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## a            2    618     309  34.333 0.0283 *
## b            2    114      57   6.333 0.1364  
## c            2    234     117  13.000 0.0714 .
## Residuals    2     18       9                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

考虑样本间的交互作用

同样，采用正交表的方法，依次计算显著性，并且计算各因素之间的交互情况，最后通过 方差分析得到合理的因素选择。

有重复试验的方差分析

方法与正交表类似。