單因素協方差分析

單因素協方差分析是單因素方差分析中包含一個或多個定量的協變量。

• 擬合單因素協方差分析

#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) ，生成的圖形將容許斜率和截距項依據組別而發生變化，這對可視化那些違揹回歸斜率同質性的實例很是有用