By

文章目录
  1. 方差分析
    1. 1、术语
    2. 2、ANOVA模型拟合
      1. 2.1、aov()函数
      2. 2.2、顺序很重要
    3. 3、单因素方差分析
      1. 3.1、多重比较
      2. 3.2、检验的假设条件
    4. 4、协方差分析
      1. 4.1、检验的假设条件
      2. 4.2、可视化
    5. 5、双因素方差分析
    6. 6、重复测量方差分析
    7. 7、多元方差分析
      1. 7.1、评估假设检验
      2. 7.2、稳健多元方差分析
    8. 8、用回归做ANOVA

方差分析

  • 当包含的因子是解释变量时，我们关注的重点通常会从预测转向组别差异的分析，这种分析法称作方差分析 （ANOVA）。

    1、术语

  • 均衡/非均衡设计：不同观测方案的观测数相同/不同
  • 单因素方差分析：仅有一个类别型变量，一般为组间两水平因子（两种药物）
  • 单因素组内方差分析：仅有一个类型形变量，一般为组内两水平因子（同种药物，作用时间不同）
  • 因素方差分析：设计包含两个甚至更多的因子
  • 混合模型方差分析：因子设计包括组内和组间因子
  • 混淆因素/干扰因子：你不关心的因素，但其能解释你关心的组间差异
  • 协方差分析（ANCOVA）：度量并将混淆因子作为协变量考虑在模型中的方差分析
  • 多元方差分析（MANOVA）：因变量不止一个
  • 多元协方差分析（MANCOVA）：因变量不止一个的协方差分析

2、ANOVA模型拟合

  • lm函数与aov函数都可以分析方差分析模型，ANOVA和回归本质上都属于广义线性模型

    2.1、aov()函数

1
2
aov(formula, data = dataframe)
#formula可使用的参数符号见回归（一）

一些常见的研究设计表达式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#单因素ANOVA
y ~ A
#含单个协变量的单因素ANCOVA 
y ~ x + A
#双因素ANOVA 
y ~ A * B
#含两个协变量的双因素ANCOVA 
y ~ x1 + x2 + A*B
#随机化区组 
y ~ B + A（B是区组因子）
#单因素组内ANOVA  Subject标识一个独特的被试者
y ~ A + Error(Subject/A)
#含单个组内因子（W）和单个组间因子（B）的重复测量ANOVA 
y ~ B * W + Error(Subject/W)

2.2、顺序很重要

  • 例如，对于双因素方差分析，不同的处理方式观测数不同，则y~AB与y~BA结果不同

假设你正使用如下表达式对数据进行建模：

  • Y ~ A + B + A:B

有三种类型的方法可以分解等式右边各效应对y所解释的方差。

  • 类型I（序贯型）：效应根据表达式中先出现的效应做调整。 A不做调整， B根据A调整， A:B交互项根据A和B调整。
  • 类型II（分层型）：效应根据同水平或低水平的效应做调整。 A根据B调整， B依据A调整，A:B交互项同时根据A和B调整。
  • 类型III（边界型）：每个效应根据模型其他各效应做相应调整。 A根据B和A:B做调整， A:B交互项根据A和B调整。

样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说，越基础性的效应越需要放在表达式前面。具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。

  • R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型III方法。
  • car包中的Anova()函数（不要与标准anova()函数混淆）提供了使用类型II或类型III方法的选项。

3、单因素方差分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
> #单因素方差分析
> library(multcomp)
> attach(cholesterol)
The following objects are masked from cholesterol (pos = 4):
    response, trt
> table(trt)#各组样本大小
trt
 1time 2times 4times  drugD  drugE 
    10     10     10     10     10 
> aggregate(response,by=list(trt),FUN=mean)#计算response的均值
  Group.1        x
1   1time  5.78197
2  2times  9.22497
3  4times 12.37478
4   drugD 15.36117
5   drugE 20.94752
> aggregate(response,by=list(trt),FUN=sd)#计算response的标准差
  Group.1        x
1   1time 2.878113
2  2times 3.483054
3  4times 2.923119
4   drugD 3.454636
5   drugE 3.345003
> 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
> #ANOVA对治疗方式（trt）的F检验非常显著（p<0.0001），说明五种疗法的效果不同
> library(gplots)
> 
> #图形展示了带有95%的置信区间的各疗法均值，可以清楚看到它们之间的差异
> plotmeans(response~trt, xlab="Treatment", ylab="Response", main="Mean Plot")
> detach(cholesterol)

3.1、多重比较

1
2
3
4
5
TukeyHSD(fit)
par(las=2)
par(mar=c(5,8,4,2))
plot(TukeyHSD(fit))
#图形中置信区间包含0的疗法说明差异不显著（p>0.05）。

1
2
3
4
5
library(multcomp)
par(mar=c(5,4,6,2))
tuk<-glht(fit,linfct=mcp(trt="Tukey"))
plot(cld(tuk, level = .05),col="lightgrey")
#有相同字母的组差异不显著，图中，1time和2times差异不显著，1time和4times差异显著

3.2、检验的假设条件

  • 样本随机独立
  • 因变量的正态性
1
2
3
#单因素方差分析中，我们假设因变量服从正态分布，各组方差相等。可以使用Q-Q图来检验正态性假设
library(car)
qqPlot(lm(response~trt,data=cholesterol),simulate=TRUE,main="Q-Q Plot", labels=FALSE)

  • 样本的方差齐性（总体方差相等）

注意：F检验对于数据的正态性非常敏感，因此在检验方差齐性的时候，Levene检验, Bartlett检验或者Brown–Forsythe检验的稳健性都要优于F检验。 F检验还可以用于三组或者多组之间的均值比较，但是如果被检验的数据无法满足均是正态分布的条件时，该数据的稳健型会大打折扣，特别是当显著性水平比较低时。但是，如果数据符合正态分布，而且alpha值至少为0.05，该检验的稳健型还是相当可靠的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#Bartlett检验
> bartlett.test(response~trt,data=cholesterol)
	Bartlett test of homogeneity of variances
data:  response by trt
Bartlett K-squared = 0.57975, df = 4,
p-value = 0.9653
#Bartlett检验表明五组的方差并没有显著不同（p=0.97）。
#方差齐性分析对离群点非常敏感
> fit<-aov(response~trt)
> outlierTest(fit) #识别离群
No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
   rstudent unadjusted p-value Bonferonni p
19 2.251149           0.029422           NA
#从输出结果来看，并没有证据说明胆固醇数据中含有离群点（当p>1时将产生NA） 。

4、协方差分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
> #协方差分析
> #怀孕小鼠分组接收不同剂量药物，因变量幼崽体重，协变量怀孕时间
> data(litter,package = "multcomp")
> attach(litter)
The following objects are masked from litter (pos = 4):
    dose, gesttime, number, weight
> table(dose)
dose
  0   5  50 500 
 20  19  18  17 
> aggregate(weight,by=list(dose),FUN=mean)#计算weight的均值
  Group.1        x
1       0 32.30850
2       5 29.30842
3      50 29.86611
4     500 29.64647
> 
> fit<-aov(weight~gesttime+dose)
> summary(fit)
            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
> #我们得到：(a)怀孕时间与幼崽出生体重相关； (b)控制怀孕时间，药物剂量与出生体重相关。
> 
> #去除协变量效应后的组均值
> library(effects)
> effect("dose",fit)
 dose effect
dose
       0        5       50      500 
32.35367 28.87672 30.56614 29.33460 
> 
> #multcomp包校验用户自定义的均值假设
> library(multcomp)
> #c(3, -1, -1, -1)设定第一组和其他三组的均值进行比较
> contrast<-rbind("no drug vs. drug"=c(3,-1,-1,-1))
> summary(glht(fit,linfct=mcp(dose=contrast)))
	 Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: User-defined Contrasts
Fit: aov(formula = weight ~ gesttime + dose)
Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)  
no drug vs. drug == 0    8.284      3.209   2.581    0.012 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
> #可以得出未用药组比其他用药条件下的出生体重高的结论

4.1、检验的假设条件

1
2
3
4
5
6
7
8
9
10
11
12
13
> #假设条件
> #正态性和同方差性同方差分析
> #回归斜率相同
> fit2<-aov(weight~gesttime*dose,data = litter)
> summary(fit2)
              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
> #可以看到交互效应不显著，支持了斜率相等的假设。

4.2、可视化

1
2
3
4
5
#可视化
library(HH)
ancova(weight~gesttime+dose,data=litter)
#可以看到，用怀孕时间来预测出生体重的回归线相互平行，只是截距项不同。随着怀
#孕时间增加，幼崽出生体重也会增加。另外，还可以看到0剂量组截距项最大， 5剂量组截距项最小。

5、双因素方差分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
> #双因素方差分析
> attach(ToothGrowth)
The following objects are masked from ToothGrowth (pos = 3):
    dose, len, supp
> 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
> aggregate(len,by=list(supp,dose),FUN=sd)
  Group.1 Group.2        x
1      OJ     0.5 4.459709
2      VC     0.5 2.746634
3      OJ     1.0 3.910953
4      VC     1.0 2.515309
5      OJ     2.0 2.655058
6      VC     2.0 4.797731
> fit<-aov(len~supp*dose)
> summary(fit)
            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
> #可以看到主效应（supp和dose）和交互效应都非常显著
#可视化 推荐interaction2wt可展示双因素、三因素等
interaction2wt(len~supp*dose)

6、重复测量方差分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
> #重复测量方差分析
> #数据为不同浓度二氧化碳环境下两种植物对其的吸收量
> w1b1<-subset(CO2,Treatment=='chilled')
> fit<-aov(uptake~conc*Type+Error(Plant/conc),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       1  888.6   888.6  215.46 0.000125 ***
conc:Type  1  239.2   239.2   58.01 0.001595 ** 
Residuals  4   16.5     4.1                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30    869   28.97               
> 
> par(las=2)
> par(mar=c(10,4,4,2))
> with(w1b1,interaction.plot(conc,Type,uptake,type="b",col=c("red","blue"),pch = c(16,18),main=
+                              "植物类型和二氧化碳浓度的交互图"))
> boxplot(uptake~conc*Type,data=w1b1,col=(c("gold","green")),main="魁北克与密西西比植物对二氧化碳吸收的交互图")
> #魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而且随着CO2浓度的升高，差异越来越明显。
> 
> #注意 该方法为传统的重复测量方差分析。该方法假设任意组内因子的协方差矩阵为球形，并且任意组内因子两水平间的方差之差都相等。


7、多元方差分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
> #研究美国谷物中的卡路
> #里、脂肪和糖含量是否会因为储存架位置的不同而发生变化。其中1代表底层货架， 2代表中层货
> #架， 3代表顶层货架。卡路里、脂肪和糖含量是因变量，货架是三水平（1、 2、 3）的自变量。
> library(MASS)
> attach(UScereal)
The following objects are masked from UScereal (pos = 7):
    calories, carbo, fat, fibre, mfr, potassium, protein, shelf, sodium, sugars, vitamins
The following objects are masked from UScereal (pos = 9):
    calories, carbo, fat, fibre, mfr, potassium, protein, shelf, sodium, sugars, vitamins
> y<-cbind(calories,fat,sugars)
> aggregate(y,by=list(shelf),FUN=mean)
  Group.1 calories       fat    sugars
1       1 119.4774 0.6621338  6.295493
2       2 129.8162 1.3413488 12.507670
3       3 180.1466 1.9449071 10.856821
> fit<-manova(y~shelf)
> summary(fit)
          Df  Pillai approx F num Df den Df  Pr(>F)   
shelf      1 0.19594    4.955      3     61 0.00383 **
Residuals 63                                          
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> summary.aov(fit)#输出单变量结果
 Response calories :
            Df Sum Sq Mean Sq F value    Pr(>F)    
shelf        1  45313   45313  13.995 0.0003983 ***
Residuals   63 203982    3238                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 Response fat :
            Df  Sum Sq Mean Sq F value   Pr(>F)   
shelf        1  18.421 18.4214   7.476 0.008108 **
Residuals   63 155.236  2.4641                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 Response sugars :
            Df  Sum Sq Mean Sq F value  Pr(>F)  
shelf        1  183.34  183.34   5.787 0.01909 *
Residuals   63 1995.87   31.68                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

7.1、评估假设检验

  • 多元正太性
1
2
3
4
5
6
7
8
center<-colMeans(y)
n<-nrow(y)
p<-ncol(y)
cov<-cov(y)
d<-mahalanobis(y,center,cov)
coord<-qqplot(qchisq(ppoints(n),df=p),d,main="Q-Q Plot",ylab="Maalanobis D2")
abline(a=0,b=1)
identify(coord$x,coord$y,labels=row.names(UScereal))

  • 方差-协方差矩阵同质性

    暂时没找到合适的方法来进行检验

  • 多元离群点

    1
    2
    library(mvoutlier)
    aq.plot(y)

7.2、稳健多元方差分析

  • 多元正态性或方差-协方差同质性不满足或担心离群点时
  • rrcov包中Wilks.test(y.dataframe,x,method=”mcp”)
  • vegan包中的adonis()
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    #稳健多元方差分析
    library(MASS)
    attach(UScereal)
    #稳健多元单因素方差分析
    library(rrcov)
    y<-cbind(calories,fat,sugars)
    Wilks.test(y,shelf,method="mcd")#性能较慢
    #非参数多元方差分析
    library(vegan)
    adonis(y~shelf)

8、用回归做ANOVA

  • 当lm()函数碰到因子时，它会用一系列与因子水平相对应的数值型对照变量来代替因子。（s1,s2,s3=>(0,0,0),(0,1,0),(0,0,1)）
创建方法 描述
contr.helmert 第二个水平对照第一个水平，第三个水平对照前两个的均值，第四个水平对照前三个的均值，以此类推
contr.poly 基于正交多项式的对照，用于趋势分析（线性、二次、三次等）和等距水平的有序因子
contr.sum 对照变量之和限制为0。也称作偏差找对，对各水平的均值与所有水平的均值进行比较
contr.treatment 各水平对照基线水平（默认第一个水平）。也称作虚拟编码
contr.SAS 类似于contr.treatment，只是基线水平变成了最后一个水平。生成的系数类似于大部分SAS过程中使用的对照变量
  • fit.lm <- lm(response ~ trt, data=cholesterol, contrasts=”contr.SAS”)的contrasts参数来改变创建方法，默认的创建方法为contr.treatment) 
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    > ##使用回归进行方差分析
    > library(multcomp)
    > levels(cholesterol$trt)
    [1] "1time"  "2times" "4times" "drugD"  "drugE" 
    > fit.aov <- aov(response ~ trt, data=cholesterol)
    > summary(fit.aov)
                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
    > 
    > fit.lm <- lm(response ~ trt, data=cholesterol)
    > summary(fit.lm) #二者一致
    Call:
    lm(formula = response ~ trt, data = cholesterol)
    Residuals:
        Min      1Q  Median      3Q     Max 
    -6.5418 -1.9672 -0.0016  1.8901  6.6008 
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)    5.782      1.021   5.665 9.78e-07 ***
    trt2times      3.443      1.443   2.385   0.0213 *  
    trt4times      6.593      1.443   4.568 3.82e-05 ***
    trtdrugD       9.579      1.443   6.637 3.53e-08 ***
    trtdrugE      15.166      1.443  10.507 1.08e-13 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    Residual standard error: 3.227 on 45 degrees of freedom
    Multiple R-squared:  0.7425,	Adjusted R-squared:  0.7196 
    F-statistic: 32.43 on 4 and 45 DF,  p-value: 9.819e-13
文章目录
  1. 方差分析
    1. 1、术语
    2. 2、ANOVA模型拟合
      1. 2.1、aov()函数
      2. 2.2、顺序很重要
    3. 3、单因素方差分析
      1. 3.1、多重比较
      2. 3.2、检验的假设条件
    4. 4、协方差分析
      1. 4.1、检验的假设条件
      2. 4.2、可视化
    5. 5、双因素方差分析
    6. 6、重复测量方差分析
    7. 7、多元方差分析
      1. 7.1、评估假设检验
      2. 7.2、稳健多元方差分析
    8. 8、用回归做ANOVA