실험계획: Fish Tank

[안재형]

Casella의 Statistical Design이라는 책에 예1.6입니다.

http://www.stat.ufl.edu/archived/casella/StatDesign/WebDataSets/FishTank.txt

3가지 먹이중 어느 것이 물고기 몸무게 증가에 효과적인지 보려는 실험으로,

어항 12개를 3가지 먹이에 무작위로 배정하는데, 

각각 어항에 6마리에 물고기가 있습니다.

여기서 중요한 것은 먹이에 무작위로 배정된 건 "어항"이지 "물고기"가 아니어서 experimental unit이 어항이라는 점입니다.

이 경우 보통 어항이 먹이에 nested되었다고 합니다.

읽어서 Diet와 Tank를 factor로 바꿔주고 (안그러면 값의 크기가 의미가 있는 numeric으로 처리됩니다)

FishTank=read.table("http://www.stat.ufl.edu/archived/casella/StatDesign/WebDataSets/FishTank.txt",header=TRUE)

FishTank$Diet=as.factor(FishTank$Diet)

FishTank$Tank=as.factor(FishTank$Tank)

보통같이 2-way anova로 처리하면 다음과 같은데, 
Diet의 F-값 206.606은 10623.9/51.4로 Residuals의 MSE로 나눈 값입니다. 
그런데 이렇게 하면 안되고...

> anova(lm(WtGain~Diet+Tank,data=FishTank))
Analysis of Variance Table

Response: WtGain
          Df  Sum Sq Mean Sq F value Pr(>F)    
Diet       2 21247.7 10623.9 206.606 <2e-16 ***
Tank       9   107.8    12.0   0.233 0.9883    
Residuals 60  3085.3    51.4    

Tank가 Diet에 무작위로 배정되었으므로, 
맞는 F-값은 Tank의 MSE인 12으로 나눈 10623.9/12=885.325를 사용해야됩니다. 
다음과 같이 하면 동일한 값을 얻습니다.

> summary(aov(WtGain~Diet+Error(Tank/Diet),data=FishTank))

Error: Tank
          Df Sum Sq Mean Sq F value   Pr(>F)    
Diet       2  21248   10624   886.8 4.62e-11 ***
Residuals  9    108      12                     

### SAS로 하려면 다음과 같이 하는데 값이 다르게 나와서 뭘 잘못이해했나 했더니 Tank(Diet)의 분산추정치가 0인걸 보니 데이터가 불안정하네요. simulation한 데이터를 사용하니 동일하게 나오네요.

proc mixed data=fishtank;

  class diet tank;

  model wtgain = diet / s;

  random tank(diet);  /*  여기에 tank가 diet에 nested되었다고 지정하면 됩니다 */

run;