单因素协方差分析

单因素协方差分析是单因素方差分析中包含一个或多个定量的协变量。

• 拟合单因素协方差分析

#multcomp包中litter数据集，怀孕小数分为四组，每组接受不一样剂量（0、五、50、500）的要药物处理。产下幼崽的体重均值为因变量
#怀孕时长为协变量
> data(litter, package="multcomp") #注意从包中直接导入数据的方式
> attach(litter)
The following objects are masked from litter (pos = 3):

    dose, gesttime, number, weight

> table(dose)   #每种剂量下幼崽数并不相同
dose
  0   5  50 500 
 20  19  18  17 
> 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


> 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 **   #ANOCOVA 的 F检验代表，怀孕时间与幼崽出生体重相关
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

因为使用了协变量，可能想获取调整的组均值，即 去除协变量效应后的组均值。

可以使用 effects 包中的 effects() 函数来计算调整的均值

> library(effects)
> effect("dose",fit)

 dose effect
dose
       0        5       50      500 
32.35367 28.87672 30.56614 29.33460

• 分析哪一种剂量对幼崽的影响

剂量 F 检验虽然代表了不一样的处理方式对幼崽的体重均值不一样，但没法告知咱们哪一种处理方式与其余方式不一样，一样的使用 multcomp包来对全部均值进行成对比（multcomp包还可用来家检验用户自定义的均值假设）

#检验未用药条件与其余三种用药条件影响是否不一样感兴趣
> library(multcomp)
> contrast <- rbind("no drug vs. drug" = c(3,-1,-1,-1))  #注意这里的设置方式，对照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 *  #p<0.05显著，接受H0，获得未用药组比其余用药条件下的出生体重要高，其余的对照可用 rbind() 函数添加（详见help(glht)）
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

 

• 评估检验的假设条件

a、与ANOVA相同，都须要验证变量正态性和方差相同假设，方法都一致

b、ANOCVA检验，还假定回归斜率相同，本例中，假定四个处理组经过怀孕时间来预测出生体重的回归斜率都相同。

ANCOVA模型包含怀孕时间*剂量的交互项时，可对 回归斜率的同质性进行检验。交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平

#检验回归斜率的同质性
> library(multcomp)
> 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   #交互项 p>0.05不显著，接受好H0，支持了斜率相等的假设
Residuals     66 1069.4   16.20                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#斜率相等的假设不成立，能够尝试切换协变量或因变量，或使用能对每一个斜率独立解释的模型，或使用不须要假设回归斜率同质性的非参数ANCOVA方法，sm 包中的 sm.ancova()函数为后者提供一个例子

• 结果可视化

HH包中的 ancova() 函数能够绘制因变量、协变量与因子之间的关系图

> library(HH)
> ancova(weight ~ gesttime + dose, data=litter)
Analysis of Variance Table

Response: weight
          Df  Sum Sq Mean Sq F value   Pr(>F)   
gesttime   1  134.30 134.304  8.0493 0.005971 **
dose       3  137.12  45.708  2.7394 0.049883 * 
Residuals 69 1151.27  16.685                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

从图可看出，用怀孕时间来预测出生体重的回归线性互相平行，只是截距项不一样。随着怀孕时间增长，幼崽出生体重也会增长，另外，还能够看到0剂量截距项最大，5剂量组截距项最小。因为上面设置的直线会保持平行，若用 ancova(weight ~ gesttime*dose) ，生成的图形将容许斜率和截距项依据组别而发生变化，这对可视化那些违背回归斜率同质性的实例很是有用