9.6 重復(fù)測量方差分析
所謂重復(fù)測量方差分析，即受試者被測量不止一次。本節(jié)重點(diǎn)關(guān)注含一個組內(nèi)和一個組間因子的重復(fù)測量方差分析（這是一個常見的設(shè)計(jì)）。
以下為R語言實(shí)戰(zhàn)示例浓恶。基礎(chǔ)安裝包中的CO2數(shù)據(jù)集包含了北方和南方牧草類植物Echinochloa crus-galli（Potvin结笨，Lechowicz包晰，Tardif，1990）的寒冷容忍度研究結(jié)果炕吸，在某濃度二氧化碳的環(huán)境中伐憾，對寒帶植物與非寒帶植物的光合作用率進(jìn)行了比較。研究所用植物一半來自于加拿大的魁北克省赫模，另一半來自美國的密西西比州树肃。因變量是二氧化碳吸收量（uptake），單位為ml/L瀑罗，自變量是植物類型Type（魁北克VS密西西比州）和七種水平（95~1000 umol/m^2 sec）的二氧化碳濃度（conc）胸嘴。另外，Type是組間因子廓脆，conc是組內(nèi)因子筛谚。
CO2
## Plant Type Treatment conc uptake
## 1 Qn1 Quebec nonchilled 95 16.0
## 2 Qn1 Quebec nonchilled 175 30.4
## 3 Qn1 Quebec nonchilled 250 34.8
## 4 Qn1 Quebec nonchilled 350 37.2
## 5 Qn1 Quebec nonchilled 500 35.3
## 6 Qn1 Quebec nonchilled 675 39.2
## 7 Qn1 Quebec nonchilled 1000 39.7
## 8 Qn2 Quebec nonchilled 95 13.6
## 9 Qn2 Quebec nonchilled 175 27.3
## 10 Qn2 Quebec nonchilled 250 37.1
## 11 Qn2 Quebec nonchilled 350 41.8
## 12 Qn2 Quebec nonchilled 500 40.6
## 13 Qn2 Quebec nonchilled 675 41.4
## 14 Qn2 Quebec nonchilled 1000 44.3
## 15 Qn3 Quebec nonchilled 95 16.2
## 16 Qn3 Quebec nonchilled 175 32.4
## 17 Qn3 Quebec nonchilled 250 40.3
## 18 Qn3 Quebec nonchilled 350 42.1
## 19 Qn3 Quebec nonchilled 500 42.9
## 20 Qn3 Quebec nonchilled 675 43.9
## 21 Qn3 Quebec nonchilled 1000 45.5
## 22 Qc1 Quebec chilled 95 14.2
## 23 Qc1 Quebec chilled 175 24.1
## 24 Qc1 Quebec chilled 250 30.3
## 25 Qc1 Quebec chilled 350 34.6
## 26 Qc1 Quebec chilled 500 32.5
## 27 Qc1 Quebec chilled 675 35.4
## 28 Qc1 Quebec chilled 1000 38.7
## 29 Qc2 Quebec chilled 95 9.3
## 30 Qc2 Quebec chilled 175 27.3
## 31 Qc2 Quebec chilled 250 35.0
## 32 Qc2 Quebec chilled 350 38.8
## 33 Qc2 Quebec chilled 500 38.6
## 34 Qc2 Quebec chilled 675 37.5
## 35 Qc2 Quebec chilled 1000 42.4
## 36 Qc3 Quebec chilled 95 15.1
## 37 Qc3 Quebec chilled 175 21.0
## 38 Qc3 Quebec chilled 250 38.1
## 39 Qc3 Quebec chilled 350 34.0
## 40 Qc3 Quebec chilled 500 38.9
## 41 Qc3 Quebec chilled 675 39.6
## 42 Qc3 Quebec chilled 1000 41.4
## 43 Mn1 Mississippi nonchilled 95 10.6
## 44 Mn1 Mississippi nonchilled 175 19.2
## 45 Mn1 Mississippi nonchilled 250 26.2
## 46 Mn1 Mississippi nonchilled 350 30.0
## 47 Mn1 Mississippi nonchilled 500 30.9
## 48 Mn1 Mississippi nonchilled 675 32.4
## 49 Mn1 Mississippi nonchilled 1000 35.5
## 50 Mn2 Mississippi nonchilled 95 12.0
## 51 Mn2 Mississippi nonchilled 175 22.0
## 52 Mn2 Mississippi nonchilled 250 30.6
## 53 Mn2 Mississippi nonchilled 350 31.8
## 54 Mn2 Mississippi nonchilled 500 32.4
## 55 Mn2 Mississippi nonchilled 675 31.1
## 56 Mn2 Mississippi nonchilled 1000 31.5
## 57 Mn3 Mississippi nonchilled 95 11.3
## 58 Mn3 Mississippi nonchilled 175 19.4
## 59 Mn3 Mississippi nonchilled 250 25.8
## 60 Mn3 Mississippi nonchilled 350 27.9
## 61 Mn3 Mississippi nonchilled 500 28.5
## 62 Mn3 Mississippi nonchilled 675 28.1
## 63 Mn3 Mississippi nonchilled 1000 27.8
## 64 Mc1 Mississippi chilled 95 10.5
## 65 Mc1 Mississippi chilled 175 14.9
## 66 Mc1 Mississippi chilled 250 18.1
## 67 Mc1 Mississippi chilled 350 18.9
## 68 Mc1 Mississippi chilled 500 19.5
## 69 Mc1 Mississippi chilled 675 22.2
## 70 Mc1 Mississippi chilled 1000 21.9
## 71 Mc2 Mississippi chilled 95 7.7
## 72 Mc2 Mississippi chilled 175 11.4
## 73 Mc2 Mississippi chilled 250 12.3
## 74 Mc2 Mississippi chilled 350 13.0
## 75 Mc2 Mississippi chilled 500 12.5
## 76 Mc2 Mississippi chilled 675 13.7
## 77 Mc2 Mississippi chilled 1000 14.4
## 78 Mc3 Mississippi chilled 95 10.6
## 79 Mc3 Mississippi chilled 175 18.0
## 80 Mc3 Mississippi chilled 250 17.9
## 81 Mc3 Mississippi chilled 350 17.9
## 82 Mc3 Mississippi chilled 500 17.9
## 83 Mc3 Mississippi chilled 675 18.9
## 84 Mc3 Mississippi chilled 1000 19.9
CO2$conc <- as.factor(CO2$conc)
w1b1 <- subset(CO2, Treatment == 'chilled')
fit16 <- aov(uptake ~ (conc*Type) + Error(Plant/conc), w1b1)
summary(fit16)
##
## Error: Plant
## Df Sum Sq Mean Sq F value Pr(>F)
## Type 1 2667.2 2667.2 60.41 0.00148 **
## Residuals 4 176.6 44.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Plant:conc
## Df Sum Sq Mean Sq F value Pr(>F)
## conc 6 1472.4 245.40 52.52 1.26e-12 ***
## conc:Type 6 428.8 71.47 15.30 3.75e-07 ***
## Residuals 24 112.1 4.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(las=2)
par=(mar=c(10, 4, 4, 2))
with(w1b1, interaction.plot(conc, Type, uptake, type="b", col=c("red", "blue"), pch=c(16, 18), main="Interaction Plot for Plant Type and Concentraction"))
image.png
boxplot(uptake ~ Type*conc, data=w1b1, col=(c("gold", "green")), main = "Chilled Quebec and Mississippi Plants", ylab = "Carbon dioxide uptake rate (umol/m^2 sec)")
image.png
百度學(xué)習(xí)了一下磁玉，重復(fù)測量方差分析的R格式：model=aov(Y ~ B * W + Error(Subject/W))停忿，其中B是組間因子，W是組內(nèi)因子蚊伞，subject是實(shí)驗(yàn)對象的ID席赂。
下面自己構(gòu)建一個數(shù)據(jù)集df7演示學(xué)習(xí)，數(shù)據(jù)包括兩個N水平(N200,N300)时迫，每個氮水平下分別在5個時期測定兩個品種5株的biomass颅停。這里N是組間因子，T是組內(nèi)因子掠拳。
df7 <- read.table(file = "D:/Documents/R wd/df7.csv", header = T, sep = ",") # 數(shù)據(jù)導(dǎo)入癞揉。
df7 # 查看數(shù)據(jù)。
## Plant N t1 t2 t3 t4 t5
## 1 Q1 200 0.1 4.5 9.9 11.0 15.2
## 2 Q1 200 0.2 4.2 8.1 12.2 15.0
## 3 Q1 200 0.1 4.0 8.3 11.8 15.1
## 4 Q1 200 0.2 4.7 7.2 12.8 16.0
## 5 Q1 200 0.1 4.5 7.1 11.5 15.4
## 6 Q2 200 0.3 4.6 8.0 12.7 15.8
## 7 Q2 200 0.1 4.8 9.6 12.5 15.4
## 8 Q2 200 0.3 5.3 9.2 13.0 14.8
## 9 Q2 200 0.4 5.5 9.5 13.5 14.9
## 10 Q2 200 0.2 5.2 10.0 13.8 15.0
## 11 M1 300 0.8 8.6 15.0 19.5 23.5
## 12 M1 300 0.9 8.8 14.7 19.6 23.8
## 13 M1 300 0.7 9.0 16.5 19.4 23.4
## 14 M1 300 0.8 8.4 15.2 18.8 25.1
## 15 M1 300 0.7 8.2 16.5 18.6 25.0
## 16 M2 300 0.6 8.1 14.8 18.0 24.8
## 17 M2 300 0.9 7.9 15.8 19.0 24.7
## 18 M2 300 1.0 7.5 16.8 18.7 24.3
## 19 M2 300 0.8 7.8 14.6 18.6 24.9
## 20 M2 300 0.9 7.7 14.2 18.7 25.0
attach(df7) # 將數(shù)據(jù)集df7加入搜索路徑溺欧。
library(reshape2) # 調(diào)用reshape2包喊熟。
w2b2 <- melt(df7,id = c("N","Plant"),variable.name = "T", value.name="Biomass") # 寬格式數(shù)據(jù)變?yōu)殚L格式。
w2b2 # 查看數(shù)據(jù)姐刁。
## N Plant T Biomass
## 1 200 Q1 t1 0.1
## 2 200 Q1 t1 0.2
## 3 200 Q1 t1 0.1
## 4 200 Q1 t1 0.2
## 5 200 Q1 t1 0.1
## 6 200 Q2 t1 0.3
## 7 200 Q2 t1 0.1
## 8 200 Q2 t1 0.3
## 9 200 Q2 t1 0.4
## 10 200 Q2 t1 0.2
## 11 300 M1 t1 0.8
## 12 300 M1 t1 0.9
## 13 300 M1 t1 0.7
## 14 300 M1 t1 0.8
## 15 300 M1 t1 0.7
## 16 300 M2 t1 0.6
## 17 300 M2 t1 0.9
## 18 300 M2 t1 1.0
## 19 300 M2 t1 0.8
## 20 300 M2 t1 0.9
## 21 200 Q1 t2 4.5
## 22 200 Q1 t2 4.2
## 23 200 Q1 t2 4.0
## 24 200 Q1 t2 4.7
## 25 200 Q1 t2 4.5
## 26 200 Q2 t2 4.6
## 27 200 Q2 t2 4.8
## 28 200 Q2 t2 5.3
## 29 200 Q2 t2 5.5
## 30 200 Q2 t2 5.2
## 31 300 M1 t2 8.6
## 32 300 M1 t2 8.8
## 33 300 M1 t2 9.0
## 34 300 M1 t2 8.4
## 35 300 M1 t2 8.2
## 36 300 M2 t2 8.1
## 37 300 M2 t2 7.9
## 38 300 M2 t2 7.5
## 39 300 M2 t2 7.8
## 40 300 M2 t2 7.7
## 41 200 Q1 t3 9.9
## 42 200 Q1 t3 8.1
## 43 200 Q1 t3 8.3
## 44 200 Q1 t3 7.2
## 45 200 Q1 t3 7.1
## 46 200 Q2 t3 8.0
## 47 200 Q2 t3 9.6
## 48 200 Q2 t3 9.2
## 49 200 Q2 t3 9.5
## 50 200 Q2 t3 10.0
## 51 300 M1 t3 15.0
## 52 300 M1 t3 14.7
## 53 300 M1 t3 16.5
## 54 300 M1 t3 15.2
## 55 300 M1 t3 16.5
## 56 300 M2 t3 14.8
## 57 300 M2 t3 15.8
## 58 300 M2 t3 16.8
## 59 300 M2 t3 14.6
## 60 300 M2 t3 14.2
## 61 200 Q1 t4 11.0
## 62 200 Q1 t4 12.2
## 63 200 Q1 t4 11.8
## 64 200 Q1 t4 12.8
## 65 200 Q1 t4 11.5
## 66 200 Q2 t4 12.7
## 67 200 Q2 t4 12.5
## 68 200 Q2 t4 13.0
## 69 200 Q2 t4 13.5
## 70 200 Q2 t4 13.8
## 71 300 M1 t4 19.5
## 72 300 M1 t4 19.6
## 73 300 M1 t4 19.4
## 74 300 M1 t4 18.8
## 75 300 M1 t4 18.6
## 76 300 M2 t4 18.0
## 77 300 M2 t4 19.0
## 78 300 M2 t4 18.7
## 79 300 M2 t4 18.6
## 80 300 M2 t4 18.7
## 81 200 Q1 t5 15.2
## 82 200 Q1 t5 15.0
## 83 200 Q1 t5 15.1
## 84 200 Q1 t5 16.0
## 85 200 Q1 t5 15.4
## 86 200 Q2 t5 15.8
## 87 200 Q2 t5 15.4
## 88 200 Q2 t5 14.8
## 89 200 Q2 t5 14.9
## 90 200 Q2 t5 15.0
## 91 300 M1 t5 23.5
## 92 300 M1 t5 23.8
## 93 300 M1 t5 23.4
## 94 300 M1 t5 25.1
## 95 300 M1 t5 25.0
## 96 300 M2 t5 24.8
## 97 300 M2 t5 24.7
## 98 300 M2 t5 24.3
## 99 300 M2 t5 24.9
## 100 300 M2 t5 25.0
str(w2b2) # 查看數(shù)據(jù)結(jié)構(gòu)芥牌。
## 'data.frame': 100 obs. of 4 variables:
## $ N : int 200 200 200 200 200 200 200 200 200 200 ...
## $ Plant : chr "Q1" "Q1" "Q1" "Q1" ...
## $ T : Factor w/ 5 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Biomass: num 0.1 0.2 0.1 0.2 0.1 0.3 0.1 0.3 0.4 0.2 ...
w2b2$N <- as.factor(w2b2$N) # df7數(shù)據(jù)集中N轉(zhuǎn)為因子。
w2b2$T <- as.factor(w2b2$T) # df7數(shù)據(jù)集中T轉(zhuǎn)為因子聂使。
fit17 <- aov(Biomass ~ T*N + Error(Plant/T), w2b2) # 構(gòu)建模型壁拉。
summary(fit17) # 返回結(jié)果谬俄。
##
## Error: Plant
## Df Sum Sq Mean Sq F value Pr(>F)
## N 1 697.0 697.0 267.9 0.00371 **
## Residuals 2 5.2 2.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Plant:T
## Df Sum Sq Mean Sq F value Pr(>F)
## T 4 4643 1160.7 1366.71 2.28e-11 ***
## T:N 4 219 54.7 64.36 4.03e-06 ***
## Residuals 8 7 0.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 80 25.16 0.3145
par(las=2) # las 參數(shù)控制x軸和y軸的刻度線上的標(biāo)簽與兩條軸的方向，可選值為0,1,2,3弃理，默認(rèn)值0表示總是平行于坐標(biāo)軸溃论；1表示總是水平方向；2表示總是垂直于坐標(biāo)軸痘昌；3表示總是垂直方向蔬芥。
par(mar=c(10,4,4,2)) # 設(shè)置圖形空白邊界行數(shù)，mar=c(bottom,left,top,right)控汉，缺省為mar=c(5.1,4.1,4.1,2.1)
with(w2b2, interaction.plot(T,N,Biomass, type = "b", col = c("red","blue"), pch = c(16,18), main = "Interaction Plot for Nitrogen and Time")) # 結(jié)果可視化笔诵。
image.png
boxplot(Biomass ~ N*T, data=w2b2, col=(c("gold","green")),
main="Nitrogen and Time",
ylab="Biomass") # 結(jié)果可視化。
image.png
參考資料：
- 《R語言實(shí)戰(zhàn)》（中文版）姑子，人民郵電出版社乎婿，2013.
