Mixed Models - RCBD

### Randomized Complete Block Design
library(nlme)

rcb <- read.csv("rcb.csv")
rcb$ingot <- as.factor(rcb$ingot)
rcb$metal <- as.factor(rcb$metal)

### lme
rcbFit <- lme(pres~ metal, random=~1|ingot, data=rcb)

### manually
t<-3; b<-7

mean_b<-tapply(rcb[,3],rcb[,1],mean)
mean_t<-tapply(rcb[,3],rcb[,2],mean)[t:1]

msa<-b*sum((mean_t-mean(mean_t))^2)/(t-1)
msb<-a*sum((mean_b-mean(mean_b))^2)/(b-1)
mse<-sum((y-rep(mean_t,b)-rep(mean_b,rep(t,b))+mean(y))^2)/(t-1)/(b-1)
F<-msa/mse
p_value <- 1-pf(F,(t-1),(t-1)*(b-1))

sb<-(msb-mse)/t # sigma_b

/* data & SAS code */
data rcb;
   input ingot metal $ pres @@;
   datalines;
1 n 67.0  1 i 71.9  1 c 72.2
2 n 67.5  2 i 68.8  2 c 66.4
3 n 76.0  3 i 82.6  3 c 74.5
4 n 72.7  4 i 78.1  4 c 67.3
5 n 73.1  5 i 74.2  5 c 73.2
6 n 65.8  6 i 70.8  6 c 68.7
7 n 75.6  7 i 84.9  7 c 69.0
run;

proc mixed data=rcb;
   class ingot metal;
   model pres=metal;
   random ingot;
run;