#' ---
#' title: "Dealing with Heteroskedasticity"
#' subtitle: "R-side Variance Modeling and different G matrices" 
#' geometry: paperheight=4.53in,paperwidth=8in,margin=0.5in
#' author: ''
#' # bibliography: bib.bib
#' output: 
#'   bookdown::pdf_document2
#' header-includes:
#'   - \usepackage{float}
#' ---   
#' \newcommand{\Var}{\mathrm{Var}}
#' \newcommand{\E}{\mathrm{E}}
#' \newcommand{\N}{\mathrm{N}}
#' \pagenumbering{arabic}
#+ include=FALSE
# https://bookdown.org/yihui/rmarkdown/tufte-figures.html
# fig.fullwidth = TRUE, with fig.width and fig.height, out.width
# https://www.statmethods.net/RiA/lattice.pdf
# 
#+ setup2, include=FALSE 
knitr::opts_chunk$set(
  fig.width=7, 
  fig.height=3.75, 
  fig.pos = 'H',
  fig.fullwidth=TRUE,
  #  echo=FALSE,
  #  include=FALSE,
  echo=TRUE,
  include=TRUE,
  comment = '  '
)
#+ setup, include=FALSE 
knitr::opts_chunk$set(comment='  ')
library <- function(...) suppressPackageStartupMessages(base::library(...))
#+
#' 
#' # Setup
#' 
library(nlme)
library(spida2)
library(lattice)
library(latticeExtra)
library(latex2exp)
#'
#' # Introduction
#' 
#' This is an example to illustrate methods to deal with heteroskedasticity. 
#' The methods
#' can be used for 'single-level' and multi-level data.
#' 
#' To model single-level heteroskedasticity or multi-level within-group 
#' heteroskedasticity, 
#' you can use the 'weights' argument to the fitting functions in the 
#' **nlme** package, namely, 'gls', 'nls', 'lme' and 'nlme'.   
#' 
#' To model heteroskedasticity in the between-group model in multi-level models,
#' i.e. heterosckedasticity in the $G$ matrix, you can use the
#' **pdBlocked** function together with other functions from the **pdClasses* 
#' family of functions.
#'  
#' For single-level models, use the 'gls' and 'nls' functions, 
#' and, for mixed models, the 'lme' and 'nlme' functions.
#' 
#' # Generating a data set
#' 
#' Pay equity data set for a hypothetical university with two 
#' faculties: Medicine and Arts with a higher level and 
#' variance in Medicine vs Arts and a different gender gap:
#' 
{
  # Starting data frame:
  dd <- expand.grid(Faculty = c("Arts", "Med"), Sex = c("F","M"), n = 1:400)
  # sample(1000000,1)
  set.seed(576530)
  dd <- within(
    dd,
    {
      Age <- 45 + 5 * (Faculty == "Arts") + 5 * (Sex == "M") + 15 * rnorm(n)
      
      ..esal <- 100 + 20 * (Faculty == 'Med') +
        (4 + .3 *(Sex == "M") + .5 * (Faculty == "Med")) * (Age - 30)
      
      ..sdsal <- 10 + 15 * (Faculty == "Med") + (.2 + .2 * (Faculty == "Med"))* (Age - 30)
      
      Base <- ..esal + ..sdsal * rnorm(n)
      keep <- Age > 28 & Age < 80
      ..sdsal <- NULL
      ..esal <- NULL
    }
  )
}
tab(dd, ~ Faculty + Sex +keep)
dd <- subset(dd, keep)
saveRDS(dd, file = 'salary.rds')
#'
#' # Analysis
#'
dd <- readRDS('salary.rds')
xqplot(dd)
xyplot(Base ~ Age | Faculty, dd, groups = Sex, alpha = .5,
       auto.key = T)

fit <- lm(Base ~ Age * Faculty * Sex, dd)
summary(fit)
for(i in c(1,2,3,5)) {
  plot(fit, which = i, add.smooth=T, mfcol = c(1,1))
}
#'
plot(fit, 1, add.smooth = T)
plot(fit, 2, add.smooth = T)
plot(fit, 3, add.smooth = T)
plot(fit, 5, add.smooth = T)
#'
xyplot(resid(fit) ~ fitted(fit)|Faculty * Sex, dd) %>% 
  useOuterStrips
#' \newpage
#'
#' # Functions in nlme to deal with heteroskedasticity:
#' 
#' Overview:
?varClasses
#' 
#' - varExp: exponential of a covariate or yhat
#' - varPower: power of a covariate or yhat
#' - varConstPower: constant + power of a covariate or yhat
#' - varConstProp: constant + proportion of a covariate or yhat
#' - varIdent: different variance in different subgroups
#' - varFixed: fixed weights given by a covariate
#' - varComb: combination of variance functions
#' - You can also build your own but count on spending a few days
#'   figuring out how to do it
#'   
#' \newpage
#'   
#' Using **nlme::gls** to do the same job as *lm*:
#' 
fit.lm <- lm(Base ~ Age * Faculty * Sex, dd) 
fit <- gls(Base ~ Age * Faculty * Sex, dd)  # same as fit with gls model
summary(fit.lm)
summary(fit)
#'
pred <- with(dd, pred.grid(Faculty, Sex, Age = seq(30,75,1)))
ww <- as.data.frame(wald(fit, pred = pred))
head(ww)
#+ fig.cap=caption
caption <- "Predicted salary for men and women in each Faculty."
xyplot(coef ~ Age | Faculty, ww, groups = Sex,
       type = 'l', auto.key = list(space='right'))
#+ fig.cap=caption
caption <- "Predicted salary for men and women in each Faculty with 95% error bands."
xyplot(coef ~ Age | Faculty, ww, groups = Sex,
       type = 'l', auto.key = list(space='right'),
       fit = ww$coef,
       lower = ww$L2,
       upper = ww$U2,
       subscripts = TRUE) +
  glayer(panel.band(...))
#'
#' \newpage
#' 
#' # Estimating the Gap
#'  
#' Estimating the gap, if it is potentially non-linear in X, requires estimating
#' a many linear combinations of coefficients -- a different linear combination
#' for each value of 'X' where the gap needs to be estimated.
#' 
#' **spida2::wald**, with **spida2::as.data.frame**, allows you to generate 
#' these combinations.
#' 
#' The **L** matrices that produce predicted values for each sex are extracted as
#' follows:
#' 
Lmale <- subset(ww, Sex == 'M')$L
Lfemale <- subset(ww, Sex == 'F')$L
#' 
#' The **L** matrix that produces an estimate of the gap is then:
#' 
Lgap <- Lmale - Lfemale
#'
#' which can then be used with the **wald** function to estimate the gap
#' and its standard error at each value in the 'pred' data frame: 
#'
wgap <- wald(fit, 
             Lgap, 
             data = subset(ww, Sex == 'F', select = c(Faculty,Age)))
#'
#' The 'data' parameter attaches the rows of its argument to the
#' output of the function so the rows can be available for further processing
#' of the result when rendered as a data frame.
#' 
#' Note that the row of the **L** matrix producing each predicted value and
#' values of predictors from the prediction data frame are available in each
#' row together with the predicted values and their standard errors:
#'  
wgap <- as.data.frame(wgap)     
head(wgap)
#'
#' The result can then be plotted as follows. Note that the 
#' **spida2::panel.band** function is designed to use the values
#' passed defined in **xyplot** as `fit`, `lower` and `upper`.
#' It is necessary to add `subscripts = TRUE` so values are
#' passed to `panel.band` in the correct positions. 
#'
#+ fig.cap=caption
caption <- paste('Estimated gaps, estimated without allowing for different', 
                 'variances in different Faculties.')
xyplot(coef ~ Age | Faculty, wgap, 
       type = 'l', auto.key = list(space='right'))
#' 
#+ fig.cap=caption
caption <- 
  paste("Estimated Male-Female wage gap with 95% error bands",
        "estimated without allowing for different variances in different Faculties.")
(
  p <- xyplot(coef ~ Age | Faculty, wgap, 
       type = 'l', auto.key = list(space='right'),
       ylab = 'Wage gap (male - female)',
       ylim = c(-15, 20), xlim = c(30, 75),
       fit = wgap$coef,
       lower = wgap$L2,
       upper = wgap$U2,
       subscripts = TRUE) +
  layer(panel.band(...)) +
  layer_(panel.grid(h=-1,v=-1)) +
  layer_(panel.abline(h=0)) 
)
#' \newpage
#' 
#' # Models with heteroskedasticity
#' 
#' We use the `weights` parameter of `gls` to fit linear models
#' for prediction as well as a number of possible models for
#' variance.  The models use functions in the **nlme** package designed
#' to work with the `weights` parameter.
#' 
#' The following illustrate a number of these functions.
#' 
#' Although these functions are called 'variance functions' 
#' they are in fact 'standard deviation' functions.  
#' 
#' For example
#' 
#' `weights = varConstPower(form = ~fitted(.)|Faculty)`  
#' 
#' 
#' models a multiplicative factor, which when multiplying the estimate
#' of the residual standard error, $\hat{\sigma}$, yields the 
#' estimated standard deviation of epsilons.
#' 
#' The function `varConstPower` fits a model for the standard deviation
#' of epsilons of the form:
#' $$
#' \hat\sigma \left( \theta_1 + |v|^{\theta_2}\right)
#' $$
#' where $\hat{\sigma}$ is the estimate of $\sigma$,
#' $v$ is the 'variance covariate' supplied in a formula
#' as the argument to the `form` parameter, and 
#' $\theta_1$ and $\theta_2$ 
#' are parameters to be estimated.
#' 
#' $v$ can be a variable in 'data' or one of two
#' special variables: **resid(.)** or **fitted(.)**.
#' 
#' If the argument of `form` has a conditioning categorical
#' variable, (following a `|`), then independent estimates
#' of $\theta_1$ and $\theta_2$ are obtained for each level
#' of the conditioning variable.
#' 
#' The following are three heteroskedastic models for base salary.
#' See help files for descriptions of `varPower` and `varIdent`.
#' 
fitconpower <- gls(Base ~ Age * Faculty * Sex, dd,
                   weights = varConstPower(form = ~fitted(.)|Faculty))
summary(fitconpower)
fitpower <- update(fit, weights = varPower(form = ~fitted(.)|Faculty))
summary(fitpower)
fitgroups <- update(fit, weights = varIdent(form = ~ 1 | Faculty))
summary(fitgroups)
#'
#' Like `lme`, `gls` uses REML by default, so that different variance
#' models can be assessed with likelihood methods (LRT, AIC, BIC),
#' provided the linear model is the same. Difference linear models
#' can be assessed with Wald tests of a containing model, or by comparisons
#' if the fitted objects have been refit with `method = 'ML'`.
#'
anova(fit , fitgroups, fitpower, fitconpower)
#'
#' When interpreting the output of `anova` be very careful to interpret
#' p-value only if the models compared are nested with respect to each
#' other. 
#' 
#' Avoid selecting models only on the basis of p-values, AIC or BIC.
#' Don't drop a parameter just because it isn't significant if it plays
#' a role in the interpretation of the model. 
#' 
#' For example don't drop the parameter of a confounding variable needed to
#' block a backdoor path.   Don't drop a parameter if dropping it will
#' affect the correct interpretation of another parameter of key interest.
#' For example a curvature parameter may have a great deal of influence on
#' the estimation of predicted values away from the points of intersection
#' of a straight line and the curved line.
#' 
#' 
#'
library(car)
compareCoefs(fit, fitgroups, fitpower, fitconpower, pvals = TRUE)
#' \newpage
#'
#' ## Revisiting the gap with heteroskedasticity
#'
wgap1 <- wald(fitgroups, 
              Lgap, 
              data = subset(ww, Sex == 'F', select = c(Faculty,Age)))

wgap2 <- wald(fitconpower, 
              Lgap, 
              data = subset(ww, Sex == 'F', select = c(Faculty,Age)))
opt <- options(max.print=30)
#' 
#' Model with different variances in different Faculties
#' 
print(wgap1)
wgap1 <- as.data.frame(wgap1)
head(wgap1)

#' 
#' Model with different variances by age in different Faculties
#' 
print(wgap2)
wgap2 <- as.data.frame(wgap2)
head(wgap2)
options(opt)    # reset to default values
#'
#' Combining results for plotting
#' 
wgap$model <- 'none'
wgap1$model <- 'by Faculty'
wgap2$model <- 'by Age and Faculty'
wg <- rbind(wgap, wgap1, wgap2)

levels <- c('none','by Faculty','by Age and Faculty')
wg$model <- factor(wg$model, levels = levels)

#+ fig.cap=caption
caption <- paste("Estimated Male-Female wage gap with 95% error bands estimated",
                 "with different variance models")

xyplot(coef ~ Age | model * Faculty, wg, 
       type = 'l', auto.key = list(space='right'),
       ylab = 'Wage gap (male - female)',
       ylim = c(-15,20),
       xlim = c(30,75),
       fit = wg$coef,
       lower = wg$L2,
       upper = wg$U2,
       subscripts = TRUE) +
  layer(panel.fit(..., alpha = .2)) +
  layer_(panel.abline(h=0)) +
  layer_(panel.grid(h=-1,v=-1)) -> plhet
useOuterStrips(plhet)
#'
#' __Question:__ Where are the bands wider and where are they narrower when incorporating
#' heteroskedasticity in the model? Do the patterns you see make sense? 
#'
#' \newpage
#' 
#' # G-side modeling and heteroskedasticity
#' 
#' - Different patterns for $G$ matrices
#' - How different subgroups have different $G$ matrices
#' 
#' # Summary
#' 
#' The **nlme** package provides tools to build different forms of $G$ matrices
#' using different 'pd' (for 'positive-definite) functions from the 
#' 'pdClasses' set of functions.
#' 
#' See `?pdClasses`
#' 
#' - `pdSymm`: general positive-definite matrix, without constraints
#' - `pdLogChol`: general positive-definite matrix, without contraints, using
#'   a log-Choleski parametrization, i.e. a Choleski lower-triangular factor
#'   $G = AA'$ with diagonal elements $= \exp(\eta_{ii})$, to ensure that
#'   diagonal elements are positive.
#' - `pdDiag`: diagonal matrix
#' - `pdIdent`: multiple of an identity so different components are uncorrelated
#' - `pdCompSymm`: compound symmetry structure: constant diagonal and constant
#'   off-diagonal elements
#' - `pdBlocked`: block diagonal with diagonal blocks consisting of matrices above
#' 
#' Example using high school data allowing different G matrices for Public
#' and Catholic schools 
#' 
library(spida2)
data(hs)
head(hs)
fit <- lme(mathach ~ ses*Sex*Sector, data = hs,
           random = ~ 1 + ses + Sex | school,
           control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = TRUE))
summary(fit)


sapply(hs, class)
hs$Sex <- factor(hs$Sex)
contrasts(hs$Sex) 
contrasts(hs$Sex) <- rbind(Female = -1/2, Male = 1/2)
contrasts(hs$Sex) 

# What does this do?
# The 0 for 'Sex' is the average of the Male and Female levels.

fit <- lme(mathach ~ ses*Sex*Sector, data = hs,
           random = ~ 1 + ses + Sex | school,
           control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = TRUE))
summary(fit)
getG(fit)

# 
# create 2 indicator variables
# 
hs <- within(
  hs,
  {
    public <- 1*(Sector == 'Public')
    catholic <- 1 - public
  }
)
head(hs)

fit2 <- lme(mathach ~ ses*Sex*Sector, data = hs,
                    random = list(
                      school = pdBlocked(
                        list(
                          pdSymm(~ 0 + public +public:(ses + Sex)),
                          pdSymm(~ 0 + catholic + catholic:(ses + Sex))
                          )
                        )
                      ),
                    control = list(msMaxIter = 1000, msVerbose = TRUE, returnObject = TRUE))
summary(fit2)

getG(fit2)