最近接觸到滿多次雙因子變異數分析,想說就記錄一下了解這個統計方法的過程吧
還是先說明什麼是變異數分析(analysis of variance, ANOVA),它是一種分析類別自變數對數值應變數影響的統計方法。
舉例說給雞添加兩種營養物質(A, B),並分析這兩種處裡各自對體增重的影響。
添加營養物質就是類別自變數; 體增重為數值應變數。
跳回雙因子變異數分析吧,就是指兩項類別自變數對數值應變數的影響。
從剛才例子中衍生,原本只有一項類別自變數(營養物質),又多加一項環境溫度(20, 25, 30℃)
也就是一隻雞同時接受兩種因子影響(如餵飼A物質且處於25℃),組合總共有2X3= 6組處理
接下用R進行演示一邊說明如何使用two-way anova
對R有興趣可以搭配程式碼,一般直接看統計結果即可。
資料呈現如下表,每組有10隻雞的體增重(由於只是演示,請忽略體增重的合理性)
20°C | 25°C | 30°C | |
---|---|---|---|
A | 4.2 | 16.5 | 23.6 |
11.5 | 16.5 | 18.5 | |
7.3 | 15.2 | 33.9 | |
5.8 | 17.3 | 25.5 | |
6.4 | 22.5 | 32.5 | |
10 | 17.3 | 32.5 | |
11.2 | 13.6 | 26.7 | |
11.2 | 14.5 | 21.5 | |
5.2 | 18.8 | 23.3 | |
7 | 15.5 | 29.5 | |
B | 15.2 | 19.7 | 25.5 |
21.5 | 23.3 | 26.4 | |
17.6 | 23.6 | 22.4 | |
9.7 | 26.4 | 24.5 | |
14.5 | 20 | 24.8 | |
10 | 25.2 | 30.9 | |
8.2 | 25.8 | 26.4 | |
9.4 | 21.2 | 27.3 | |
16.5 | 14.5 | 29.4 | |
9.7 | 27.3 | 23 |
一、輸入資料
setwd("C:/R/biostatistics")
data <- read.csv("data_twoway_anova.csv", header=T)
data$temp <- as.factor(data$temp)
```
二、描述性統計(盒鬚圖、平均值及標準差)
```R
#boxplot
boxplot(gain ~ supp + temp, col=c("red", "orange"), ylab="Body weight gain", xlab="supplement and temperature", data=data)
libraray(ggplot2)
ggplot(data, aes(x=supp, y=gain, col=temp)) + geom_boxplot()
```
boxplot可使用內建函式(boxplot)或用ggplot繪圖
![](boxplot.png)
```R
# class by supp and temp, and then count means and standard error
library(reshape); library(agricolae)
mystats <- function(x) {c(n=length(x), mean=mean(x, na.rm=T), sd=sd(x))}
dfm <- melt(data, measure.vars="gain", id.vars=c("supp", "temp"))
cast(dfm, supp + temp + variable ~ ., mystates)
六個組別各自的平均值(mean)與標準差(sd)
supp temp variable n mean sd
1 A 20 gain 10 7.98 2.746634
2 A 25 gain 10 16.77 2.515309
3 A 30 gain 10 26.14 4.797731
4 B 20 gain 10 13.23 4.459709
5 B 25 gain 10 22.70 3.910953
6 B 30 gain 10 26.06 2.655058
三、主效果變異數分析
先不考慮交互作用,只分析營養物質(supp)和環境溫度(temp)個別對體增重(gain)的影響效果
# main effect ANOVA with no interaction
model_aov <- aov(data=data, gain ~ supp + temp)
summary(model_aov)
model_aov$cofficients
# if we don't care supp and temp, mean of gain is 12.455
# if suppOJ = 0. and then suppVC = -3.7