第四章多个样本均数比较的方差分析

x2yline

2017年10月18日

知识清单

完全随机设计的方差分析(Completely Randomized Design)
随机区组设计的方差分析(Randomized block design)
多个样本均值间的多重比较

课外

拉丁方设计资料的方差分析
两阶段交叉设计的方差分析
多样本方差比较的Bartlett检验和Levene检验

主要是aov函数的使用，可以直接对aov的结果直接plot

1. 完全随机设计的方差分析（completely random design）¹

*前提：来自正态分布的独立样本；各样本方差相等

测试数据：例04-02.sav

> # install.packages("memisc")
> library(memisc)
> group_df <- data.frame(as.data.set(spss.system.file('/mnt/e/医学统计学（第4版）/各章例题SPSS数据文件/例04-02.sav')))

> head(group_df)
    group ldl_c
1 placebo  3.53
2 placebo  4.59
3 placebo  4.34
4 placebo  2.66
5 placebo  3.59
6 placebo  3.13

> with(group_df, summary(aov(ldl_c ~group)))
             Df Sum Sq Mean Sq F value   Pr(>F)    
group         3  32.16  10.719   24.88 1.67e-12 ***
Residuals   116  49.97   0.431                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

2. 随机区组设计的方差分析（randomized block design）

测试数据：例04-04.sav

> group_df <- data.frame(as.data.set(spss.system.file('E:/医学统计学（第4版）/各章例题SPSS数据文件/例04-04.sav')))
> head(group_df)
  group treat weight
1     1   A药   0.82
2     2   A药   0.73
3     3   A药   0.43
4     4   A药   0.41
5     5   A药   0.68
6     1   B药   0.65

> with(group_df, summary(aov(weight ~ factor(treat)+factor(group))))
              Df Sum Sq Mean Sq F value  Pr(>F)   
factor(treat)  2 0.2280 0.11400  11.937 0.00397 **
factor(group)  4 0.2284 0.05709   5.978 0.01579 * 
Residuals      8 0.0764 0.00955                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

3. 多个样本均值间的多重比较

若用两样本均数比较的t检验进行多重比较，将会加大犯Ⅰ类错误（把本无差别的两个总体均数判为有差别）的概率。

当方差分析的结果为拒绝H0，接受H1时，说明g个总体均数不全相等。

若想进一步了解哪些两个总体均数不等，需进行多个样本均数间的两两比较或称多重比较。

3.1 LSD-t检验 (least significant difference, LSD)

查表为t分布表，使用条件：一对或几对比较

检验统计量：

\[ LSD-t = \frac{\bar{X_{i}}-\bar{X_{j}}}{S_{\bar{X_{i}}-\bar{X_{j}}}}, \ \nu=\nu_{误差}\\ S_{\bar{X_{i}}-\bar{X_{j}}} = \sqrt{MS_{误差}(\frac{1}{n_{i}}+\frac{1}{n_{j}})} \]

LSD-t检验公式与两样本均数比较的t检验公式区别在于检验统计量和自由度ν的计算上。

示例数据：例04-02.sav

方法1：

> group_df <- data.frame(as.data.set(spss.system.file('E:/医学统计学（第4版）/各章例题SPSS数据文件/例04-02.sav')))
> r_mean_sq <- summary(aov(group_df$ldl_c ~ factor(group_df$group)))[[1]][-1, 3]
> v_lsd <- summary(aov(group_df$ldl_c ~ factor(group_df$group)))[[1]][-1, 1]
> ns <- table(group_df$group)
> S_lsd <- sqrt(r_mean_sq*(1/ns[c("2.4g")]+1/ns[c("placebo")]))
> delta_x <- mean(group_df$ldl_c[group_df$group=="placebo"])-mean(group_df$ldl_c[group_df$group=="2.4g"])
> pt(abs(delta_x/S_lsd), df=v_lsd, lower.tail=F)*2
        2.4g 
4.887216e-05

方法2：

> library(agricolae)
> lsd_result <- LSD.test(aov(group_df$ldl_c ~ (group_df$group)), "group_df$group", group=T, alpha=0.01)
> lsd_result # 分组为相同说明无显著差异，如2.4g与4.8g
$statistics
    MSerror  Df   Mean      CV  t.value       LSD
  0.4307502 116 2.7025 24.2855 2.618878 0.4437949

$parameters
        test p.ajusted         name.t ntr alpha
  Fisher-LSD      none group_df$group   4  0.01

$means
        group_df$ldl_c       std  r      LCL      UCL  Min  Max    Q25   Q50    Q75
2.4g          2.715333 0.6381586 30 2.401523 3.029144 1.56 4.32 2.3175 2.665 2.9650
4.8g          2.698000 0.4971671 30 2.384190 3.011810 1.68 3.68 2.3375 2.655 2.9800
7.2g          1.966333 0.7464421 30 1.652523 2.280144 0.89 3.71 1.3350 1.905 2.4225
placebo       3.430333 0.7151247 30 3.116523 3.744144 1.37 4.59 2.9650 3.530 3.9825

$comparison
NULL

$groups
        group_df$ldl_c groups
placebo       3.430333      a
2.4g          2.715333      b
4.8g          2.698000      b
7.2g          1.966333      c

attr(,"class")
[1] "group"

方法3：

> pairwise.t.test(group_df$ldl_c, group_df$group)

    Pairwise comparisons using t tests with pooled SD 

data:  group_df$ldl_c and group_df$group 

     placebo 2.4g    4.8g   
2.4g 0.00013 -       -      
4.8g 0.00013 0.91871 -      
7.2g 2.1e-13 0.00011 0.00013

P value adjustment method: holm

参考：
https://stackoverflow.com/questions/11454521/r-t-test-and-pairwise-t-test-give-different-results

3.2 Dunnett-t检验

适用条件：g-1个实验组与一个对照组均数差别的多重比较，检验统计量Dunnett-t 有专门的界值表，是LSD-t的特殊化

> library(multcomp)
载入需要的程辑包：mvtnorm
载入需要的程辑包：survival
载入需要的程辑包：TH.data

载入程辑包：‘TH.data’

The following object is masked from ‘package:MASS’:

    geyser

> cht <- glht(aov(ldl_c ~ group, data=group_df), linfct = mcp(group= "Tukey"))
> summary(cht, test = univariate())

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = ldl_c ~ group, data = group_df)

Linear Hypotheses:
                    Estimate Std. Error t value Pr(>|t|)    
2.4g - placebo == 0 -0.71500    0.16946  -4.219 4.89e-05 ***
4.8g - placebo == 0 -0.73233    0.16946  -4.322 3.29e-05 ***
7.2g - placebo == 0 -1.46400    0.16946  -8.639 3.57e-14 ***
4.8g - 2.4g == 0    -0.01733    0.16946  -0.102    0.919    
7.2g - 2.4g == 0    -0.74900    0.16946  -4.420 2.23e-05 ***
7.2g - 4.8g == 0    -0.73167    0.16946  -4.318 3.34e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Univariate p values reported)

> cht <- glht(aov(ldl_c ~ group, data=group_df), linfct = mcp(group= "Dunnett"), alternative="two.sided")
> summary(cht)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = ldl_c ~ group, data = group_df)

Linear Hypotheses:
                    Estimate Std. Error t value Pr(>|t|)    
2.4g - placebo == 0  -0.7150     0.1695  -4.219 0.000131 ***
4.8g - placebo == 0  -0.7323     0.1695  -4.322  < 1e-04 ***
7.2g - placebo == 0  -1.4640     0.1695  -8.639  < 1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

3.3 SNK-q检验

适用条件：多个样本均数两两之间的全面比较（比较所有的组）

> # 新用法
> agricolae::SNK.test(aov(ldl_c ~ group, data=group_df), "group",  console = T)

Study: aov(ldl_c ~ group, data = group_df) ~ "group"

Student Newman Keuls Test
for ldl_c 

Mean Square Error:  0.4307502 

group,  means

           ldl_c       std  r  Min  Max
2.4g    2.715333 0.6381586 30 1.56 4.32
4.8g    2.698000 0.4971671 30 1.68 3.68
7.2g    1.966333 0.7464421 30 0.89 3.71
placebo 3.430333 0.7151247 30 1.37 4.59

Alpha: 0.05 ; DF Error: 116 

Critical Range
        2         3         4 
0.3356368 0.4023275 0.4417253 

Means with the same letter are not significantly different.

           ldl_c groups
placebo 3.430333      a
2.4g    2.715333      b
4.8g    2.698000      b
7.2g    1.966333      c
> print(agricolae::SNK.test(aov(ldl_c ~ group, data=group_df), "group", alpha=0.01))
$statistics
    MSerror  Df   Mean      CV
  0.4307502 116 2.7025 24.2855

$parameters
  test name.t ntr alpha
   SNK  group   4  0.01

$snk
     Table CriticalRange
2 3.703652     0.4437949
3 4.202736     0.5035983
4 4.500339     0.5392590

$means
           ldl_c       std  r  Min  Max    Q25   Q50    Q75
2.4g    2.715333 0.6381586 30 1.56 4.32 2.3175 2.665 2.9650
4.8g    2.698000 0.4971671 30 1.68 3.68 2.3375 2.655 2.9800
7.2g    1.966333 0.7464421 30 0.89 3.71 1.3350 1.905 2.4225
placebo 3.430333 0.7151247 30 1.37 4.59 2.9650 3.530 3.9825

$comparison
NULL

$groups
           ldl_c groups
placebo 3.430333      a
2.4g    2.715333      b
4.8g    2.698000      b
7.2g    1.966333      c

attr(,"class")
[1] "group"

# 旧用法
SNK.test(y, treat, DFerror, MSerror)

4. 拉丁方设计资料的方差分析

拉丁方的前提条件是横行、直列单位组因素与试验因素间不存在交互作用²

在临床试验中，拉丁方设计的行区组通常是受试对象，列区组通常是重复观测次数，并要求前一个处理的作用不会影响到下一个处理³

拉丁方实验设计与区组实验设计一样，都是为了平衡额外变量，以防止这些额外变量成为混淆因子，破坏实验研究的内部效度。如果简化点来解释，一般来说，区组实验设计多用于对一个额外变量的平衡，如被试因素、时间顺序因素、空间位置因素等；拉丁方实验设计则可以看成是区组设计的扩展，即扩展到可以平衡两个额外变量（当然，如果设计巧妙，也可以达到对多于两个额外变量的平衡，但那也是在二维平衡模式上变化出来的）⁴⁵

拉丁方实验设计既有优点也有缺点。其优点是，在许多研究情境中，这种设计比完全随机和随机区组设计更加有效，它可以使研究者平衡并分离出两个额外变量的影响，因而减小实验误差，可获得对实验处理效应的更精确的估价。另外，通过对方格单元内误差与残差的F检验，可以检验额外变量与自变量是否有交互作用，以检验采用拉丁方设计是否合适。

拉丁方设计的缺点是，它的关于自变量与额外变量不存在交互作用的假设在很多情况下都难以保证，尤其当实验中含有多个自变量的时候。因此，拉丁方实验设计在多因素实验中不常用。另外，拉丁方实验设计要求每个额外变量的水平数与实验处理数必须相等，这也在一定程度上限制了拉丁方实验设计的使用。

数据准备

> library("memisc")
> group_df <- data.frame(as.data.set(spss.system.file('E:/医学统计学（第4版）/各章例题SPSS数据文件/例04-05_1.sav')))
> head(group_df)
  row column treat result
1   1      1     C     87
2   1      2     B     75
3   1      3     E     81
4   1      4     D     75
5   1      5     A     84
6   1      6     F     66

方差分析

> group_df$row <- factor(group_df$row)
> group_df$col <- factor(group_df$column)
> model<- aov(result~row + col + treat, data=group_df)
> summary(model)
            Df Sum Sq Mean Sq F value Pr(>F)  
row          5  250.5   50.09   1.424 0.2584  
col          5   75.1   15.03   0.427 0.8242  
treat        5  657.3  131.47   3.738 0.0149 *
Residuals   20  703.4   35.17                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5. 两阶段交叉设计的方差分析

即按个体分解，按阶段分解和按处理分解进行方差分析，代码与拉丁方的方法大致相同

6. Bartlett检验（需要正态整体）⁶

> library("memisc")
> group_df <- data.frame(as.data.set(spss.system.file('E:/医学统计学（第4版）/各章例题SPSS数据文件/例04-02.sav')))
> head(group_df)
    group ldl_c
1 placebo  3.53
2 placebo  4.59
3 placebo  4.34
4 placebo  2.66
5 placebo  3.59
6 placebo  3.13
> bartlett.test(ldl_c~group, data=group_df) # bartlett.test(group_df$ldl_c,group_df$group)

    Bartlett test of homogeneity of variances

data:  ldl_c by group
Bartlett's K-squared = 5.2192, df = 3, p-value = 0.1564

7. Levene检验（不必要正态整体）【Fligner-Killeen test也很常用】⁷⁸

> library(car)

载入程辑包：‘car’

The following object is masked from ‘package:memisc’:

    recode

> leveneTest(ldl_c~group, data=group_df)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   3   1.493 0.2201
      116

Reference:

Completely randomized design, http://www.r-tutor.com/elementary-statistics/analysis-variance/completely-randomized-design ↩
拉丁方设计之R篇, http://blog.sciencenet.cn/blog-1114360-745402.html ↩
拉丁方设计, https://wenku.baidu.com/view/c8cb8ad9a58da0116c1749e2.html ↩
如何理解“拉丁方实验设计”(《实验心理学》问答009), http://blog.sina.com.cn/s/blog_49b5655001000d66.html ↩
心理与教育研究中的多因素实验设计 (北京师范大学出版社) (1994)↩
Bartlett test of homogeneity of variances, https://stat.ethz.ch/R-manual/R-devel/library/stats/html/bartlett.test.html ↩
Homogeneity of variance, http://www.cookbook-r.com/Statistical_analysis/Homogeneity_of_variance/↩
R语言: 用于多个比较的方差分析（ANOVA）, http://blog.sina.com.cn/s/blog_5cd2f1e20101979p.html ↩