library(nlme)
library(multcomp)
#############
#Oneway Random Effects Model
########
#The Pulp Example
# Paper Brightness Depending on Operator?

#REML Estimation
D=read.table(file="pulp.txt")
plot(lm1)
lme1=lme(bright~1,random=~1|operator,data=D)
summary(lme1)
intervals(lme1)
anova(lme1)
plot(lme1)
ICC1=0.2609286^2/(0.2609286^2+0.32596^2)
#ML Estimation
lme2=lme(bright~1,random=~1|operator,data=D,method="ML")
summary(lme2)
intervals(lme2)
anova(lme2)
plot(lme2)
ICC2= 0.2138924^2/(0.2138924^2+0.3259601^2)

# Fixed Effects Model
lm1=lm(bright~operator,data=D)
summary(lm1)

anova(lme1,lm1)

###############
#Ranomized Block Design 
####################
#The Ergonomics Experiment (Sigle Fixed Factor)
plot(ergoStool)
plot.design(ergoStool)
attach(ergoStool)
interaction.plot(Type,Subject,effort)
detach()
ergoStool$Type=factor(ergoStool$Type,levels=c("T4","T1","T2","T3"))
lmeb=lme(effort~Type,data=ergoStool,random=~1|Subject)
summary(lmeb)
anova(lmeb)
#Alternative EQUIVALENT form of specifiying the variance component
lmeb1=lme(effort~Type,data=ergoStool, random=list(Subject=pdIdent(~1)))
summary(lmeb1)
#Testing the significance of the sbjuect variance component.
lmb=lm(effort~Type,data=ergoStool)
anova(lmeb,lmb)

# Multiple Pairwise Comparisons
mpcb=glht(lmeb,linfct=mcp(Type="Tukey"))
summary(mpcb)
plot(confint(mpcb))

##################
# Replicated Block Design
#########
#Two Way Mixed Effects Model with Interaction (Replicated Block Design)
# The Machines Example
#Machines Example---Four EQUIVALENT Specification of the variance components
plot(Machines)
lme3=lme(score~Machine, data=Machines, random=list(Worker=pdBlocked(list(pdIdent(~1),pdIdent(~Machine-1)))))
lme4=lme(score~Machine,data=Machines, random=~1|Worker/Machine)
lme5=lme(score~Machine, data=Machines, random=list(Worker=pdCompSymm(~Machine-1))) #avoid this parametrization for the variance components
lme6=lme(score~Machine,data=Machines, random=list(Worker=~1, Machine=~1))




#########
#Nested Design
#########
#The Wafer Example
lme2=lme(Thickness~1,data=Oxide, random=list(Lot=~1, Wafer=~1))
summary(lme2)
intervals(lme2)


#Two way design two Fixed Factors with a random block effects and two random factors nested within block
#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)