# Notes and code for practical parts of home assignment 2

# Part 4
library(nlme)
data(Pixel)
fm1Pixel <- lme(pixel ~ day + I(day^2), data=Pixel, random=list(Dog=~day, Side=~1))
summary(fm1Pixel)
# estimates agree well with MLwiN
residuals(fm1Pixel) #lowest level residuals (level=2)
ranef(fm1Pixel,level=2)
ranef(fm1Pixel,level=1)
# predicted random effects (residuals) also agree well with MLwiN

# Part 5
residuals(fm1Pixel)/residuals(fm1Pixel,type="pearson")
ranef(fm1Pixel,level=2)/ranef(fm1Pixel,level=2,standard=TRUE)
ranef(fm1Pixel,level=1)/ranef(fm1Pixel,level=1,standard=TRUE)
# standard errors at all 3 levels are constant (across obs)
# don't agree with any of the standard errors in MLwiN

# Part 6
fm1Pixel<- lme(pixel~day+I(day^2),data=Pixel,random=~1|Dog/Side)
# simpler model used
Dz1<- matrix(0,nrow=102,ncol=10)
  for (i in 1:10) Dz1[,i]<-Pixel$Dog==i
numPixel<- (Pixel$Side=="R")+2*(Pixel$Side=="L")
Dz2<- matrix(0,nrow=102,ncol=20)
  for (i in 1:10) {
     for (j in 1:2) Dz2[,i+(j-1)*10]<- (Pixel$Dog==i)*(numPixel==j)
     }
# The fixed effects
fixed<- Pixel$pixel-residuals(fm1Pixel,level=0)

S<-c()
S1<-c()
S2<-c()
MRR<-c()
MRR1<-c()
MRR2<-c()
SDRR<-c()
SDRR1<-c()
SDRR2<-c()
set.seed(666)

for (i in 1:100) {
  simPixel<- fixed+rnorm(102)*12.91646 + Dz1%*%rnorm(10)*  22.82149+ Dz2%*%rnorm(20)*15.70170 
  fmPixSim<-lme(simPixel~day+I(day^2),data=Pixel,random=~1|Dog/Side)
  res<-residuals(fmPixSim,level=0,type="pearson")
  S<-c(S,shapiro.test(res)$p.value)
  MRR<-c(MRR,mean(res))
  SDRR<-c(SDRR,sd(res))
  res1<-ranef(fmPixSim,level=1,standard=TRUE)$"(Intercept)"
  res2<-ranef(fmPixSim,level=2,standard=TRUE)$"(Intercept)"
  SDRR1<-c(SDRR1,sd(res1))
  MRR1<-c(MRR1,mean(res1))
  SDRR2<-c(SDRR2,sd(res2))
  MRR2<-c(MRR2,mean(res2))
  S1<-c(S1,shapiro.test(res1)$p.value)
  S2<-c(S2,shapiro.test(res2)$p.value)
  }

summary(MRR)
summary(SDRR)
summary(MRR1)
summary(SDRR1)
summary(MRR2)
summary(SDRR2)
summary(S>0.05)
summary(S1>0.05)
summary(S2>0.05)