#' ---
#' title: "Coding Examples"
#' author: "MATH 4939"
#' date: "2026-03-15"
#' geometry: top=.5in,left=.5in,right=.5in,bottom=.5in,paperwidth=6in,paperheight=4in
#' output: pdf_document
#'   toc: false
#'     number_sections: false
#'     toc_depth: 3
#'   highlight: tango
#' header-includes:
#'   - \pagenumbering{arabic}
#' ---
#' \pagenumbering{gobble}
#' \pagenumbering{arabic}
#' \raggedright
#' 
#' # ML, REML, using dvar or raw variable in RE model
#'
knitr::opts_chunk$set(error=TRUE)
#' 
library(spida2)
library(latticeExtra)
library(nlme)

hsp <- subset(hs, Sector == "Public")
hsp <- hsp[order(hsp$ses), ]

model_1 <- lme(mathach ~ ses, hsp, random = ~ 1 + ses | school)

model_1 <- lme(mathach ~ ses, hsp, random = ~ 1 + ses | school, 
               control = list(msVerbose = TRUE))

model_1 <- lme(mathach ~ ses, hsp, random = ~ 1 + ses | school, 
               control = list(msMaxIter = 1000, msVerbose = TRUE))

model_1 <- lme(mathach ~ ses, hsp, random = ~ 1 + ses | school, 
               control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = T))
summary(model_1)
#' 
#' pt. of minimum variance is at
#' $$-g_{01}/g_{11}$$ 
#' or using $r_{01}$, $\sqrt{g_{00}}, \sqrt{g_{11}}$:    
#' $$r_{01}\frac{\sqrt{g_{00}}}{\sqrt{g_{11}}}
#'
1 *  2.2066885 / 0.6588577


model_1b <- lme(mathach ~ ses, hsp, random = ~ 1 + I(ses + 3.34) | school, 
               control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = T))
summary(model_1b)

hsp$model_1b.blup <- predict(model_1b, level = 1)
xyplot(model_1b.blup ~ ses, hsp, group = school, type = 'l')

# 
# what if we had something other than 3.34
# 
model_1c <- lme(mathach ~ ses, hsp, random = ~ 1 + I(ses + 3.6) | school, 
                control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = T))
model_1c <- lme(mathach ~ ses, hsp, random = ~ 1 + I(ses + 3.5) | school, 
                control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = T))

hsp$model_1c.blup <- predict(model_1c, level = 1)
xyplot(model_1c.blup ~ ses, hsp, group = school, type = 'l')  

AIC(model_1b, model_1c)  # equivalent models except numerically

#' 
#' Does it make sense to suppose that random variation in school lines
#' amounts to radiation from a point?
#' 

model_2 <- lme(mathach ~ ses, hsp, random = ~ 1 + dvar(ses, school) | school)

model_2 <- lme(mathach ~ ses, hsp, random = ~ 1 + dvar(ses, school) | school,
               control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = T))
summary(model_2)

model_2b <- lme(mathach ~ ses, hsp, random = ~ 1 + I(dvar(ses, school)+1.5) | school,
               control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = T))
summary(model_2b)

hsp$model_2b.blup <- predict(model_2b, level = 1)
xyplot(model_2b.blup ~ ses, hsp, group = school, type = 'l')  

#' 
#' Compare raw with dvar models
#' 
#' Same FE models so we can compare with AIC using REML (default) models
#' 
AIC(model_1c, model_2b)

#' 
#' Different FE
#' 

model_3 <- lme(mathach ~ ses + cvar(ses, school), hsp, random = ~ 1 + ses | school,
               control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = T))

#' surprise

summary(model_3)
hsp$model_3.blup <- predict(model_3, level = 1)
xyplot(model_3.blup ~ ses, hsp, group = school, type = 'l')  

xyplot(model_3.blup ~ ses, hsp, group = school, type = 'l')  +
  xyplot(cvar(mathach, school) ~ cvar(ses,school), hsp, group = school, pch = 16)

#'
#' have a look at BLUEs
#'
model_3lm <-  lm(mathach ~ ses * factor(school),hsp)
hsp$model_3.blue <- predict(model_3lm,hsp)

xyplot(model_3.blup ~ ses, hsp, group = school, type = 'l')  +
  xyplot(cvar(mathach, school) ~ cvar(ses,school), hsp, group = school, pch = 16) +
  xyplot(model_3.blue ~ ses, hsp, group = school, type = 'l', lty = 3)

#' 
#' small example of OOP:
#' 
form <- function(fit,...) UseMethod('form')
form.lme <- function(fit,...)  {
  paste(c(as.character(formula(fit))[c(2,1,3)], '  random = ',
          as.character(fit$call$random)[c(1,2)]), collapse = ' ')
}
form.lm <- function(fit,...) {
  paste(c(as.character(formula(fit))[c(2,1,3)]), collapse = ' ')
} 

form(model_1b)
form(model_1c)
form(model_2b)
form(model_3)

clist <- list(msMaxIter = 1000, msVerbose = TRUE, returnObject = TRUE)

model_4 <- lme(mathach ~ ses + cvar(ses, school), hsp, random = ~ 1 + dvar(ses, school) | school,
               control = clist)

summary(model_4)


form(model_1b)
form(model_1c)
form(model_2b)
form(model_3)
form(model_4)

#' 
#' How to compare model_3 and model_4
#' 
#' - Can't use Wald tests
#' - REML fits, same FE different RE
#' - So we can use Likelihood
#' - Models not nested
#' - So compare with AIC, or BIC
#' 
AIC(model_3, model_4)
#' 
#' How to compare model_3 and model_1c
#' 
#' - Different FEs, equivalent REs
#' - Could use Wald tests
#' - Can's use Likelihood for FE since REML fits
#' - Let's compare Wald with LRT on ML fits
#' 
AIC(model_3, model_1c)
wald(model_3)
wald(model_3, 'cvar')
#
model_1c.ml <- update(model_1c, method = 'ML')
model_3.ml <- update(model_3, method = 'ML')

# LRT

anova(model_1c.ml, model_3.ml)

#'
#' Finally 2 models equivalent to 3 and 4
#' 


model_5 <- lme(mathach ~ dvar(ses,school) + cvar(ses, school), hsp, random = ~ 1 + ses | school)
model_6 <- lme(mathach ~ dvar(ses,school) + cvar(ses, school), hsp, random = ~ 1 + dvar(ses, school) | school)

AIC(model_5, model_6)
AIC(model_3, model_4)