library(nlme)
###########
# Split-plot---Whole plot completely Randomized
#The semiconductor experiment
######
D=read.table("semiconductor.dat",header=T)
D$ET=factor(D$ET)
D$WAFER=factor(D$WAFER)
D$POS=factor(D$POS)
GD=groupedData(RESISTANCE~1|WAFER,data=D,outer=~ET*POS)
plot(GD,outer=T)
formula(GD)
gsummary(GD,invar=TRUE)
GD$M=factor(interaction(GD$ET,GD$WAFER))
lmefit1=lme(RESISTANCE~ET*POS,data=GD,random=~1|M)
anova(lmefit1)
lmfit1=lm(RESISTANCE~ET*POS,data=GD)
anova(lmfit1)
summary(lmefit1)
summary(lmfit1) 
anova(lmefit1,lmfit1)
lmefit2=update(lmefit1, RESISTANCE~ ET + POS)
summary(lmefit2)
anova(lmefit2)
lmefit3=update(lmefit1,RESISTANCE~ POS)
summary(lmefit3)
anova(lmefit3)
lmfit2=lm(RESISTANCE~POS,data=GD)
summary(lmfit2)
anova(lmfit2)
anova(lmefit3,lmfit2) 
plot(lmefit3,formula=~resid(.)*fitted(.))

###########
# Split-Plot---Whole plot Randomized Complete Block
#The Agricultural Experiment
######
gsummary(Oats,invar=T)
plot(Oats,inner=~Variety)
Oats$nitro=factor(Oats$nitro)
lme1Oats=lme(yield~nitro*Variety,data=Oats,random=~1|Block/Variety)
anova(lme1Oats)
lme2Oats=update(lme1Oats, yield~nitro + Variety)
anova(lme2Oats)
summary(lme2Oats)
lme3Oats=update(lme1Oats, yield ~nitro)
summary(lme3Oats)
VarCorr(lme3Oats)
intervals(lme3Oats)

################
# Strip-Split Plot Design
#Cell Culture Bioassay Example
########
lme7=lme(logDens~sample*dilut,data=Assay,random=list(Block=pdBlocked(list(pdIdent(~1),pdIdent(~sample-1),pdIdent(~dilut-1)))))
intervals(lme7)
anova(lme7)
lm7=lm(logDens~sample*dilut,data=Assay)
anova(lme7)
anova(lme7,lm7)