## R-code for Random Effects Notes
## AHW3
## April 28, 2008
rm(list = ls())

## one-factor random effects model
y <- c(14.0, 14.1, 14.2, 14.0, 14.1,
	   13.9, 18.8, 13.9, 14.0, 14.0,
	   14.1, 14.2, 14.1, 14.0, 13.9, 
	   13.6, 13.8, 14.0, 13.9, 13.7, 
	   13.8, 13.6, 13.9, 13.8, 14.0)
	   
loom <- factor(rep(1:5, each=5))

## anova
fit.lm <- lm(y~ as.factor(loom))
anova(fit.lm)

sigma.sq.hat <-  anova(fit.lm)[2,3]
sigma.sq.tau.hat <- (anova(fit.lm)[1,3] - anova(fit.lm)[2,3])/5

sigma.sq.hat
sigma.sq.tau.hat

## lme
library(nlme)
fit.me<-lme(fixed=y~ +1, random=~1|as.factor(loom))
fit.me


## two-factor mixed effects model
## pos is random
## temp is fixed
rm(list = ls())

y <- c(570, 1063, 565,
       565, 1080, 510,
	   583, 1043, 590,
	   528, 988, 526,
	   547, 1026, 538,
	   521, 1004, 532)

pos <- as.factor(rep(1:2, each=9))
temp <- as.factor(rep(c(800,825, 850), 6)) 

## anova
fit.lm <- lm(y ~ as.factor(temp) +  as.factor(pos) + as.factor(temp):as.factor(pos))

## get the Mean Squares
MSA <- anova(fit.lm)[1,3]
MSB <- anova(fit.lm)[2,3]
MSAB <- anova(fit.lm)[3,3]
MSE <- anova(fit.lm)[4,3]

MSA
MSB
MSAB
MSE

## calculate the observed F-statistics
F.A.obs <- MSA/MSAB
F.B.obs <- MSB/MSE
F.AB.obs <- MSAB/MSE

F.A.obs
F.B.obs
F.AB.obs

## calculate the p-values
p.value.A <- 1-pf(F.A.obs,  3-1, (3-1)*(2-1))
p.value.B <- 1-pf(F.B.obs, 2-1, 3*2*(3-1))
p.value.AB <- 1-pf(F.AB.obs,  (3-1)*(2-1), 3*2*(3-1))

p.value.A
p.value.B
p.value.AB