R action 9

20180502（从有道迁移）

方差分析

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

1. 术语

   1. 在每种方案下观测数相等，这种设计也称为均衡设计（balanced design）；若观测数不同，则称作非均衡设计（unbalanced design）
   2. 仅有一个类别型变量
      • 称为单因素方差分析（one-way ANOVA），或进一步称为单因素组间方差分析；
      • 在这种情况下，方差分析主要通过F检验来进行效果评测，若治疗方案的F检验显著，则说明方案的因变量均值不同，不同方案之间存在差异；
   3. 当同一方案在不同时间内进行比较时，时间（time）是两水平（统一方案下不同时间）的组内因子。
      • 因为对所有评估对象在所有水平下都进行了测量，所以这种统计设计称单因素组内方差分析，又由于每个评估对象都不止一次被测量，也称作重复测量方差分析。
      • 当时间的F检验显著时，说明测量对象的因变量在不同时间的测量下发生了改变
   4. 当设计包含两个甚至更多的因子时，便是因素方差分析设计，比如两因子时称作双因素方差分析，三因子时称作三因素方差分析，以此类推。若因子设计包括组内和组间因子，又称作混合模型方差分析
   5. 示例：随机分配五位患者到CBT，另外五位到EMDR，在五周和六个月后分别评价他们的STAI结果
      • 含组间和组内因子的双因素方差分析
      • 将做三次F检验：
        • 疗法因素一次，若疗法结果显著，说明CBT和EMDR对焦虑症的治疗效果不同
        • 时间因素一次，若时间结果显著，说明焦虑度从五周到六个月发生了变化
        • 两者交互因素一次，若两者交互效应显著，说明两种疗法随着时间变化对焦虑症治疗影响不同（也就是说，焦虑度从五周到六个月的改变程度在两种疗法间是不同的）
   6. 混淆因素（confounding factor），当对这个因素不感兴趣时，它也被称作干扰变数（nuisance variable）。
   7. 当存在协变量，该设计为协方差分析（ANCOVA）。
   8. 当因变量不止一个时，设计被称作多元方差分析（MANOVA），若协变量也存在，那么就叫多元协方差分析（MANCOVA）

2. ANOVA 模型拟合

   1. aov()函数，语法为aov(formula, data=dataframe)；

      • R表达式中的特殊符号

        符 号	用 法
        ~	分隔符号，左边为响应变量，右边为解释变量。</br>例如，用 A、B 和 C 预测 y，代码为 y ~ A + B + C
        :	表示变量的交互项。</br>例如，用 A、B 和 A 与 B 的交互项来预测 y，代码为 y ~ A + B + A:B
        *	表示所有可能交互项。</br>代码 y ~ A * B * C 可展开为 y ~ A + B + C + A:B + A:C + B:C + A:B:C
        ^	表示交互项达到某个次数。</br>代码 y ~ (A + B + C)^2 可展开为 y ~ A + B + C + A:B + A:C + B:C
        .	表示包含除因变量外的所有变量。</br>例如，若一个数据框包含变量 y、A、B 和 C，代码 y ~ .可展开为 y ~ A + B + C

      • 常见研究设计的表达式

        设 计	表 达 式
        单因素	ANOVA y ~ A
        含单个协变量的单因素	ANCOVA y ~ x + A
        双因素	ANOVA y ~ A * B
        含两个协变量的双因素	ANCOVA y ~ x1 + x2 + A*B
        随机化区组	y ~ B + A（B 是区组因子）
        单因素组内	ANOVA y ~ A + Error(Subject/A)
        含单个组内因子（W）和单个组间因子（B）的重复测量	ANOVA y ~ B * W + Error(Subject/W)

   2. 表达式中各项的顺序

      表达式中效应的顺序在两种情况下会造成影响：(a)因子不止一个，并且是非平衡设计；(b)存在协变量。出现任意一种情况时，等式右边的变量都与其他每个变量相关。

      image

      car包中的Anova()函数（不要与标准anova()函数混淆）提供了使用类型Ⅱ或类型Ⅲ方法的选项，而aov()函数使用的是类型I方法

3. 单因素方差分析

   1. 单因素方差分析，是比较分类因子定义的两个或多个组别中的因变量均值，可以使用：
      • 各组均值：aggregate(response, by=list(trt), FUN=mean)
      • 各组标准差： aggregate(response, by=list(trt), FUN=sd)
      • 检验组间差异（ANOVA）： fit <- aov(response ~ trt)
      • gplots包中的plotmeans()可以用来绘制带有置信区间的组均值图形
   2. 多重比较，ANOVA对各方案的F检验表明方案的效果不同，但是并没有说明哪种方案与其他疗法不同

      • TukeyHSD()函数提供了对各组均值差异的成对检验

      • multcomp包中的glht()函数提供了多重均值比较更为全面的方法，既适用于线性模型，也适用于广义线性模型。

        library(multcomp)
        par(mar=c(5,4,6,2))
        tuk <- glht(fit, linfct=mcp(trt="Tukey"))
        plot(cld(tuk, level=.05),col="lightgrey")
        par(opar)
        

   3. 评估检验的假设条件

      • 单因素方差分析中，我们假设因变量服从正态分布，各组方差相等。可以使用Q-Q图来检验正态性假设

        library(car)
        qqPlot(lm(response ~ trt, data=cholesterol), 
               simulate=TRUE, main="Q-Q Plot", labels=FALSE)
        

      • R提供了一些可用来做方差齐性检验的函数。例如，可以通过做Bartlett检验, Fligner-Killeen检验（fligner.test()函数）和Brown-Forsythe检验（HH包中的hov()函数

        bartlett.test(response ~ trt, data=cholesterol)
        
                Bartlett test of homogeneity of variances
        
        data:  cholesterol$response by cholesterol$trt
        Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653
        
        # Bartlett检验表明五组的方差并没有显著不同（p=0.97）
        

      • 方差齐性分析对离群点非常敏感。可利用car包中的outlierTest()函数来检测离群点

        library(car)
        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）
        
        

      • 因此根据Q-Q图、Bartlett检验和离群点检验，该数据似乎可以用ANOVA模型拟合得很好

4. 单因素协方差分析

   单因素协方差分析（ANCOVA）扩展了单因素方差分析（ANOVA），包含一个或多个定量的协变量

   1. 分析：

      • 单因素协方差分析可以先对自变量进行table和mean的分析，接着对自变量和协变量使用aov进行F校验；
      • 使用了协变量，想要获取调整的组均值，即去除协变量效应后的组均值。可使用effects包中的effects()函数来计算调整的均值
      • 对用户定义的对照的多重比较

      # Listing 9.3 - One-way ANCOVA
      data(litter, package="multcomp")
      attach(litter)
      table(dose) 
      aggregate(weight, by=list(dose), FUN=mean)
      fit <- aov(weight ~ gesttime + dose)                             
      summary(fit)
      
      # Obtaining adjusted means
      library(effects)
      effect("dose", fit)
      
      
      #  Listing 9.4 - Multiple comparisons using user supplied contrasts
      library(multcomp)
      contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1))
      summary(glht(fit, linfct=mcp(dose=contrast)))
      
      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.0120 *
      ---
      Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
      
      # 对照c(3, -1, -1, -1)设定第一组和其他三组的均值进行比较。假设检验的t统计量（2.581）在p<0.05水平下显著，因此，可以得出未用药组比其他用药条件下的出生体重高的结论
      
      

   2. 评估检验的假设条件

      • ANCOVA与ANOVA相同，都需要正态性和同方差性假设，另外，ANCOVA还假定回归斜率相同

   3. 结果可视化

      • HH包中的ancova()函数可以绘制因变量、协变量和因子之间的关系图

        library(HH) 
        ancova(weight ~ gesttime + dose, data=litter) 
        

5. 双因素方差分析

   1. 需要关注样本是否是均衡设计的，如果不是均衡设计，则需要考虑因素间顺序的问题；
   2. 对双因素方差分析有几种可视化方法：
      • interaction.plot()函数来展示双因素方差分析的交互效应；
      • gplots包中的plotmeans()函数来展示交互效应；
      • HH包中的interaction2wt()。
   3. 模型假设检验和均值比较，双因素方差与单因素方差的分析相似。

6. 重复测量方差分析

   1. 所谓重复测量方差分析，即受试者被测量不止一次；
   2. 传统的重复测量方差分析。该方法假设任意组内因子的协方差矩阵为球形，并且任意组内因子两水平间的方差之差都相等。但在现实中这种假设不可能满足，于是衍生了一系列备选方法：
      • 使用lme4包中的lmer()函数拟合线性混合模型（Bates，2005）；
      • 使用car包中的Anova()函数调整传统检验统计量以弥补球形假设的不满足（例如Geisser-Greenhouse校正）；
      • 使用nlme包中的gls()函数拟合给定方差-协方差结构的广义最小二乘模型（UCLA，2009）；
      • 用多元方差分析对重复测量数据进行建模（Hand，1987）。

7. 多元方差分析

   1. 当因变量（结果变量）不止一个时，可用多元方差分析（MANOVA）对它们同时进行分析
   2. 评估假设检验
   3. 稳健多元方差分析

8. 用回归来做 ANOVA