library(nlme)

# Section 2.4
data(Orthodont)
OrthoFem <- Orthodont[Orthodont$Sex == "Female",]
fm1OrthF <- lme(distance ~ age, data=OrthoFem, random=~1|Subject)
fm2OrthF <- update(fm1OrthF, random= ~ age|Subject)
orthLRTsim <- simulate.lme(fm1OrthF, fm2OrthF, nsim=1000)
plot(orthLRTsim, df=c(1,2))

data(Machines)
fm1Machine <- lme(score~Machine, data=Machines, random=~1|Worker)
fm2Machine <- update(fm1Machine, random=~1|Worker/Machine)
machineLRTsim <- simulate.lme(fm1Machine, fm2Machine, nsim=1000)
plot(machineLRTsim, df=c(0,1), layout=c(4,1), between=list(x=c(0,0.5)))

PBIB <- groupedData( y ~ as.factor(tx) | block,
data=read.table("c:/data.avc/teaching/vhm881/pinbat11.dat", header=TRUE))
PBIB$tx <- as.factor(PBIB$tx)
pbibLRTsim <- simulate.lme(m1=list(fixed=y~1, data=PBIB,
random=~1|block), m2=list(fixed=y ~ tx), method="ML", nsim=1000)
plot(pbibLRTsim, df=c(14,16,18), weights=FALSE)

# analyses on different datasets to see effect of df method for fixed methods
# see also corresponding SAS proc mixed code
anova(fm2Machine)
MachinesUnbal <- Machines[-c(2,3,6,8,9,12,19,20,27,33),]
fm2MachineU <- lme(score~Machine, data=MachinesUnbal, random=~1|Worker/Machine)
fm2MachineU
anova(fm2MachineU)

fm1Rail.lme <- lme(travel ~1, data=Rail, random= ~ 1|Rail)
anova(fm1Rail.lme)
# last example seems very strange