Week9_LMM2.qmd

---
title: "PSY4210: Linear Mixed Models (LMMs) 2"
author: 
  - name: Pei Hwa Goh
  - name: Joshua F. Wiley
    email: joshua.wiley@monash.edu
  - name: Michelle Byrne
date: "`r Sys.Date()`"
format:
  html:
    toc: true
    number-sections: true
    embed-resources: true
    anchor-sections: true
    smooth-scroll: true
---

These are the `R` packages we will use.

```{r setup}
options(digits = 4)

## new packages are lme4, lmerTest, and multilevelTools

library(data.table)
library(JWileymisc)
library(extraoperators)
library(lme4)
library(lmerTest)
library(multilevelTools)
library(visreg)
library(ggplot2)
library(ggpubr)
library(haven)

## Set your working directory to the folder that has the datafiles.
## "Merged" is a a merged long dataset of both baseline and daily diary
dm <- as.data.table(read_sav("[2021] PSY4210 merged.sav"))

dm[, sex := factor(
  sex,
  levels = c(1,2),
  labels = c("male", "female"))]

## While we are at it, let's also tell R that relationship 
## status is a factor too
dm[, relsta := factor(
  relsta, levels = c(1,2,3), 
  labels = c("single", "in a committed exclusive relationship", "in a committed nonexclusive relationship"))]

dm[, c("Bstress", "Wstress") := meanDeviations(dStress), by = ID]

```

# Linear Mixed Models 

## ML and REML

Two common estimators used with linear mixed mdoels are 
**M**aximum **L**ikelihood (ML) and 
**Re**stricted **M**aximum **L**ikelihood (REML). You can think of
this somewhat the same as the formula for calculating a population's
variance versus estimating the population variance from a 
sample.
Also see: https://stats.stackexchange.com/questions/16008/what-does-unbiasedness-mean/16009#16009.

$$
\sigma^2_{pop} = \frac{\sum_{i=1}^{n} (X_i - \bar{X})^{2}}{n}
$$

$$
\sigma^2_{sample} = \frac{\sum_{i=1}^{n} (X_i - \bar{X})^{2}}{n - 1}
$$

Generally speaking, REML, is less biased so we prefer it as an
estimator and it is the usual default in LMMs. However, for some
comparisons between models to be valid and some statistical tests to
be valid, we need true maximum likelihood (ML) estimates so you will
sometimes see the default REML option "turned off" to get pure ML
estimates or a model be refit using ML instaed of REML. This is done
in `R` using `REML = FALSE` to get ML.

## Random Intercept Model

There are two main uses of intercept only models:

- To calculate the intraclass correlation coefficient (ICC)
- As a comparison to see how much better a more complex model
  fits. 

To calculate the ICC, we use this equation:

$$ICC = \frac{\sigma^{2}_{intercept}}{\sigma^{2}_{intercept} +
\sigma^{2}_{residual}}$$

Following is an example of an intercept only model, where there is
both a fixed effects intercept and a random intercept.
The outcome variable is `stress`.  All predictors come after the
tilde, `~`. In this case, the only "predictors" are the fixed and
random intercept, represented by `1`. The random intercept is random
by `ID`. The function to fit linear mixed models is `lmer()` and
comes from the `lme4` package. It also requires a dataset be
specified, here `dm`. We can get a summary using `summary()`.

```{r}

ri.m <- lmer(dStress ~ 1 + (1 | ID),
            data = dm,
            REML = TRUE)

summary(ri.m)

iccMixed("dStress", id = "ID", data = dm)

``` 

There are four main "blocks" of output from the summary.

1. A repetition of the model options, formula we used, and dataset
   used. This is for records so you know exactly what the model was.
   In *this* model, it shows use that we fit a LMM using restricted
   maximum likelihood (REML) and that the degrees of freedom were
   approximated using Satterthwaite's method. The outcome variable is
   stress (`dStress`) and there are only intercept predictors,
   `1`. The REML criterion at convergence is kind of like the log
   likelihood (LL), but unfortunately cannot be readily used to
   compare across models as easily as the actual LL (e.g., in AIC or
   BIC, which we'll talk about more later).
2. Scaled Pearson residuals. These are raw residuals divided by the
   estimated standard deviation, so that they can be roughly
   interpretted as z-scores. The minimum and maximum are useful for 
   identifying whether there are outliers present in the model
   residuals.
   In *this* model, we can see that the lowest residual is 
   `r min(residuals(ri.m, type = "pearson", scaled=TRUE))` 
   and the maximum residual is 
   `r max(residuals(ri.m, type = "pearson", scaled=TRUE))`
   which are not too large. Large residuals imply that they are 
   extremely unlikely by chance alone and likely represent outliers.
3. Random effects. These show a summary of the random effects in the
   model. Random effects are basically always also fixed effects, so
   the random effects only shows the standard deviation and variance
   of random effects, plus, if applicable, their correlations. The
   means are shown in the fixed effects section. In the case of a
   random intercept only model like this one, there are only two
   random effects: (1) the random intercept and (2) the random
   residual. We have both the standard deviation and variance of
   both. We will use the variances to calculate ICCs.
   In *this* model, the standard deviation of the random intercept, 
   tells us that the average or typical difference between an
   individual's average stress, and the population average
   stress is
   `r as.data.frame(VarCorr(ri.m))[1, "sdcor"]`.
   The standard deviation of the residuals
   tells us that the average or typical difference between an
   individual stress score and the predicted stress
   score is
   `r as.data.frame(VarCorr(ri.m))[2, "sdcor"]`.
   The random effects section also tells us how many observations and
   unique people/groups went into the analysis. 
   In *this* model we can see that we had `r as.integer(ngrps(ri.m))` 
   people providing `r nobs(ri.m)` unique observations.
4. Fixed effects. This section shows the fixed effects. It is a
    table, where each row is for a different effect / predictor and
    each column gives a different piece of information. 
  The "Estimate" is the actual parameter estimate (i.e., THE fixed
    effect, the regression coefficient, etc.). The "Std. Error" is the
    standard error of the estimate, which captures uncertainty in the
    coefficient due to sampling variation. The "df" is the
    Satterthwaite estimated degrees of freedom. As an estimate, it may
    have decimals. The "t value" is the ratio of the coefficient to
    its standard error, that is: $t = \frac{Estimate}{StdError}$. 
	The "Pr(>|t|)" is the p-value, the probability that by chance
    alone one would obtain as or a larger absolute t-value. The
    vertical bars indicate absolute values and the "Pr" stands for
    probability value. Note that `R` uses 
	[scientific E notation](https://en.wikipedia.org/wiki/Scientific_notation).
	The number following the "e" indicates how many places to the
    right (if positive) or left (if negative) the decimal point should
    be moved. For example, 0.001 could be written 1e-3. 0.00052 could
    be written 5.2e-4. These often are used for p-values which may be
    numbers very close to zero.
	In *this* model, we can see that the fixed effect for the
    intercept is `r fixef(ri.m)[["(Intercept)"]]` which is like the
    the mean of the random intercept and tells us the average
    level of stress, in this instance since there are no
    other predictors in the model.

Profile likelihood confidence intervals can be obtained using the 
`confint()` function. These confidence intervals capture the
uncertainty in parameter estimates for both the fixed and random
effects due to sampling variation. They do not capture individual
differences directly. Note that you only get confidence intervals for
random effects when using the profile method, not when
`method = "Wald"` although the Wald method is much faster.
Profile likelihood confidence intervals are the default and are
generally more accurate than Wald based confidence intervals, but also
can be slower to calculate. Here both ways are shown.

```{r}

## Profile Confidence Intervals
ri.ci <- confint(ri.m, method = "profile", oldNames = FALSE)
ri.ci

## Wald Confidence Intervals
confint(ri.m, method = "Wald", oldNames = FALSE)

```

### Diagnostics

Typical diagnostics and checks include checking for outliers,
assessing whether the distributional assumptions are met, checking for
homogeneity of variance and checking whether there is a linear
association between predictors and outcome. With only an intercept,
there is no need for checking whether a linear association is
appropriate.

As for linear regressions, we can use `modelDiagnostics()` to
calculate diagnostics and `plot()` to get a plot of the diagnostics.

The top left plot of the residuals helps check for outliers on the
residuals.
These plots show one relatively extreme value on the residuals. In
this case, using the scaled pearson residuals, which are roughly like
z scores, the size of the residual outliers are unlikely to be an 
issue (-3.13 is not that large in z scores). The bottom left plot 
shows the distribution of the random intercept. No outliers are 
observed.

The density plots (and QQ deviates plot for residuals) indicate 
mild non-normality, but it is not too extreme and close enough for
inference.

Finally, on the top right, we check the homogeneity of variance. 
There is no clear trend in the residuals indicating that the
homogeneity of variance assumption is reasonably met.
The residuals show a characteristic banding when dealing with
"continuous" variable that have a handful of possible values.
This is not a problem per se, although if too extreme may indicate
that treating the data as continuous is not a great idea.

```{r}

## When I initially tried to run the model diagnostics,
## an error about haven_labelled class appeared, so let's
## get rid of the haven_labelled class type for stress
## before running the model and model diagnostics
dm[, dStress := as.numeric(dStress)]

## Now we re-run the model
ri.m <- lmer(dStress ~ 1 + (1 | ID),
            data = dm,
            REML = TRUE)

## And re-run the model diagnostics
plot(modelDiagnostics(ri.m), ncol = 2, nrow = 2, ask = FALSE)

``` 

With diagnostics reasonably met, we proceed with a write up.

### Sample Write Up

An intercept only linear mixed model was fit to 
`r nobs(ri.m)` stress scores from 
`r as.integer(ngrps(ri.m))` people. The intraclass correlation
coefficient was 
`r as.data.frame(VarCorr(ri.m))[1, "vcov"] / sum(as.data.frame(VarCorr(ri.m))[, "vcov"])` 
indicating that about 40% of the total variance in stress
was between people and the other 60% is within person due to
fluctuations across days. The fixed effect intercept revealed that the
average [95% CI] stress was 
`r fixef(ri.m)[["(Intercept)"]]`
`r sprintf("[%0.2f, %0.2f]", ri.ci[3, 1], ri.ci[3, 2])`.
However, there were individual differences, with the standard
deviation for the random intercept being
`r as.data.frame(VarCorr(ri.m))[1, "sdcor"]`
indicating that there are individual differences in the mean
stress. Assuming the random intercepts follow a normal distribution,
we expect most people to fall within one standard deviation of the
mean, which in these data would be somewhere between:
`r fixef(ri.m)[["(Intercept)"]] + c(-1, 1) *  as.data.frame(VarCorr(ri.m))[1, "sdcor"]`. 


## Random Intercept and Fixed Effects Models

Following is an example of a LMM with fixed effects and a random
intercept (no random slopes). Although we did not explicitly add a
fixed effects intercept by adding `1` to the equation, it is there by
default. We still have a random intercept.

```{r}

fp.m <- lmer(dStress ~ dEnergy + (1 | ID),
            data = dm,
            REML = TRUE)

summary(fp.m)

``` 

There are four main "blocks" of output from the summary.

1. A repetition of the model options, formula we used, and dataset
   used. This is for records so you know exactly what the model was.
   In *this* model, it shows use that we fit a LMM using restricted
   maximum likelihood (REML) and that the degrees of freedom were
   approximated using Satterthwaite's method. The outcome variable is
   stress (`stress`) and energy is a predictor.
   The REML criterion at convergence is kind of like the log
   likelihood (LL), but unfortunately cannot be readily used to
   compare across models as easily as the actual LL (e.g., in AIC or
   BIC).
2. Scaled Pearson residuals. These are raw residuals divided by the
   estimated standard deviation, so that they can be roughly
   interpretted as z-scores. The minimum and maximum are useful for 
   identifying whether there are outliers present in the model
   residuals.
   In *this* model, we can see that the lowest residual is 
   `r min(residuals(fp.m, type = "pearson", scaled=TRUE))` 
   and the maximum residual is 
   `r max(residuals(fp.m, type = "pearson", scaled=TRUE))`
   which are acceptable. Absolute residuals of 10, for example, would be
   large enough that they are extremely unlikely by chance alone and likely
   represent outliers. We can see there are some more extreme negative
   than positive residuals. That means that predictions are sometimes
   too (extremely) high rather than too (extremely) low.
3. Random effects. These show a summary of the random effects in the
   model. Random effects are basically always also fixed effects, so
   the random effects only shows the standard deviation and variance
   of random effects, plus, if applicable, their correlations. The
   means are shown in the fixed effects section. In the case of a
   model where the only random effect is the intercept, the
   random effects show: (1) the random intercept and (2) the random
   residual. We have both the standard deviation and variance of
   both. 
   In *this* model, the standard deviation of the random intercept, 
   tells us that the average or typical difference between an
   individual's estimated stress when energy is 0, 
   and the population average estimated stress when energy is
   0 is
   `r as.data.frame(VarCorr(fp.m))[1, "sdcor"]`.
   The standard deviation of the residuals
   tells us that the average or typical difference between an
   individual stress score and the predicted stress
   score is
   `r as.data.frame(VarCorr(fp.m))[2, "sdcor"]`.
   The random effects section also tells us how many observations and
   unique people/groups went into the analysis. 
   In *this* model we can see that we had `r as.integer(ngrps(fp.m))` 
   people providing `r nobs(fp.m)` unique observations.
4.  Fixed effects. This section shows the fixed effects. It is a
    table, where each row is for a different effect / predictor and
    each column gives a different piece of information.
	The "Estimate" is the actual parameter estimate (i.e., THE fixed
    effect, the regression coefficient, etc.). The "Std. Error" is the
    standard error of the estimate, which captures uncertainty in the
    coefficient due to sampling variation. The "df" is the
    Satterthwaite estimated degrees of freedom. As an estimate, it may
    have decimals. The "t value" is the ratio of the coefficient to
    its standard error, that is: $t = \frac{Estimate}{StdError}$. 
	The "Pr(>|t|)" is the p-value, the probability that by chance
    alone one would obtain as or a larger absolute t-value. The
    vertical bars indicate absolute values and the "Pr" stands for
    probability value. Note that `R` uses 
	[scientific E notation](https://en.wikipedia.org/wiki/Scientific_notation).
	The number following the "e" indicates how many places to the
    right (if positive) or left (if negative) the decimal point should
    be moved. For example, 0.001 could be written 1e-3. 0.00052 could
    be written 5.2e-4. These often are used for p-values which may be
    numbers very close to zero.
	In *this* model, we can see that the fixed effect for the
    intercept is `r fixef(fp.m)[["(Intercept)"]]` which is like
    the mean of the random intercept and tells us the average
    estimated stress score when energy = 0.
	The fixed effect (regression coefficient) for energy is 
	`r fixef(fp.m)[["dEnergy"]]` which tells us how much on average
    (fixed effect) lower stress is expected to be when energy
    is one unit higher. 

Profile likelihood confidence intervals can be obtained using the 
`confint()` function. These confidence intervals capture the
uncertainty in parameter estimates for both the fixed and random
effects due to sampling variation. They do not capture individual
differences directly. Note that you only get confidence intervals for
random effects when using the profile method, not when
`method = "Wald"` although the Wald method is much faster.

```{r}

fp.ci <- confint(fp.m, method = "profile", oldNames = FALSE)
fp.ci

```

### Diagnostics and Checks

Typical diagnostics and checks include checking for outliers,
assessing whether the distributional assumptions are met, checking for
homogeneity of variance and checking whether there is a linear
association between predictors and outcome. 

Results look about the same as for the intercept only model.
 
```{r}

plot(modelDiagnostics(fp.m),
     ncol = 2, nrow = 2, ask = FALSE)

```

### Sample Write Up

To examine the association of stress and energy, a linear
mixed model was fit. The final model included `r nobs(fp.m)` stress
scores from `r as.integer(ngrps(fp.m))` people.  
The fixed effect intercept revealed that the
average [95% CI] stress when energy is 0 was 
`r fixef(fp.m)[["(Intercept)"]]`
`r sprintf("[%0.2f, %0.2f]", fp.ci[3, 1], fp.ci[3, 2])`.
However, there were individual differences, with the standard
deviation for the random intercept being
`r as.data.frame(VarCorr(fp.m))[1, "sdcor"]`
indicating that there are individual differences in the mean
stress. Assuming the random intercepts follow a normal distribution, 
we expect most people to fall within one standard deviation of the
mean, which in these data would be somewhere between:
`r fixef(fp.m)[["(Intercept)"]] + c(-1, 1) *  as.data.frame(VarCorr(fp.m))[1, "sdcor"]`. 
Using Satterthwaite's approximation for degrees of freedom revealed
that energy was statistically significantly associated with stress 
(p = .005). On average across people, a one unit higher energy score
was associated with 0.16 lower stress scores.


## Between and Within Effects

When we have a time-varying predictor variable, we can separate it
into a between and within portion and both of these new variables can
be included in a LMM as predictors. After creating a between and
within version of energy, `Benergy` and `Wenergy` we just enter both
as fixed effects predictors.

```{r}

dm[, c("Benergy", "Wenergy") := meanDeviations(dEnergy), by = ID]

fp.m2 <- lmer(dStress ~ Benergy + Wenergy + (1 | ID),
              data = dm)

summary(fp.m2)

``` 

There are four main "blocks" of output from the summary.

1. A repetition of the model options, formula we used, and dataset
   used. This is for records so you know exactly what the model was.
   In *this* model, it shows use that we fit a LMM using restricted
   maximum likelihood (REML) and that the degrees of freedom were
   approximated using Satterthwaite's method. The outcome variable is
   stress (`dStress`) and `Benergy` and `Wenergy` are fixed effects
   predictors. 
   The REML criterion at convergence is kind of like the log
   likelihood (LL), but unfortunately cannot be readily used to
   compare across models as easily as the actual LL (e.g., in AIC or
   BIC).
2. Scaled Pearson residuals. These are raw residuals divided by the
   estimated standard deviation, so that they can be roughly
   interpretted as z-scores. The minimum and maximum are useful for 
   identifying whether there are outliers present in the model
   residuals.
   In *this* model, we can see that the lowest residual is 
   `r min(residuals(fp.m2, type = "pearson", scaled=TRUE))` 
   and the maximum residual is 
   `r max(residuals(fp.m2, type = "pearson", scaled=TRUE))`
   which while a bit large, 
   are not so extreme if interpretted as z-scores as to be
   concerning. Absolute residuals of 10, for example, would be large enough
   that they are extremely unlikely by chance alone and likely
   represent outliers. We can see there are slightly higher negative
   than positive residuals. That means that predictions are sometimes
   too (extremely) high rather than too (extremely) low relative to
   actual values.
3. Random effects. These show a summary of the random effects in the
   model. Random effects are basically always also fixed effects, so
   the random effects only shows the standard deviation and variance
   of random effects, plus, if applicable, their correlations. The
   means are showon in the fixed effects section. In the case of a
   model where the only random effect is the intercept, the
   random effects show: (1) the random intercept and (2) the random
   residual. We have both the standard deviation and variance of
   both. 
   In *this* model, the standard deviation of the random intercept, 
   tells us that the average or typical difference between an
   individual's estimated stress when `Benergy` and `Wenergy` are 0, 
   and the population average estimated stress when the predictors are
   0 is
   `r as.data.frame(VarCorr(fp.m2))[1, "sdcor"]`.
   The standard deviation of the residuals
   tells us that the average or typical difference between an
   individual stress score and the predicted stress
   score is
   `r as.data.frame(VarCorr(fp.m2))[2, "sdcor"]`.
   The random effects section also tells us how many observations and
   unique people/groups went into the analysis. 
   In *this* model we can see that we had `r as.integer(ngrps(fp.m2))` 
   people providing `r nobs(fp.m2)` unique observations.
4.  Fixed effects. This section shows the fixed effects. It is a
    table, where each row is for a different effect / predictor and
    each column gives a different piece of information.
	The "Estimate" is the actual parameter estimate (i.e., THE fixed
    effect, the regression coefficient, etc.). The "Std. Error" is the
    standard error of the estimate, which captures uncertainty in the
    coefficient due to sampling variation. The "df" is the
    Satterthwaite estimated degrees of freedom. As an estimate, it may
    have decimals. The "t value" is the ratio of the coefficient to
    its standard error, that is: $t = \frac{Estimate}{StdError}$. 
	The "Pr(>|t|)" is the p-value, the probability that by chance
    alone one would obtain as or a larger absolute t-value. The
    vertical bars indicate absolute values and the "Pr" stands for
    probability value. Note that `R` uses 
	[scientific E notation](https://en.wikipedia.org/wiki/Scientific_notation).
	The number following the "e" indicates how many places to the
    right (if positive) or left (if negative) the decimal point should
    be moved. For example, 0.001 could be written 1e-3. 0.00052 could
    be written 5.2e-4. These often are used for p-values which may be
    numbers very close to zero.
	In *this* model, we can see that the fixed effect for the
    intercept is `r fixef(fp.m2)[["(Intercept)"]]` which is like
    the mean of the random intercept and tells us the average
    estimated stress score when both the between and within energy
    = 0. The fixed effect (regression coefficient) for `Benergy` is 
	`r fixef(fp.m2)[["Benergy"]]` which tells us how different on average
    (fixed effect) stress is expected to be when average energy
    is one unit higher. That is, for people who, in general (on
    average) have higher energy, how much lower stress do they have on
    average across days?
	The fixed effect (regression coefficient) for `Wenergy` is 
	`r fixef(fp.m2)[["Wenergy"]]` which tells us how much on average
    (fixed effect) stress is expected to change when energy
    is one unit higher than an individuals' own average. 
	That is, on days when someone has one more unit higher energy than
    usual (than their own mean), how much lower stress do they have on
    that same day?
	The approximate p-values indicate that the average energy scores are
	  statistically significant, indicating that people who have more 
	  energy in general have significantly lower stress in general. 
	  The daily energy scores are NOT significant, but if they were, would suggest
	  that within people, higher energy days would also be associated with
	  lower stress days. 

Profile likelihood confidence intervals can be obtained using the 
`confint()` function. These confidence intervals capture the
uncertainty in parameter estimates for both the fixed and random
effects due to sampling variation. They do not capture indivdiual
differences directly. Note that you only get confidence intervals for
random effects when using the profile method, not when
`method = "Wald"` although the Wald method is much faster.

```{r}

fp.ci2 <- confint(fp.m2, method = "profile", oldNames = FALSE)
fp.ci2

```
### Diagnostics and Checks

Typical diagnostics and checks include checking for outliers,
assessing whether the distributional assumptions are met, checking for
homogeneity of variance and checking whether there is a linear
association between predictors and outcome. 

Results look about the same as for the intercept only model.
Homogeneity of variance still appears approximately met.
 
```{r}

plot(modelDiagnostics(fp.m2),
     ncol = 2, nrow = 2, ask = FALSE)

```

### Sample Write Up

To examine the association of stress and energy, a linear
mixed model was fit. To understand the association of energy and
stress at both the between person and within person level, daily
energy ratings were separated into average energy across the five days
of the study (between energy) and daily deviations of energy from
individuals' own mean energy level (within energy).
The final model included `r nobs(fp.m2)` stress
scores from `r as.integer(ngrps(fp.m2))` people.  
The fixed effect intercept revealed that the
average [95% CI] stress when between and within energy are 0 was 
`r fixef(fp.m2)[["(Intercept)"]]`
`r sprintf("[%0.2f, %0.2f]", fp.ci2[3, 1], fp.ci2[3, 2])`.
However, there were individual differences, with the standard
deviation for the random intercept being
`r as.data.frame(VarCorr(fp.m2))[1, "sdcor"]`
indicating that there are individual differences in the mean
stress. Assuming the random intercepts follow a normal distribution, 
we expect most people to fall within one standard deviation of the
mean, which in these data would be somewhere between:
`r fixef(fp.m2)[["(Intercept)"]] + c(-1, 1) *  as.data.frame(VarCorr(fp.m2))[1, "sdcor"]`. 
Using Satterthwaite's approximation for degrees of freedom revealed
that only between person energy was statistically significantly 
associated with stress. On average across people, at the between person level, 
a one unit higher average energy score was associated with 0.30 
lower stress scores. Although not significant, on average across people, 
at the within person level, a one unit higher than usual daily energy score was 
associated with 0.12 lower stress scores that same day. The magnitude 
suggests that a one unit higher average energy level has about 
twice the magnitude of association with average stress as does 
a one unit higher daily energy level.

## Statistical Inference 
An extra note about how we got the df's above.
There is ambiguity in terms of how best to calculate degrees of
freedom (df) for LMMs. By default `R` does not calculate the df and so
does not provide p-values for the regression coefficients (fixed
effects) from LMMs.

One easy, albeit imperfect, solution is to use the `lmerTest`
package. `lmerTest` use Satterthwaite's method to calculate
approximate degrees of freedom and use these for the t-tests and
p-values for each regression coefficient. To use `lmerTest` simply
make sure that **both** `lme4` and `lmerTest` packages are installed
and that you load the `lmerTest` package after `lme4`, by using:
`library(lmerTest)`. This is shown in the examples above.
Once that is done, all regular calls to `lmer()` function used to fit
LMMs will automatically have df estimated and p-values.

Although this is a relatively common approach, it is not without
limitations or debate. Other approaches include:

- assuming the sample size is large enough that the degrees of freedom
  are large enough that we can assume the parameters follow a normal
  distribution instead of a t-distribution.
- Relying on things like profile-based confidence intervals
- Bootstrapping, which is time intensive but probably the most robust
  of these options.

Here is a short example bootstrapping. Here we bootstrap the 95%
confidence intervals, and if the interval does not include 0, we can
conclude that $p < .05$, it is statistically significant.
Although bootstrapping is robust, we often do not use it as it can be
quite time intensive for more complex models, especially those with
more random effects.

```{r}

confint(fp.m, method = "boot")

``` 

# Random Slopes

## Theory / Conceptual

Beyond random intercepts, we can allow 
allow regression slopes to differ by person. In LMMs, we always
include a random intercept. That means if we also include a random
slope, we will have two random effects. With two (or more) random
effects, because they are assumed to be a random variable with a
normal distribution, we can also look at how the random effects
correlate with each other. That is, just as how you could look at how
two random variables, say age and stress, correlate with each other,
you can look at how any two random effects correlate with each other
(e.g., how a random intercept and random slope correlate).

When working with random slopes, we also need to slightly modify our
understanding of the assumptions and notation for LMMs.

Let $\mathbf{G}$ be a square, $q$ x $q$ covariance matrix, where $q$
is the number of random effects in the model. In our simplest of
examples, there is only one random effect, the random intercept, so $q
= 1$ and the random effect covariance matrix is:

$$
\mathbf{G} =
\begin{bmatrix}
\sigma^{2}_{int} \\
\end{bmatrix}
$$

If we had a random intercept only, then we assume:

$$ u_{0j} \sim \mathcal{N}(0, \mathbf{G}) $$

This is the same assumption we covered in the introduction to LMMs,
however the subtle change in notation sets us up for more complexity.
By using $\mathbf{G}$ we are replacing a single standard
deviation/variance with a covariance matrix. The benefit of the matrix
is that a covariance matrix can be small, like a 1 x 1 matrix or
bigger like a 2 x 2 or $q$ x $q$ matrix for any number of random
effects.

So what is the implication of this change? We do not *just* assume
that each random effect follows a normal distribution. In fact, for
LMMs, the assumption is that the random effects follow a
**multivariate** normal distribution^[for more info, see here:
https://en.wikipedia.org/wiki/Multivariate_normal_distribution]. 
The multivariate normal distribution (MVN) is a distribution for
multiple variables. If $q$ variables follow a MVN, then each $q_i$
variable will itself follow a univariate normal distribution.
However, just because each $q_i$ variable individually follows a
univariate normal distribution does not mean that the $q$ variables
together follow a MVN.  In other words, MVN implies univariate normal,
but univariate normal does not imply MVN.

Like the univariate normal distribution, the MVN has two basic
parameters, the mean and variance. However, unlike a univariate normal
distribution, the MVN the means will be a vector of means, one for
each variable, and the variances, will be a $q$ x $q$ covariance
matrix. 

Back to LMMs, what this means is that if we have a random intercept
and a random slope, $q = 2$ and we'd have:

$$
\mathbf{G} =
\begin{bmatrix}
\sigma^{2}_{int} & \sigma_{int,slope} \\
\sigma_{int,slope} & \sigma^{2}_{slope} \\
\end{bmatrix}
$$

The variances are on the diagonal and the covariance (the
unstandardized correlation between the intercept and the slope) is on
the off diagonal.

The usual practice in LMMs is to freely estimate both the variances
and covariances (correlations) for any random effects. We can,
however, also assume that the random effects are uncorrelated, that
is, assume that the random effects follow a multivariate normal
distribution like this:

$$
\mathbf{G} =
\begin{bmatrix}
\sigma^{2}_{int} & 0 \\
0 & \sigma^{2}_{slope} \\
\end{bmatrix}
$$

Although random slopes are different in some respects from a random
intercept, in most ways they are comparable. In both cases we assume a
normal distribution, with a mean (the "fixed effect" part, the average
slope across people) and a standard deviation (how much variability
there is in the slope across people). We estimate the mean and
standard deviation, but we *can* get BLUPs for the individual slopes
as well, which are predictions about what the slope is for any
specific individual. While a random intercept only impacts the level
of lines, a random slope will impact both the level and the slope of
the lines, shown graphically in the following figure.

```{r, fig.width = 7, fig.height = 10, fig.cap = "Figure showing a random intercept or random intercept and slope"}

ggarrange(
ggplot(expand.grid(x = c(0, 10), y = c(0, 10)), aes(x, y)) +
  geom_point(colour = "white") + 
  geom_abline(intercept = 1, slope = 1, colour = "black", linewidth = 1) +
  geom_abline(intercept = 2, slope = 1, colour = "yellow", linewidth = 1) +
  geom_abline(intercept = 3, slope = 1, colour = "blue", linewidth = 1) +
  geom_abline(intercept = 4, slope = 1, colour = "purple", linewidth = 1) +
  theme_pubr() +
  coord_cartesian(xlim = c(0, 9), ylim = c(0, 10), expand=FALSE) +
  ggtitle("Random Intercepts"),
ggplot(expand.grid(x = c(0, 10), y = c(0, 10)), aes(x, y)) +
  geom_point(colour = "white") + 
  geom_abline(intercept = 1, slope = 1, colour = "black", linewidth = 1) +
  geom_abline(intercept = 2, slope = .5, colour = "yellow", linewidth = 1) +
  geom_abline(intercept = 3, slope = 1.5, colour = "blue", linewidth = 1) +
  geom_abline(intercept = 4, slope = .5, colour = "purple", linewidth = 1) +
  theme_pubr() +
  coord_cartesian(xlim = c(0, 9), ylim = c(0, 10), expand = FALSE) +
ggtitle("Random Slopes"),
ncol = 1, nrow = 2)

```

To evaluate the distributional assumptions of LMMs with multiple
random effects, the ideal test is to evaluate whether the random
effects come from a MVN. Visualizing MVN distributions is not easy,
especially with more than 2 dimensions, so instead we use another
approach.

The Mahalanobis distance measures the distance between a point and a
space defined by a vector of means and a covariance matrix. The
"point" can be unidimensional or multidimensional and if multidimensional
is not limited to only two dimensions. That is, the Mahalanobis
distance can compute the distance between a multidimensional point and
a multidimensional distribution. This is very convenient as it scales
arbitrarily to however complex (many random effects) or simple (one
random effect) we may have. If for each row of data, we calculate its
Mahalanobis distance, those squared distances will follow a $\chi^{2}$
(chi-squared) distribution with degrees of freedom equal to the number
of dimensions ($p$), that is a $\chi^2$ distribution with $df =
p$ if the raw data were MVN. This allows us to compare our observed
Mahalanobis distances to a chi-squared distribution and if they are a
close match, it indicates the data on which the Mahalanobis distances
were calculated were MVN. Thus, we can use Mahalanobis distances to
evaluate whether a set of variables have a MVN distribution.
They also can help identify multivariate outliers. We will see
examples of this when we look at diagnostics for LMMs with multiple
random effects.

## Random Slopes in `R`

Random slopes can only be estimated for predictors that vary within
units. Using the merged daily data collection dataset, we can
look at energy predicting stress. First, we will take a look at a model
with a random intercept and `dEnergy` as a fixed effect only.

```{r}

## get rid of the haven_labeled class type for stress
dm[, dStress := as.numeric(dStress)]

m0 <- lmer(dStress ~ dEnergy + (1 | ID),
           data = dm)

summary(m0)

plot(modelDiagnostics(m0), ncol = 2, nrow = 2, ask = FALSE)

```

From the results, we can see that on average, higher energy is
associated with lower stress. Energy has not been separated into a
between or within component, so we cannot say whether its driven by
people who have higher energy on average having lower stress or days
with higher than usual energy having lower than usual stress, or both.
We can see that the residuals (top left graph panel) and random
intercept by ID (bottom left graph panel) are about normally
distributed with only modest if any outliers and we can see that the
homogeneity of variance assumption is approximately met at least in
the residuals vs fitted/predicted values (top right).

Next, we add dEnergy as a random slope by adding it inside the
parentheses before ID, `(dEnergy | ID)`.

```{r}

m1 <- lmer(dStress ~ dEnergy + (dEnergy | ID),
           data = dm)

summary(m1)

```

The output is fairly familiar to the fixed effect slope and random
intercept only model but there is another random effect under the
groups ID for `dEnergy`. We get the variance and standard deviation of
`dEnergy` as well as the correlation between the random intercept and the
random slope. The correlation is 
`r as.data.table(VarCorr(m1))[3, sdcor]` indicating that people who
have a higher level of stress when energy is 0 (the intercept) also tend
to have a more negative slope of the association between stress and energy. Note that although
the fixed effect slope estimate for energy is about the same between the
two models, the standard error is larger so the p-value is higher
(further away from 0) for the model with a fixed and random slope of
`dEnergy` than the model with only a fixed slope of `dEnergy`. This is a
fairly common pattern.

Now we can look at model diagnostic plots. We use 2 columns and 3 rows
as there are more plots now, but otherwise we use the usual 
`modelDiagnostics()` function.


```{r, fig.width = 7, fig.height = 10}

plot(modelDiagnostics(m1), ncol = 2, nrow = 3, ask = FALSE)

```

In the figure, we see our familiar density and QQ deviates plot for
the residuals (top left). Univariate density plots (black lines) and
univariate normal distributions (blue dashed lines) for the random
intercept alone (`ID : (Intercept)`) and the random slope of energy
alone (`ID : dEnergy`). The naming convention is to include the
coefficient name, intercept or energy for the energy slope, prefixed by
what the variable name is for the random effect, here ID. Both of the
univariate random effect distribution plots can be interpreted as
usual and are shown on the middle row, left and right panels. What is
technically graphed are the BLUPs and at the bottom we can get a five
number summary to help see the minimum, 25th percentile, 50th
percentile (median), 75th percentile, and maximum, which can give us
some sense of the spread of the random intercept and slope. For
example, for the random slope of energy, the maximum BLUP is 0.07,
indicating that the highest predicted individual slope for anyone is
0.07. This tells us that not many people are expected to have *higher* stress if
they have higher energy. We can also see that 50% of people have BLUPs
between -0.20 and -0.07, giving us a sense of the spread of the random
slopes.

The new graph is on the bottom row, left side and it evaluates whether
the random effects by ID follow a MVN or not by using the Mahalanobis
distances. The observed density is the black line. The theoretical
density again is in dashed blue line. Here the theoretical density
**is not** a normal density, but a chi-squared density with $df = p =
2$. In this case, we can see that its not a terrible fit between the
observed and theoretical chi-squared density, suggesting the MVN
assumption is reasonably well met. We also can see, however, that
there are some relatively extreme distances, with the maximum at 13.72
being quite a bit higher than the next nearest. Although its not
flagged as an extreme value, partly because there are only about 88
people.

In practice, I may not actually make any changes here. For the sake of
illustration, however, and to show how to do it, we will use a more
stringent criteria for extreme values. Rather than defining an extreme
value as the top and bottom 0.5% of the theoretical distribution,
let's use the top and bottom 1% of the theoretical distribution.
Because I know this will yield some extreme values to remove, I am
saving the results of `modelDiagnostics()` in an object, `m1.diag`
which I can plot but I can also use to identify which rows / IDs in
the dataset are extreme and should / could be dropped.
The new figure shows a few extreme values in the random effect
BLUPs. Maybe a low value on the intercept, definitely a couple too low
slopes, and maybe one MVN extreme value, the 13.72 Mahalanobis
distance.

```{r, fig.width = 7, fig.height = 10}

m1.diag <- modelDiagnostics(m1, ev.perc = .01)
plot(m1.diag, ncol = 2, nrow = 3, ask=FALSE)

``` 

To consider removing these cases, we need to identify them.
The `m1.diag` object in `R` has several subparts, but we are
particularly interested here in the `extremeValues` subpart, which is
itself a little data table, shown in the following.

```{r}

m1.diag$extremeValues

``` 

In this extreme values data table, we can see a few columns. The most
important columns are:

- `ID` this is the ID variable from the dataset and will be helpful
  later.
- `Index` this is the row number in the dataset used in the LMM and
  can be used to identify specific rows of data that are extreme.
- `EffectType` this indicates what type of effect the extreme value
  was identified on.
  
The first column, `dstress` just shows the outcome variable score for
reference, which may help but is not necessarily that useful always.
The name of the first column will depend on the name of the outcome
variable.

In this dataset, we can see that there are three extreme values
identified on the residuals. For all the random effects, because the
BLUPs are calculated per person, not per observation if a person is
classified as an extreme value, then all observations belonging to
that specific ID will be shown. That is because at the random effect
level, a person as a whole unit is either extreme or not and if
extreme the only "choice" would be to exclude / modify that entire person.

Because there are multiple types of extreme values, we could have some
choice in how to address them. We could remove any extreme values, or
remove one at a time and re-run the model to see if anything else
remained extreme or not. For example, we can see that ID 23 is an
extreme value on the random effect of dEnergy as well as the random
intercept. Dropping ID 23 might change the rest enough that say some
of the other residuals are no longer extreme. Conversely, we could
drop some of the extreme residuals, specific 'weird' days and see if
that happens to address any of the random effects. In this case, it is
not too likely as the extreme residuals come from different
participants (IDs 29, 64, 110) than do the extreme values on the random
effects.

A relatively intense approach would be to drop all of these extreme
values and re-estimate the LMM in the dataset with these rows / IDs
excluded.  Here are different ways that we could do that.

The first approach uses the `Index` column from the extreme values
data table and then we pass that to our dataset, `dm` and use the
minus sign to indicate we want to drop those rows of data.
Note that because the same rows are extreme on a few different
measures, we use the `unique()` function so that we only drop each row
once. This never hurts to use, as if all rows are already unique its
fine, but it helps if there are duplicate extreme values (e.g., Index
29 is extreme on the random intercept and random effect of dEnergy).

```{r}

dm.noev1 <- dm[-unique(m1.diag$extremeValues$Index)]

``` 

Another approach, supposing we only wanted to drop selected IDs. For
example, we could decide that we want to only get rid of ID 23. 
Here we use the `%nin%` operator which takes a variable
on the left and returns `TRUE` if it is not in (nin) the value/vector
of values on the right hand side. This operator comes from the
`extraoperators` package so we need to have that loaded.

```{r}

dm.noev2 <- dm[ID %nin% 23]

#Here I'm going to create an alternate version of this with a different ID (29) to show you something later, forget it for now:
dm.noev2a <- dm[ID %nin% 29]

``` 

Supposing we wanted to remove IDs 23, 109 and 143, we could use the
same approach but instead of listing one ID, we'd list all three as a
vector, like this:

```{r}

dm.noev3 <- dm[ID %nin% c(23, 109, 143)]

``` 

If we wanted to only exclude one type of effect, say only the MVN
extreme values, we could use the row indices but pre-select specific
effect types. We can do this because the extreme values dataset is a
regular data table, so we can operate on it the same as we would on
any dataset. Let's just confirm that we can subset the extreme values
dataset to only give us the MVN extreme values.


```{r}

m1.diag$extremeValues[EffectType == "Multivariate Random Effect ID"]

```

Now that that works, we can use the same approach we did to exclude
all extreme values, but using this subsetted dataset. We replace:
`m1.diag$extremeValues` with 
`m1.diag$extremeValues[EffectType == "Multivariate Random Effect ID"]`
and then we use just the `Index` column as before by writing:
`$Index`.


```{r}

dm.noev4 <- dm[-unique(m1.diag$extremeValues[EffectType == "Multivariate Random Effect ID"]$Index)]

``` 

If you wanted even more customized options, you could achieve those by
making the subsetting of the extreme values dataset fancier (e.g.,
picking multiple but not all effect types, etc.)^[This is one of the
benefits of `R` and an approach where everything is an object. We can
use the output from one function, `modelDiagnostics()`, and because it
returns an object, a dataset of extreme values, we can operate on it
as we want and then use those results to subset our main dataset for
analysis.].

At this point, we can re-run our analysis using the revised dataset.
For the sake of example, I will just use our dataset that
excluded only one person, `dm.noev2`.

```{r, fig.width = 7, fig.height = 10}

m1noev <- lmer(dStress ~ dEnergy + (dEnergy | ID),
           data = dm.noev2)

summary(m1noev)

m1noev.diag <- modelDiagnostics(m1noev, ev.perc = .01)
plot(m1noev.diag, ncol = 2, nrow = 3, ask=FALSE)

```

The results look fairly similar, although the correlation between the
random intercept and slope has dropped from 
`r as.data.table(VarCorr(m1))[3, sdcor]`
to 
`r as.data.table(VarCorr(m1noev))[3, sdcor]`
and the average slope of energy on stress, the fixed effect part, is
somewhat weaker going from 
`r fixef(m1)[["dEnergy"]]` to 
`r fixef(m1noev)[["dEnergy"]]`.

The diagnostics also look improved.

A shorter trick to re-running the same model on a new dataset is to
use the `update()` function. The `update()` function takes an existing
model, and you can update that model in different ways. Here, we will
update the model by just using a new dataset, here the dataset where
we only excluded a different ID, `dm.noev2a`. This is just to
illustrate how easy it is even without knowing the formula for a model
to just update an old model on a different dataset. This is very
helpful for running models with and without extreme values or with and
without some participants who perhaps did not follow the study
protocol correctly, etc. The benefit of using `update()` is that even
if you have lots of predictors in your model, you don't have to type
them all again and you are guaranteed to have the same model as
before, just with a new dataset, no chance for typos and forgetting a
predictor or covariate.

```{r}

m1noev2 <- update(m1, data = dm.noev2a)
summary(m1noev2)
# Note - you could do this theoretically with any of the dm.noev models we created.
# It just so happens that when we remove too much in our dataset, the models don't coverge!
```

In the previous topic, we added both between and within person
variants of a predictor into a LMM as fixed effects. Now let's take a
look at doing that with both fixed and random slopes.
First, we need to create a between and within person version of our
predictor, `dEnergy`. Note that this only works because energy was measured
each day (which is why it has the d at the front of the variable name).
It would not work if `dEnergy` was measured only once, it would
already be a between person variable. Note also that mean deviations
only makes sense for continuous predictors. If `dEnergy` was binary, it
would not make sense to look at the average and the deviations from
the average, probably.


```{r}

dm[, c("Benergy", "Wenergy") := meanDeviations(dEnergy), by = ID]

```

Now we can fit a LMM with `Benergy` and `Wenergy` as fixed effects and a
random intercept and random slope for `Wenergy`. Note that the random
intercept is included automatically, without us having to ask for it
explicitly. We also look at the model diagnostics and confidence
intervals.

```{r, fig.width = 7, fig.height = 10, warning = FALSE}

iccMixed("dStress",id="ID",data=dm) ## Let's get the ICC of both
iccMixed("dEnergy", id ="ID", data=dm)

m2 <- lmer(dStress ~ Benergy + Wenergy +
             (Wenergy | ID),
           data = dm)

summary(m2)

m2.diag <- modelDiagnostics(m2)
plot(m2.diag, ncol = 2, nrow = 3, ask=FALSE)

m2.ci <- confint(m2, oldNames = FALSE)
m2.ci

```

In the output, we now see the random slope is for `Wenergy` and under
fixed effects we have both `Benergy` and `Wenergy`. In this instance both
the sign and the magnitude of the fixed effects for the slope of
`Benergy` and `Wenergy` on stress are about the same. However, this does
not have to be true. In fact, its even possible for the between person
effect to have one sign and the within person effect to have the
opposite sign for the slope^[A conceptual example of why the between
and within might differ. Imagine that you look at the relationship
between exercise and well being. You might find that people that
exercise more on average have better average wellbeing. However, if
people exercise more than usual for them, they might over exercise and
that be associated with worse wellbeing. Because overexercising is
relative to an individuals' own fitness level, it only shows up at the
within person only level. Although this is relatively rare, it often
happens that the magnitude of effects varies across the between person
and within person effects.].

Visually, all the diagnostics look relatively OK, although there are a
few extreme values (now using the default top and bottom 0.5% of the
theoretical distribution definition) on the residuals and on the
random intercept. In practice, it would be worthwhile to consider
whether these extreme values should be addressed somehow or if you are
comfortable proceeding as is.

The profile likelihood confidence intervals take a few seconds to
complete and generate a few warnings related to the fact that the
lower bound for the intercept-slope correlation is -1, the lowest
possible correlation (MLB: I have suppressed these in the markdown output).

### Sample Write Up

A linear mixed model was fit to 
`r nobs(m2)` stress scores from 
`r as.integer(ngrps(m2))` people. The intraclass correlation
coefficient of stress, the outcome, was 
`r iccMixed("dStress", id = "ID", data = dm)$ICC[1]` and of energy, the
predictor was 
`r iccMixed("dEnergy", id = "ID", data = dm)$ICC[1]`,
indicating that about 40% of the total variance in stress and about
28% of the total variance in energy exists between people with the
remaining due to daily fluctuations within people.
The fixed effect intercept revealed that the
estimated average [95% CI] stress was 
`r fixef(m2)[["(Intercept)"]]`
`r sprintf("[%0.2f, %0.2f]", m2.ci[5, 1], m2.ci[5, 2])`, when `Benergy`
and `Wenergy` are zero.
However, there were individual differences, with the standard
deviation for the random intercept being
`r as.data.table(VarCorr(m2))[1, sdcor]`
indicating that there are individual differences. Assuming the random
intercepts follow a normal distribution, 
we expect most people to fall within one standard deviation of the
mean, which in these data would be somewhere between:
`r fixef(m2)[["(Intercept)"]] + c(-1, 1) *  as.data.table(VarCorr(m2))[1, sdcor]`. 
There was also a significant fixed effect of average energy on stress,
such that each one unit higher average energy people had was associated
with `r fixef(m2)[["Benergy"]]` 
95% CI = `r sprintf("[%0.2f, %0.2f]", m2.ci[6, 1], m2.ci[6, 2])` lower stress.

There was no significant fixed effect of within person energy on
stress,`r fixef(m2)[["Wenergy"]]` 
95% CI = `r sprintf("[%0.2f, %0.2f]", m2.ci[7, 1], m2.ci[7, 2])`.

(If this would have been significant, you could report that: On 
days where people were one unit higher on energy than their 
own average, people were expected to have 
`r fixef(m2)[["Wenergy"]]` 
95% CI = `r sprintf("[%0.2f, %0.2f]", m2.ci[7, 1], m2.ci[7, 2])`
lower stress that same day, on average. However, there were
individual differences, with the standard deviation for the
random slope being
`r as.data.table(VarCorr(m2))[2, sdcor]`
indicating that there are individual differences in the slope of
within person energy on stress. 
Assuming the random slope follow a normal distribution,
we expect most people to fall within one standard deviation of the
mean, which in these data would indicate that most people are expected
to have a within person slope of stress on energy between:
`r fixef(m2)[["Wenergy"]] + c(-1, 1) *  as.data.table(VarCorr(m2))[2, sdcor]`.)

Finally, there was a negative correlation between the random intercept
and slope, 
`r as.data.table(VarCorr(m2))[3, sdcor]`
indicating that compared to the population averages, people who had a
relatively higher stress when energy was 0 also tended to have a
more negative within person association between energy and stress.

# Graphing LMMs

We can graph results from LMMs much as we do for linear regressions or
GLMs. There are a few nuances, however. In LMMs, people sometimes
refer to marginal or conditional effects. In some ways, these
correspond to the idea of fixed or random effects.

Marginal predictions are based only on the averages, basically the fixed
effects portion of the model. Conditional predictions incorporate
**both** the fixed and random effects portion of the model.

Because the conditional predictions incorporate random effects, they
can only be made for the existing data because it is impossible to
know what a new person would look like (e.g., if you recruited one
more participant, would their BLUP be like ID 1, ID 2, or ID XX?).
Therefore, conditional predictions are only made off the data used to
build the model. Marginal predictions could be made for new data,
though to answer a question like, what would the model predict, on
average, the stress score would be if mood was 2 points above average
on a day? That would not tell you how any single individual is
predicted to be, but does indicate what, on average, the prediction
would be from the model. Additionally, because conditional predictions
will vary by ID, we can get different conditional predictions for each
participant in the dataset.
We will look at examples of both of these briefly.

## Marginal

Graphs based off the marginal predictions from the model are probably
the default and most common approach. These graphs show the average
association or predictions from the model and basically graph the
fixed effects part of the model. They are interpreted nearly
identically to linear regression graphs.

To start, we will use our model with between person and within person
energy predicting stress with a random intercept and random slope for
within person energy, called `m2`. We use the `visreg()` function to
graph the results with the x axis being `Benergy` and the y axis the
predicted stress scores.

```{r, error = TRUE}

m2 <- lmer(dStress ~ Benergy + Wenergy +
             (Wenergy | ID),
           data = dm)

visreg(m2, xvar = "Benergy",
       partial = FALSE, rug = FALSE,
       gg = TRUE,
       xlab = "Average Energy",
       ylab = "Predicted Stress",
       line = list(col = "black", linewidth = 1)) + 
  theme_pubr()

``` 

The graph shows us how average energy (between person energy) and stress
are associated on average. We made a number of customizations to the
graph, including not showing partial residuals, no rug plot, defining
our own x and y axis labels, and changing the line colour to black
(the default is blue), and use the `pubr` theme.


We can make the same plot for within person energy. We just change the
`xvar` option to be `Wenergy`. We also make a bit fancier x axis
label. Note the use of "\n" that tells `R` that we want a line break
in the x axis label.

```{r}

visreg(m2, xvar = "Wenergy",
       partial = FALSE, rug = FALSE,
       gg = TRUE,
       xlab = "Within Energy\n(deviations from own average)",
       ylab = "Predicted Stress",
       line = list(col = "black", linewidth = 1)) + 
  theme_pubr()

``` 

This figure shows the negative fixed effect coefficient for `Wenergy` graphically.

## Conditional

We can graph conditional predictions as
well. We will make a few modifications compared to the marginal
graph.

We specify specific IDs should be plotted.
To see available IDs, we look at `unique(dm$ID)`. I've
chosen a few IDs to show for example. You could pick others.

For this we generate predicted values using the `predict()` function.
We specify the `re.form` argument. It defaults to including 
all random effects, but more specifics can be chosen. See:
`?predict.merMod` for more details.

To make predictions, we have to pick values for every predictor.
To start, we use a sequence of values for `Wenergy` from the minimum to the
maximum of `Wenergy` in the dataset. We also need to pick a value for
`Benergy`. We can use the mean of `Benergy` in the dataset.
If there were other predictors or covariates, we would need to 
pick values for those as well. Often, the mean/median or the 
most frequent value is used for categorical predictors.

We use the `expand.grid()` function to create a combinations
of our predictor values and IDs.
First we predict using all random effects, and save that in a variabled called
`YhatAllRE`. We then use the `re.form` argument to specify
we only want the random intercept, which we save in a variable 
called `YhatIntRE`. Finally, we set `re.form = NA` to
indicate we do not want any random effects, which we save in
a variable called `YhatNoRE`.

```{r}

unique(dm$ID)

newdata <- expand.grid(
  Wenergy = seq(from = min(dm$Wenergy), to = max(dm$Wenergy), length.out = 10),
  Benergy = mean(dm$Benergy),
  ID = c(10, 35, 54, 64, 99, 101))

## default is include all random effects
newdata$YhatAllRE <- predict(
  m2, newdata = newdata,
  type = "response")

## only random intercept
newdata$YhatIntRE <- predict(
  m2, newdata = newdata,
  re.form = ~ (1 | ID),
  type = "response")
## no random effects
newdata$YhatNoRE <- predict(
  m2, newdata = newdata,
  re.form = NA,
  type = "response")

```

Next we can plot the results using `ggplot2` predicting using all 
random effects.
Note that these are not raw regression lines per person, they are
conditional predictions from the LMM, so they will include shrinkage
of the random intercept and slopes to the overall sample average
intercept and slope. Nevertheless, we can see that there are quite
large differences in the within person association between energy and
stress for say ID 10 compared to ID 101. They are not the same, and
although in the marginal plots section earlier, we saw that on average
there is a negative association between within person energy and stress,
not everyone exactly has a negative predicted slope (10 looks a little positive).

```{r, fig.fullwidth = TRUE}

ggplot(newdata, aes(x = Wenergy, y = YhatAllRE, colour = factor(ID))) +
  geom_line() +
  labs(x = "Within Energy\n(deviations from own average)",
       y = "Predicted Stress") +
  ggtitle("Conditional Random Intercept and Slope") +
  theme_pubr() + facet_wrap(~ID)
```

Here we plot using the random intercept only.
What this figure showed is each individual IDs own predicted
intercept, but it used the average marginal slope for within person
mood and stress for all IDs.

```{r, fig.fullwidth = TRUE}
ggplot(newdata, aes(x = Wenergy, y = YhatIntRE, colour = factor(ID))) +
  geom_line() +
  labs(x = "Within Energy\n(deviations from own average)",
       y = "Predicted Stress") +
  ggtitle("Conditional Random Intercept Only") +
  theme_pubr() + facet_wrap(~ID)
```

Finally we plot using no random effects.
We get the marginal intercept and marginal slope, which is just the same as
our marginal graph but repeated for each ID, a fairly useless figure,
but it highlights what the marginal (or fixed effects only) predictions mean.

```{r, fig.fullwidth = TRUE}
ggplot(newdata, aes(x = Wenergy, y = YhatNoRE, colour = factor(ID))) +
  geom_line() +
  labs(x = "Within Energy\n(deviations from own average)",
       y = "Predicted Stress") +
  ggtitle("No Random Effects") +
  theme_pubr() + facet_wrap(~ID)
```


# Optional Only: Solving Convergence/Estimation Issues

We have talked before some about convergence, where algorithms iterate
through cycles and try to find the best parameter estimates. Sometimes
this process fails or runs into problems, which we may broadly call
convergence or estimation issues.

Sometimes these do not make too big of a difference, sort of false
positives, but sometimes these issues may severely impact a models'
results. If you see warnings about convergence or estimation, its best
to address and resolve them before being confident in the model
results.

Let's create a between and within version of the daily variable,
`dMood` and look at a fixed and random slope LMM predicting stress
from mood

```{r}

dm[, c("Bmood", "Wmood") := meanDeviations(dMood), by = ID]

# Please note, for demonstration purposes, I'm using a different outcome variable, "stress"
# rather than dStress (as that model converged fine) - stress is just the average across days
mest <- lmer(stress ~ Bmood + Wmood +
             (Wmood | ID),
             data = dm)

``` 

When we estimate the model, we get a message about a boundary
(singular) fit with the suggestion to see `?isSingular`.
This is a helpful suggestion and provides more details on potential
approaches to resolving it. If we run 
`?isSingular` in the `R` console will get a help page with some more
information, which you **should** do and go read now.

Armed with that knowledge, we can run a `summary()` on our model to
learn a bit more. Its not always clear, but in this case we can see a
correlation of 1 between the random intercept and slope that is
causing the singularity.

```{r}

summary(mest)

```

Our options for resolving this are listed in `?isSingular`. In this
case, a simple path forward is to remove some random effects. We
basically always keep a random intercept in our LMMs, so the only
candidate to remove is the random slope for `Wmood`. 
Sometimes when doing that we no longer get the singularity warning, but
in this case, the model still failed to converge, so we should not trust 
the results.

```{r}

mest2 <- lmer(stress ~ Bmood + Wmood +
             (1 | ID),
             data = dm)

summary(mest2)

``` 

If you have multiple random slopes, it can be tricky to decide which
one to drop.

Now let's look at a convergence issue. For this, we'll use another
dataset, the `aces_daily` data. First we setup the data (from the
`JWileymisc` package) and then fit a model with a random slope of
stress predicting negative affect.

```{r}

data(aces_daily)

mconv <- lmer(NegAff ~ STRESS + (1 + STRESS | UserID),
          data = aces_daily)

``` 

In this case, we get a warning message that the model did not
converge. This is not a singularity issue. This means that the
algorithm that tries to find the best parameter estimates tried and
failed to find estimates it's confident are 'best'. For now, we won't
worry too much about what the 'best' actually means. In this case,
sometimes using a different optimizer, basically the algorithm that
goes about trying to find the best estimates, or changing its options
can help. We can control the underlying algorithms using the
`lmerControl()` function. In this unit, we are not going to talk much
about the different options but just show one alternate you could try
if you run into convergence problems.


```{r}

## use the control function to pick a  different algorithm
## and adjust the options to be a bit stricter
## may take longer but also may have a better chance of converging
strictControl <- lmerControl(optCtrl = list(
   algorithm = "NLOPT_LN_NELDERMEAD",
   xtol_abs = 1e-12,
   ftol_abs = 1e-12))

## now fit our model
mconv2 <- lmer(NegAff ~ STRESS + (1 + STRESS | UserID),
               data = aces_daily,
               control = strictControl)

```

In this case, the change in algorithm and options did it and we now
get our model converging without warnings. Sometimes it may still have
warnings and none of what you know how to try will solve those. In
those cases, you'd either want to consult an expert or simplify the
model. Very often, much as with singularity issues, convergence issues
surround the random effects in a model. If you remove some of the
random effects (at the extreme a random intercept only) you are likely
to resolve the convergence issues, although you may give up aspects of
the model you wanted to keep.

Sometimes as well, the convergence issues may not make much
difference. Here we can look at the model results that did not
converge quite and those that did.

```{r}

summary(mconv)

summary(mconv2)

``` 

In this instance, the only apparent difference is in the 4th decimal
point of the standard deviation of the random intercept.
While the first algorithm did not quite achieve estimates it was
confident are the best, it was really very close. Sometimes however,
the estimates could be far off. If you cannot get a model that
converges, though, you have no comparison so its hard to know if the
non convergence was an issue or not. That is why its generally the
best idea to always resolve a non convergence issue if possible and be
very cautious about interpreting a model that did not converge.