# Extra exercise for Pinheiro & Bates, Chapter 2

library(nlme)
library(lattice)
data(Machines)
fm1Machine <- lme(score~Machine, data=Machines, random=~1|Worker)
fm2Machine <- update(fm1Machine, random=~1|Worker/Machine)

#a
set.seed(251)
machineLRTsim <- simulate.lme(fm1Machine, fm2Machine, nsim=1000, method="ML")
plot(machineLRTsim, df=c(0,1), between=list(x=c(0,0,5)))
fm1<-machineLRTsim$alt$ML[,2]
fm1
names(fm1)<-NULL
fm2<-machineLRTsim$null$ML[,2]
names(fm2)<-NULL
min2logL <- 2*(fm1-fm2)

#b
summary(min2logL)
min2logL <- pmax(0,2*(fm1-fm2))
summary(min2logL)
mean(min2logL>0)

#c
n<-c(1:1000)
empprob<- (n-0.5)/1000
min2logL <- sort(min2logL)
chi1perc <- qchisq(empprob,1)
xyplot(chi1perc~min2logL)
chi1prob <- pchisq(min2logL,1)
xyplot(chi1prob~empprob)
xyplot(1-chi1prob~1-empprob)

#d
mixchiperc <- (empprob>0.5)*(qchisq(2*(abs(empprob-0.5)),1))
xyplot(mixchiperc~min2logL)
mixchiprob <- 0.5 + 0.5*pchisq(min2logL,1)
xyplot(mixchiprob~empprob)
xyplot(1-mixchiprob~1-empprob)
# another try with a .25:.75 split - closer fit but not perfect
mixchiperc <- (empprob>0.75)*(qchisq((pmin(.25,abs(empprob-0.75))/.25000001),1))
plot(mixchiperc,min2logL)
lines(chi1perc,chi1perc)


# Continuation
# Simulation of datasets from model and estimation in null and alternative models

fm1Machine <- lme(score~Machine, data=Machines, random=~1|Worker, method="ML")
fm2Machine <- update(fm1Machine, random=~1|Worker/Machine)
z1<-matrix(0,nrow=54,ncol=6)
for(i in 1:6) z1[,i] <- Machines$Worker==i
#z2<-matrix(0,nrow=54,ncol=18)
#numM <- 1 + (Machines$Machine=="B") + 2*(Machines$Machine=="C")
#for(i in 1:6){
#  for(j in 1:3) z2[,i+(j-1)*6] <- (Machines$Worker==i)*(numM==j)}

fixed1 <- Machines$score-residuals(fm1Machine, level=0)
# fixed2 <- Machines$score-residuals(fm2Machine, level=0)

set.seed(811)
simMachines <- Machines
lrtest <- c()
lmeControl (tolerance=1e-10,msTol=1e-10)
lmeControl (niterEM=0) # makes a substantial difference
for(i in 1:10)
{
#simMachines$score <- fixed2 + rnorm(54)*0.96158 + z1 %*% rnorm(6)*4.3645 + z2 %*% rnorm(18)*3.397
simMachines$score <- fixed1 + rnorm(54)*4.6833 + z1 %*% rnorm(6)*3.0951
fm1sim <- lme(score~Machine, data=simMachines, random=~1|Worker, method="ML", control=list(niterEM=0))
fm2sim <- update(fm1sim, random=~1|Worker/Machine)
lrtest <- c(lrtest,-2*(logLik(fm1sim)-logLik(fm2sim)))
names(lrtest)<-NULL
}
summary(lrtest)