R语言_方差分析

R语言_学习笔记

Part 1. 前言介绍

1、常见研究设计的表达式

• 单因素ANOVA : y ~ A
• 含单个协变量的单因素ANCOVA : y ~ x + A
• 双因素ANOVA : y ~ A * B
• 含两个协变量的双因素ANCOVA : y ~ x1 + x2 + A * B
• 单因素组内ANOVA : y ~ A + Error(Subject/A)
• 含单个组内因子(W)和单个组间因子(B)的重复测量ANOVA : y ~ B * W + Error(Subject/W)

其中，小写字母表示定量变量，大写字母表示组别因子，Subject表示对被试者独有的表示变量（name)

2、表达式中各项的顺序

• 因子不止一个，且非平衡设计（不同组别因子下的观测样本数相等）
• 存在协变量

(1)、出现上述任意一种情况，等式右边的变量都与其他每个变量相关，表达式中的效应顺序就会造成影响。
(2)、当样本大小越不平衡，效应项的顺序对结果的影响越大。
(3)、一般来说，越基础性的效应越需要放在表达式前面。具体而言，首先是“协变量”，然后是“主效应”，接着是“双因素的交互项”，再接着是“三因素的交互项”。对于主效应，越基础性的变量越应放在表达式前面

Part 2. 单因素方差分析

1.1、单因素协方差分析 - R语言实现的方法

#一条一条Run，不要Source.
> library(mvtnorm)
> library(survival)
> library(splines)
> library(multcomp)
> library(gplots)
> attach(cholesterol)
> aggregate(response,by=list(trt),FUN=mean)   #用aggregate()函数的好处是输出的形式非常漂亮。
  Group.1        x
1   1time  5.78197
2  2times  9.22497
3  4times 12.37478
4   drugD 15.36117
5   drugE 20.94752
> fitav <- aov(response~trt)
> summary(fitav)
            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 
> plotmeans(response~trt,xlab="Treatment",ylab="Response",
            main="Mean Plot\nwith 95% CI",col="red")

1.2、单因素方差分析 - 假设检验

• 因变量应服从正态分布（正态性）

#用qqPlot()函数检验
> library(car)
> qqPlot(lm(response~trt,),
         simulate=TRUE,
         main="Q-Q Plot",
         labels=FALSE)

• 各组方差相等（同方差性）

> bartlett.test(response~trt)
    Bartlett test of homogeneity of variances
data:  response by trt
Bartlett's K-squared = 0.5797, df = 4, p-value = 0.9653

p > 0.05,符合同方差性。

2.1、单因素协方差分析 - R语言实现的方法

> library(colorspace)
> library(effects)
> library(multcomp)
> attach(litter)
> aggregate(weight,by=list(dose),FUN=mean)
  Group.1        x
1       0       32.30850
2       5       29.30842
3      50       29.86611
4     500       29.64647
> fitkl <- aov(weight~gesttime+dose)
> summary(fitkl)
            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  
#用effect()函数去除协变量效应的组均值。
> effect("dose",fitkl)
 dose effect
dose
       0        5       50      500 
32.35367 28.87672 30.56614 29.33460 
#在实际情况中，这个结果和直接求均值的结果会相差较远。

2.2 单因素协方差分析 - 假设检验

#正态性和同方差性检验，与ANOVA相同，此处不赘述。
#回归斜率的同质性检验
> fitkl2 <- aov(weight~gesttime*dose)
> summary(fitkl2)
              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   

3、单因素方差分析 - 成对组间比较

#此处使用glht()函数，因为它比TukeyHSD()函数更高大上！
> library(multcomp)
> par(mar=c(5,4,6,2))   #增大顶部边界面积，能更好地书写a、b、c、d
> tuk <- glht(fitav,linfct=mcp(trt="Tukey"))
> plot(cld(tuk,level=.05),col="light blue")

结论分析：

• 1time和2times的差异不显著（因为含有相同的字母a），1time和4times的差异显著（因为含有完全不同的字母）。
• 一天四次5mg剂量比一天两次10mg和一天一次20mg的效果更好。drugE比所有其它药物和疗法都要优。

4、稳健单因素方差分析

> kruskal.test(response~trt)
    Kruskal-Wallis rank sum test
data:  response by trt
Kruskal-Wallis chi-squared = 36.5421, df = 4, p-value = 2.238e-07

p < 0.05，表明不同组之间的中位数有显著差异。

Part 3. 双因素方差分析

1、R语言实现方法

#一条一条Run，不要Source.
> library(lattice)
> library(grid)
> library(latticeExtra)
> library(RColorBrewer)
> library(multcomp)
> library(mvtnorm)
> library(survival)
> library(splines)
> library(HH)
> attach(ToothGrowth)
> interaction2wt(len~supp*dose)

结论分析：

• 主效应：
  
  • “右上”图表明，不同剂量的抗坏血酸(dose)对牙齿长度的影响是不同的。
  • “左下”图表明，不同的喂食方式(supp)对牙齿长度的影响是不同的。
• 交互效应：
  
  • “右下”图表明，不同的喂食方式(supp)在不同剂量的抗坏血酸(dose)下，对牙齿长度的影响是不同的。

2、混合模型方差分析
又被称为“重复测量方差分析”。

> library(lattice)
> attach(CO2)
> wed <- subset(CO2,Treatment=='chilled')   #取出数据集CO2的子集
> attach(wed)
> fitavc <- aov(uptake~conc*Type+Error(Plant/conc),wed)   #在子集wed中做混合模型方差分析
> summary(fitavc)
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                   
---
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                     
---
Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30    869   28.97    

0.00148和0.000125表明主效应显著，0.001595表明交互效应显著。

#用interaction.plot()函数展示交互效应，分析交互效应是否显著。
> par(las=2,mar=c(10,4,4,2))
> interaction.plot(conc,Type,uptake,type="b",col=c("orange","royalblue3"),pch=c(16,18),
                   main="Interaction Plot for Plant Type and Concentration")

结论分析：

• 组间因子(Type)显著，说明uptake在Quebec和Mississippi上的吸收量显著不同。
• 组内因子(conc)显著，说明uptake在不同的CO2浓度下的吸收量显著不同。
• 交互效应显著，说明uptake在95~1000浓度下的吸收量在Quebec和Mississippi两个地区间是显著不同的。（两条线不平行，所以存在交互作用）

#若交互效应显著，用boxplot()函数展示交互效应的不同侧面。
boxplot(uptake~Type*conc,col=c("palegreen1","purple3"),
        main="Chilled Quebec andd Mississippi Plants",
        ylab="Carbon dioxide uptake rate (umol/m^2 sec)")

同一浓度下的箱线图，不处于同一水平线上，这就表明存在交互效应。

Part 4. 多元方差分析

1.1、R语言实现的方法

#一条一条Run，不要Source.
> library(MASS)
> attach(UScereal)
> y <- cbind(calories,fat,sugars)
> aggregate(y,by=list(shelf),FUN=mean)
> cov(y)    #生成各因变量方差与协方差。
  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
> fituk <- manova(y~shelf)
> summary(fituk)
          Df Pillai approx F num Df den Df    Pr(>F)    
shelf      2 0.4021   5.1167      6    122 0.0001015 ***
Residuals 62    

p=0.0001015 < 0.05，表明卡路里(calories)、脂肪(fat)、糖(sugars)均会因为储存架位置的不同而变化。

1.2、均值比较

#以sugars为例。
> library(multcomp)
> par(mar=c(5,4,6,2))   #增大顶部边界面积，能更好地书写a、b、c、d
> fitsk <- aov(sugars~shelf)
> tuk2 <- glht(fitsk,linfct=mcp(shelf="Tukey"))
> plot(cld(tuk2,level=.05),col="light blue")

1.3、假设检验

#检验多元正态性
> 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 > Assessing MUlutivariate Normality",
                  ylab="Mahalanobis D2")
> abline(a=0,b=1)
> identify(coord$x,coord$y,labels=row.names(UScereal))

2、稳健多元方差分析

#稳健多元方差分析,其对离群点和违反“多元正态性”和“方差——协方差矩阵同质性”的情况不敏感。
> library(robustbase)
> library(rrcov)
> Wilks.test(y,shelf,method="mcd")
    Robust One-way MANOVA (Bartlett Chi2)
data:  x
Wilks' Lambda = 0.5107, Chi2-Value = 23.883, DF = 4.868, p-value = 0.0002019
sample estimates:
  calories       fat    sugars
1 119.8210 0.7010828  5.663143
2 128.0407 1.1849576 12.537533
3 160.8604 1.6524559 10.352646

p=0.0002019 < 0.05，表明存储在货架顶部、中部和底部的谷物营养成分含量显著不同。