#' ---
#' title: "Contextual Variables: Contextual and Compositional effects"
#' author: "MATH 4939"
#' fontsize: 12pt
#' header-includes:
#' - \usepackage{amsmath}
#' - \usepackage{geometry}
#' - \geometry{papersize={9in,6in},left=.4in,right=.4in,top=.4in,bottom=.5in}
#' - \newcommand{\var}{\mathrm{Var}}
#' - \usepackage{mathrsfs}
#' - \raggedright
#' output:
#'   pdf_document:
#'     toc: false
#'     number_sections: false
#'     toc_depth: 3
#'     highlight: tango
#' # bibliography: math4939.bib
#' # link-citations: yes
#' ---
#+ include=FALSE
knitr::opts_chunk$set(fig.asp=.8) 
#' \raggedright
#' 
#' # Estimating within-group and between-group effects with mixed models
#' 
#' Let $\hat{\beta}_X^{MM}$ be the mixed-model estimate of the effect of $X$.    
#' 
#' Then $\hat{\beta}_X^{MM}$ is an unbiased estimate of $\beta_X^{Within}$ only if
#' $\beta_X^{Within} = \beta_X^{Between}$ in the population. This is
#' an assumption that is generally unlikely to be valid.
#' 
#' We will see:
#' 
#' - Using the group-mean^[which can generated using `spida2::cvar(X,group)`]
#'   of $X$: $\bar{X}_j$ 
#'   and the within-group deviation^[which can generated using `spida2::dvar(X,group)`]
#'   of $X$: $X_{ij} - \bar{X}_j$
#'   (also referred by other terms, such as CWG: centered-within-group variable)/   
#'   to estimate to estimate $\beta_X^{Within}$ and $\beta_X^{Between}$.
#' - The 'parallelogram law':   
#'   $$\beta_X^{Between} = \beta_X^{Within} + \beta_X^{Contextual}$$
#' - How the fixed-effects model can use $X_{ij} - \bar{X}_j$ or $X_{ij}$
#'   equivalently with $\bar{X}_j$, but which one you use in the random
#'   effects portion of the model can make a large difference.
#' 
#' Beware:
#' 
#'  - The terms 'contextual' and 'compositional' are not used consistently
#'    by different authors -- or even by the same author.
#'    I like 'between' as a synonym for 'compositional' because it isn't 
#'    ambiguous, but the world prefers words that sound more serious. 
#'  - Adding to the confusion (or creating it) is that the 'effect' of a 
#'    variable (a predictor in a model) depends entirely on the other
#'    variables in the model. That is, the meaning of 
#'    $\frac{\partial}{\partial X}E(Y)$ depends on the other variables because
#'    the choice of other variables determines the path along which predictors
#'    change as $X$ changes. The choice determines the direction of
#'    $\frac{\partial}{\partial X}E(Y)$ viewed as a directional derivative.
#'  - The *effect* of a *contextual variable* for a variable $X$ 
#'    is the *contextual effect* **if** the *raw* value of the variable is
#'    another variable in the model. 
#'  - The *effect* of a *contextual variable*
#'    is the *compositional effect* **if** the *centered within groups* value 
#'    of the variable is another variable in the model. 
#'  - If you change a variable, $X$, in a model (e.g. centre it), you can change 
#'    the meaning of the effects
#'    of all the other variables, unless $X$ and its modified value remain within
#'    the orthogonal complement of the other variables.
#'  - Try to retain the ideas, not just the terminology.  It's easier to apply
#'    the ideas in a new context.
#' \newpage  
#' This is an example using the Public schools in the `hs` data set in the **spida2**
#' package.
#'  
#+ include=FALSE
knitr::opts_chunk$set(comment='    ')
#+
 
library(spida2)
library(nlme)     # this is a version of nlme tailor-made for a specific purpose 
library(lattice)
library(latticeExtra)
library(car)

hsp <- subset(hs, Sector == "Public")   # Public schools from the hs data set

#+ fig1
#|   fig.cap=caption,
#|   fig.pos='h'

caption <- 'Public schools with data ellipses and centroids'

p <- xyplot(mathach ~ ses, hsp, groups = school, pch = 1, cex = .2) +
  layer(panel.ellipse(..., center.pch = 21, center.cex = 1, fill = '#444444'))
p
#' \newpage
#' 
#' Pooled, within and between models:
#' 
fit.pooled <- lm(mathach ~ ses, hsp)
fit.within <- lm(mathach ~ ses + factor(school), hsp) # What happens if we drop 'factor'?
fit.between <- lm(mathach ~ cvar(ses, school), hsp)

fit.list <- list(pooled = fit.pooled, within = fit.within, between = fit.between)
fit.list %>% lapply(coef) %>% lapply(`[`, 2) %>% cbind %>% 
  {colnames(.) <- 'ses' ; .}

hsp <- within(
  hsp,
  {
    pooled <- predict(fit.pooled)
    within <- predict(fit.within)
    between <- predict(fit.between)
  }
)
#+ fig2, fig.cap=caption
caption <-"Within school regression lines using OLS"
p + xyplot(within ~ ses, hsp, groups = school, type = 'l') +
  layer(panel.abline(v = 0))
#+ fig3, fig.cap=caption
caption = "Pooled regression line using OLS"
p + xyplot(within ~ ses, hsp, groups = school, type = 'l') +
  xyplot(pooled ~ ses, hsp, type = 'l', lwd = 3) +
  layer(panel.abline(v = 0))
  
#+ fig4, fig.cap=caption
caption = "Between school regression line using OLS"
p + xyplot(within ~ ses, hsp, groups = school, type = 'l') +
  xyplot(pooled ~ ses, hsp, type = 'l', lwd = 3) +
  xyplot(between ~ cvar(ses,school), hsp, type = 'l', lwd = 4, col = 'red') +
  layer(panel.abline(v = 0))
#'
#' ## What does a mixed model do?
#'  
fit.mm <- lme(mathach ~ ses, hsp, random = ~ 1 | school)
fit.mm
fit.list <- c(fit.list, list(mixed = fit.mm))

(
  fit.list %>% 
  lapply(
    function(fit) {
      if(class(fit)=='lm') coef(fit) else fixef(fit)
    }
  ) %>% 
  sapply(`[`, 2) %>% cbind %>% 
  {colnames(.) <- 'ses' ; .} %>% 
  as.data.frame %>% 
  round(3) %>% 
  sortdf -> coef.table 
)
#'
#' ## Using mixed models to estimate within, between and contextual effects
#'
#' We get equivalent models by using 
#' 
#' - raw ses and the within-group mean of ses, or
#' - centered-within-group ses and the within-group mean of ses
fit.contextual <- lme(mathach ~ ses + cvar(ses, school), hsp, random = ~ 1 | school)
fit.compositional <- lme(mathach ~ dvar(ses, school) + cvar(ses, school), 
                          hsp, random = ~ 1 | school)
wald(fit.contextual)
wald(fit.compositional)
#'
#' Note that variance estimates are the same, although the X matrix of
#' the contextual model is more collinear and, in extreme cases,
#' the fitting algorithm for the mixed model is more likely to
#' produce numerically flawed results.
#'
summary(fit.contextual)
summary(fit.compositional)
#'
#' Easier comparisons:
#' 
wald(fit.contextual)
wald(fit.compositional)
#' 
#' You can get the within and between estimates from both models
#'
 
L.cntx <- rbind(
  'compositional (between) effect' = c(0, 1, 1),
  'within effect' = c(0,1,0),
  'contextual effect' = c(0,0,1)
)
L.comp <- rbind(
  'compositional (between) effect' = c(0, 0, 1),
  'within effect' = c(0,1,0),
  'contextual effect' = c(0,-1,1)
)
L.cntx
L.comp

wald(fit.contextual,    L.cntx)
wald(fit.compositional, L.comp)
#' Comparing with the OLS estimates:
rbind(
  waldf(fit.within)[2, 1:2],
  waldf(fit.between)[2,1:2]
)
#'
#' The contextual model and the compositional model ^[Even models can be
#' contextual or compositional -- from a computer science perspective this
#' is an example of *overloading* or *ad hoc polymorphism*, closely related
#' to the use of *generic functions* and *methods* in R] are **equivalent**
#' since they have X matrices that span the same space. (**Exercise:** Prove it.)
#'
#' Using the contextual model:
#' 
#' $$E(Y) = \psi_0 + \psi_1 X_{ij} + \psi_2 \bar{X}_j$$
#' 
#' With the compositional model:
#' 
#' $$E(Y) = \phi_0 + \phi_1 \left(X_{ij} - \bar{X}_j \right) + \phi_2 \bar{X}_j$$
#' 
#' Rearranging terms:
#' 
#' $$E(Y) = \phi_0 + \phi_1 X_{ij} + (\phi_2 - \phi_1) \bar{X}_j$$
#' 
#' Thus, the two expressions yield identical functions of $X_{ij}$ and
#' $\bar{X}_j$ iff
#' $$\begin{aligned}
#'    \psi_0 &= \phi_0 \\
#'    \psi_1 &= \phi_1 \\
#'    \psi_2 &= \phi_2 - \phi_1
#' \end{aligned}$$ 
#' Equivalently:
#' $$\begin{aligned}
#'    \phi_2 &= \psi_2 + \psi_1
#' \end{aligned}$$ 
#' Note that:
#' $$\begin{aligned}
#'    \beta_0 &= \psi_0 = \phi_0 \\
#'    \beta_{Within} &= \psi_1 = \phi_1 \\
#'    \beta_{Contextual} &= \psi_2 = \phi_2 - \phi_1\\
#'    \beta_{Compositional} &=\beta_{Between} = \phi_2 = \psi_1 + \psi_2 
#' \end{aligned}$$ 
#' 
#' I think of this as the 'parallelogram law'^[I just call it that because I think that giving ideas names that evoke
#' a geometric mnemonic makes them harder to forget.]. Imagine two schools with mean ses
#' one unit apart, for example, 1 and 2. 
#' 
#' Suppose that $\beta_{Within} = 1/2$ and $\beta_{Between} = 3$. 
#' Let $\beta_0 = 1$.
#' 
#' \newpage
#' \phantom{ }  
#' \newpage  
#' 
#' The compositional and the contextual models have equivalent fixed effects models:
#' 
#' $$\begin{aligned}   
#'    {\bf X}_{Contextual} &= \begin{bmatrix} 1 & X_{ij} & \bar{X}_.j\end{bmatrix} \\
#'    {\bf X}_{Compositional} &= \begin{bmatrix} 1 & (X_{ij} -\bar{X}_.j ) & \bar{X}_.j\end{bmatrix} \\
#' \end{aligned}$$ 
#' 
#' **Exercises:** 
#' 
#' - Show that ${\mathscr L}({\bf X}_{Contextual}) = {\mathscr L}({\bf X}_{Contextual})$
#' - Compare the within- and between-group estimates obtained with OLS and with 
#'   mixed models. Compare the SEs of estimates. Are the differences small or large?
#'    
#' ## CWG vs raw variable for random effects
#' 
#' The models we used above did not include a random effect for $X$.  
#' As far as possible mixed models should include random effects for
#' level-1 variables. However, the number of parameters in the $G$
#' matrix grows rapidly with additional random effects and the number
#' of effects that can be included is limited unless one has an
#' enormous amount of data.
#' 
#' Consider the following models where we use the verbose mode to see
#' how the models converge: 
#'  
fit1 <- lme(mathach ~ ses + cvar(ses,school), hsp, random = ~ 1 + ses | school,
            control = list(msVerbose = TRUE))
summary(fit1)
fit2 <- lme(mathach ~ ses + cvar(ses,school), hsp, random = ~ 1 + dvar(ses,school) | school,
            control = list(msVerbose = TRUE))
summary(fit2)
#'
#' There wasn't a large difference in the number of iterations in this case 
#' but with larger more complex models there can be a considerable difference.
#' Sometimes, one or the other model will fail to converge.
#' 
#' Notice that there is a sensible difference in the estimate of the 
#' between effect.
#' 
#' We will see some reasons for this when we take up a deeper dive into the
#' interpretation of the $G$ matrix in the next section.