Week10_LMM_Moderation.qmd

---
title: "PSY4210: Linear Mixed Models (LMMs) - Moderations & Comparisons"
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)

## emmeans is a new package

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

## load data collection exercise data
## merged is a a merged long dataset of baseline and daily
dm <- as.data.table(read_sav("[2021] PSY4210 merged.sav"))

## Remind R which of our variables are factors

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

dm[, relsta := factor(
  relsta, levels = c(1,2,3), 
  labels = c("single", "in a committed exclusive relationship", "in a committed nonexclusive relationship"))]

```

# LMM Notation

Let's consider the formula for a relatively simple LMM:

$$
y_{ij} = b_{0j} + b_1 * x_{1j} + b_2 * x_{2ij} + \varepsilon_{ij}
$$

Here as before, the *i* indicates the *i*th observation for a specific
unit (e.g., person but the unit could also be classrooms, doctors,
etc.) and the *j* indicates the *j*th unit (in psychology usually
person).

Regression coefficients, the $b$s with a *j* subscript indicate a
fixed and random effect. That is, the coefficient is allowed to vary
across the units, *j*. As before, these coefficients in practice are
decomposed into a fixed and random part:

$$
b_{0j} = \gamma_{00} + u_{0j}
$$

and we estimate in our LMM the fixed effect part, $\gamma_{00}$, and
the variance / standard deviation of the random effect or the
covariance matrix if there are multiple random effects, $\mathbf{G}$: 

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

Regression coefficients without any *j*
subscript indicate fixed only effects, effects that do not vary across
units, *j*. These are fixed effects and get estimated directly.

Predictors / explanatory variables, the $x$s with an *i* subscript
indicate that the variable varies within a unit. Note that the
outcome, $y$ **must** vary within units to be used in a LMM.

In this case, the notation tells us the following:

- $y_{ij}$ the outcome variable, which varies both within and between
  people
- $b_{0j}$ the intercept regression coefficient, which is both a fixed
  and random effect
- $b_1$ the regression coefficient, slope, for the first predictor,
  which is a fixed effect only
- $x_{1j}$ the first predictor/explanatory variable, this is a between
  unit variable only, as the lack of an *i* subscript indicates it
  does not vary within units. It could not have a random slope.
- $b_2$ the regression coefficient, slope, for the second predictor,
  which is a fixed effect only
- $x_{2ij}$ the second predictor/explanatory variable, this variable
  varies within individuals as shown by its *i* subscript. It could
  have a random slope, although in this model, it only has a fixed
  effect slope.
- $\varepsilon_{ij}$ the residuals, these vary within and between
  units.

The following decision tree provides some guide to when a predictor /
explanatory variable can be a fixed and random effect.

```{r, echo = FALSE, fig.cap = "Type of effect decision tree"}

DiagrammeR::grViz('
digraph "Type of effect decision tree" {
  graph [overlap = true, fontsize = 12]
  node [fontname = Helvetica, shape = rectangle]

  variable [ label = "What level does your variable vary at?" ];
  between [ label = "Between Variable" ];
  within [ label = "Within Variable" ];

  fixed [ label = "Fixed Effect Only" ];
  random [ label = "Fixed & Random Effect" ];
  type [ label = "Do you want a fixed effect only?" ];

  variable -> between [ label = "only between units" ];
  variable -> within [ label = "varies within (+/- between) units" ];
  between -> fixed ;
  within -> type ;
  type -> fixed [ label = "yes" ];
  type -> random [ label = "no" ];
}
')

```

Let's see two examples of putting this basic model into practice.

$$
energy_{ij} = b_{0j} + b_1 * loneliness_{j} + b_2 * stress_{ij} + \varepsilon_{ij}
$$

The corresponding `R` code is:

```{r}

summary(lmer(dEnergy ~ loneliness + dStress + (1 | ID), data = dm))

``` 

Here is another example decomposing stress into a between and within component.

$$
energy_{ij} = b_{0j} + b_1 * Bstress_{j} + b_2 * Wstress_{ij} + \varepsilon_{ij}
$$

```{r}

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

summary(lmer(dEnergy ~ Bstress + Wstress + (1 | ID), data = dm))

``` 

We can make more effects random effects. For example, taking our
earlier example and just changing $b_2$ into $b_{2j}$:

$$
energy_{ij} = b_{0j} + b_1 * loneliness_{j} + b_{2j} * stress_{ij} + \varepsilon_{ij}
$$

The corresponding `R` code is:

```{r}

summary(lmer(dEnergy ~ loneliness + dStress + (dStress | ID), data = dm))

``` 

Now with two random effects, we assume that the random effects,
$u_{0j}$ and $u_{2j}$, which we collectively denote
$\mathbf{u}_{j}$ follow a multivariate normal distribution with
covariance matrix $\mathbf{G}$.

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

Based on the little decision chart, between unit only variables, like
$loneliness_j$ and $Bstress_j$ *cannot* be random effects. In the
data collection exercise, we measured loneliness at baseline and also
in the daily diary questionnaires. In this example we are using the 
baseline (trait) loneliness and not the daily one. Also, while it is
technically possible for something to only be a random effect without
a corresponding fixed effect, its not common and not recommended as it
would be equivalent to assuming that the fixed effect, the mean of the
distribution, is 0, which is rarely appropriate.

# Interactions in LMMs

Interactions in LMMs work effectively the same way that interactions
in GLMs do, although there are a few nuances in options and possible
interpretations.
Using the notation from above, let's consider a few different possible
interactions. 

## Cross Level (Between and Within Unit) Interactions

First, let's take our model with loneliness and stress and
include an interaction. Here is the model without an interaction.

$$
energy_{ij} = b_{0j} + b_1 * loneliness_{j} + b_2 * stress_{ij} + \varepsilon_{ij}
$$

The corresponding `R` code is:

```{r}

summary(lmer(dEnergy ~ loneliness + dStress + (1 | ID), data = dm))

``` 

Now let's add the interaction, as a fixed effect.

$$
energy_{ij} = b_{0j} + b_1 * loneliness_{j} + b_2 * stress_{ij} + 
  b_3 * (loneliness_{j} * stress_{ij}) + 
  \varepsilon_{ij}
$$

The corresponding `R` code is:

```{r}

## long way
summary(lmer(dEnergy ~ loneliness + dStress + loneliness:dStress + (1 | ID), data = dm))

## short hand in R for simple main effect + interaction
## identical, but shorter to the above
summary(lmer(dEnergy ~ loneliness * dStress + (1 | ID), data = dm))

``` 

The relevant, new, part is the interaction term, $b_3$, a fixed effect
in this case. If we focus just on that one term, we see that the
coefficient, $b_3$ is applied to the arithmetic product of two
variables, here loneliness and stress. As it happens, one of them, loneliness,
varies between units whereas the other, stress, varies within
units. You will sometimes see this termed as "cross level" interaction
between it involves a between and within varying variable.

$$
b_3 * (loneliness_{j} * stress_{ij})
$$

As with interactions for regular GLMs, interactions in LMMs can be
interpretted in different ways. The two common interpretations are
easiest to see by factoring the regression equation.
Here are three equal equations that highlight different ways of
viewing the interaction.

In the latter two formats, it highlights how the simple effect of
stress varies by loneliness and how the simple effect of loneliness varies by
stress. 

$$
\begin{align}
energy_{ij} &= b_{0j} + b_1 * loneliness_{j} + b_2 * stress_{ij} + b_3 * (loneliness_{j} * stress_{ij}) + \varepsilon_{ij} \\
          &= b_{0j} + b_1 * loneliness_{j} + (b_2 + b_3 * loneliness_j) * stress_{ij} + \varepsilon_{ij} \\
          &= b_{0j} + (b_1 + b_3 * stress_{ij}) * loneliness_{j} + b_2 * stress_{ij} + \varepsilon_{ij} \\
\end{align}
$$

The nuance in LMMs comes in because some variables vary only between
units and others within units. For example, when interpretting the
interaction with respect to the simple effect of stress, we could say
that the association between daily stress and energy on the same day
depends on the loneliness of a participant. Conversely, when interpretting
with respect to the simple effect of loneliness, we could say that the
association of participant loneliness and average energy depends on how
stressed someone is feeling on a given day. Loneliness varies between
people, stress varies within people, so that must be taken into
account in the interpretation.

## Between Unit Interactions

The same approach would work with other type of variables in LMMs. For
example, here we have a model with loneliness and sex as predictors. Both
vary only between units.

$$
\begin{align}
energy_{ij} &= b_{0j} + b_1 * loneliness_{j} + b_2 * sex_{j} + b_3 * (loneliness_{j} * sex_{j}) + \varepsilon_{ij} \\
          &= b_{0j} + b_1 * loneliness_{j} + (b_2 + b_3 * loneliness_j) * sex_{j} + \varepsilon_{ij} \\
          &= b_{0j} + (b_1 + b_3 * sex_{j}) * loneliness_{j} + b_2 * sex_{j} + \varepsilon_{ij} \\
\end{align}
$$

When interpretting the interaction with respect to the simple effect of
sex, we could say that the association between participant sex and 
average energy depends on the loneliness of a participant. Conversely, 
when interpretting with respect to the simple effect of loneliness, we
could say that the association of participant loneliness and average 
energy depends on participant's sex.

## Within Unit Interactions

Finally, both variables could vary within units.

$$
\begin{align}
energy_{ij} &= b_{0j} + b_1 * selfesteem_{ij} + b_2 * stress_{ij} + b_3 * (selfesteem_{ij} * stress_{ij}) + \varepsilon_{ij} \\
          &= b_{0j} + b_1 * selfesteem_{ij} + (b_2 + b_3 * selfesteem_{ij}) * stress_{ij} + \varepsilon_{ij} \\
          &= b_{0j} + (b_1 + b_3 * stress_{ij}) * selfesteem_{ij} + b_2 * stress_{ij} + \varepsilon_{ij} \\
\end{align}
$$

When interpretting the interaction with respect to the simple effect
of stress, we could say that the association between daily stress and
energy on the same day depends on same day self-esteem level. 
Conversely, when interpretting with respect to the simple effect of
self-esteem, we could say that the association of daily self-esteem 
and same day energy depends on how stressed someone is feeling on a 
given day. 


# Continuous Interactions in `R`

Aside from the notes about some minor interpretation differences, in
general interactions in LMMs are analysed, graphed, and interpreted
the same way as for GLMs.

First to avoid any issues around diagnostics etc. from haven labeled
type data, we will convert the variable we are going to work with to
numeric. Then we fit a LMM with an interaction between stress and
neuroticism, energy as the outcome and a random intercept as the only
random effect.

```{r}

dm[, dSE := as.numeric(dSE)]
dm[, dMood := as.numeric(dMood)]
dm[, dEnergy := as.numeric(dEnergy)]
dm[, dStress := as.numeric(dStress)]
dm[, neuroticism := as.numeric(neuroticism)]

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

```

A quick check of the model diagnostics suggests that the
data look fairly good. The intercepts do not appear to 
follow a normal distribution that closely, partly due to 
the long left tail, but for now we will leave it. If you're interested
in how to fix left skews, please see the bonus material at the end of 
this markdown (optional).

```{r}

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

```

No extreme values are present. If there were any, we can remove that, 
as we have discussed in depth in previous lectures, update the model,
and re-run diagnostics. In practice it could take a few rounds of 
extreme value removal or you may decide to stop at one round.

Let's look at the summary of our model, *m*. Although we have used 
`summary()` a lot in the past, we'll introduce another function to 
help look at `lmer()` model results, `modelTest()`. In this lecture,
we will only learn and interpret part of its output, with the rest 
of the output from `modelTest()` covered later. In addition to get
nicely formatted results rather than a set of datasets containing 
the results, we use the `APAStyler()` function.

```{r}

APAStyler(modelTest(m))

``` 

The results show the regression coefficients, asterisks for p-values,
and 95% confidence intervals in brackets for the fixed effects, the
standard deviations of the random effects, the model degrees of
freedom, which is how many parameters were estimated in the model
total, and the number of people and observations. For now, we will
ignore all the output under row 9, N (Observations). In the case of
this model we can see the following.

A LMM was fit with 283 observations from 88 people. There was no
significant neuroticism x stress interaction, b [95% CI] = 0.07
[-0.03, 0.16], p = .160.

We can also pass multiple model results in a list together, which puts the
results side by side. This is particularly helpful for comparing
models with and without covariates, to evaluate whether removing
extreme values changed the results substantially, or to compare models
with different outcomes.

```{r}
# Remember: m <- lmer(dEnergy ~ neuroticism * dStress + (1 | ID), data = dm)

# Here's a new model with mood as the outcome to compare:

mtest1 <- lmer(dMood ~ neuroticism * dStress + (1 | ID), data = dm)

APAStyler(list(
  Energy = modelTest(m),
  Mood = modelTest(mtest1) ))

``` 

These results show us that we have similar results when predicting
daily energy and daily mood from stress and neuroticism. The 
relationship between daily stress and daily mood, and daily stress
and daily energy did not vary by individual differences in 
neuroticism. In this case, it would make sense to re-run these
models without the interaction term to test the main effects
of daily stress and neuroticism.

```{r}
mmain <- lmer(dEnergy ~ neuroticism + dStress + (1 | ID), data = dm)
mtest1main <- lmer(dMood ~ neuroticism + dStress + (1 | ID), data = dm)

APAStyler(list(
  Energy = modelTest(mmain),
  Mood = modelTest(mtest1main) ))

```

Results are indeed similar for daily mood and energy as outcomes. 
There are significant negative associations between daily stress
and both outcome variables, and also neuroticism and both outcome
variables.


## Plotting

Typically, to plot our significant interaction, a few exemplar 
lines are graphed showing the slope of one variable with the
outcome at different values of the moderator. 
As with GLMs, we can use the `visreg()` function.
Here, we'll use neuroticism as the moderator. A common 
approach to picking level of the moderator is to use the 
Mean - 1 SD and Mean + 1 SD. To do that, we first need the mean
and standard deviation of neuroticism, which we can get using 
`egltable()` after excluding duplicates by ID, since 
neuroticism only varies between units. Note that we are doing
this for the sake of example only since our interaction was
not signficant.

```{r}

egltable(c("neuroticism"), data = dm[!duplicated(ID)])

visreg(m, xvar = "dStress",
       by = "neuroticism", overlay=TRUE,
       breaks = c(3.46-1.11, 3.46+1.11),
       partial = FALSE, rug = FALSE)

```

The results show, in an easier to interpret way, what the positive
interaction coefficient of $b = 0.07$ means, people with higher levels
of neuroticism are less sensitive to the (negative) effects of stress. 
People higher in neuroticism are relatively less sensitive to the 
effects of stress, although in both cases, higher stress is associated 
with lower energy levels. Keep in mind again that we are interpreting
this as though the interaction was signficant for the sake of example
only.

Another common way of picking some exemplar values is to use the 25th
and 75th percentiles. These work particularly well for very skewed
distributions where the mean +/- SD could be outside the observed
range of the data. Again we exclude duplicates by ID and then use the
`quantile()` function to get the values, 3, and 4.5 for the 25th and
75th percentiles.

```{r}

quantile(dm[!duplicated(ID), neuroticism], na.rm = TRUE)

visreg(m, xvar = "dStress",
       by = "neuroticism", overlay=TRUE,
       breaks = c(3, 4.5),
       partial = FALSE, rug = FALSE)

``` 

## Simple Effects

When working with models that have interactions, a common aid to
interpretation is to test the simple effects / slopes from the
model. For example, previously we graphed the association between
stress and energy at M - 1 SD and M + 1 SD on neuroticism. 
However, although visually both lines appeared to have a
negative slope, we do not know from the graph alone whether there is a
significant association between stress and energy at both the
low (M - 1 SD) and high (M + 1 SD) levels of neuroticism. To answer
that, we need to test the simple slope of stress at specific values of
neuroticism. Again this is for demonstration only given our non-
significant moderation. We would not need to compute simple slopes
or effects for non-significant moderations in reality.

Our default model does actually give us one simple slope:
it is the simple slope of stress when $neuroticism = 0$. However, as
we can tell from the mean and standard deviation of neuroticism, 0 is
very far outside the plausible range of values so that simple slope
given to us by default from the model is not too useful. We could
either center neuroticism and re-run the model, which would get us a
different simple slope, or use post hoc functions to calculate simple
slopes.

We will use the `emtrends()` function from the `emmeans` package to
test the simple slopes. This function also works with GLMs, for your
reference.

The `emtrends()` function take a model as its first argument, then the
variable that you want to calculate a simple slope for, here `stress`,
the argument `at` requires a list of specific values of the moderator,
and then we tell it how we want degrees of freedom calculated (note
this only applies to `lmer` models). We store the results in an `R`
object, `mem` and then call `summary()` to get a summary table. The
`infer = TRUE` argument is needed in `summary()` if you want
p-values. 

```{r}

mem <- emtrends(m, var = "dStress",
                at = list(neuroticism = c(3.46-1.11, 3.46+1.11)),
                lmer.df = "satterthwaite")

summary(mem, infer=TRUE)

```

The relevant parts of the output, for us, are the columns for
`stress.trend` which are the simple slopes, the values of
`neuroticism` which tell us at what values of neuroticism we have
calculated simple slopes, the confidence intervals, `lower.CL` and
`upper.CL`, 95% by default, and the p-value. From these results, we
can see that when $neuroticism = 2.35$ there is a significant
negative association between stress and energy, but not when 
$neuroticism = 4.57$.

## Sample Write Up

With all of this information, we can plan out some final steps for a
polished write up of the results. First, let's get exact p-values for
all our results. We can do this through options to `pcontrol` in
`APAStyler()`. We also re-print the simple slopes here.

```{r}

APAStyler(modelTest(m),
  pcontrol = list(digits = 3, stars = FALSE, includeP = TRUE,
                  includeSign = TRUE, dropLeadingZero = TRUE))

summary(mem, infer=TRUE)

``` 

Now we will make a polished, finalized figure. I have customized the
colours, and turned off the legends. In place of legends, I have
manually added text annotations including the simple slopes and
confidence intervals and p-values for the simple slopes^[For your
reference, it took about 8 trial and errors of different x and y
values and angles to get the text to line up about right. I did not
magically get the values to use to get a graph that I thought looked
nice. That is why I think sometimes it is easier to add this sort of
text after the fact in your slides or papers rather than building it
into the code.].

```{r}

visreg(m, xvar = "dStress",
       by = "neuroticism", overlay = TRUE,
       breaks = c(3.46-1.11, 3.46+1.11),
       partial = FALSE, rug = FALSE, gg = TRUE,
       xlab = "Daily Stress",
       ylab = "Predicted Daily Energy") +
  scale_color_manual(values = c("2.35" = "black", "4.57" = "grey70")) +
  theme_pubr() +
  guides(colour = "none", fill = "none") +
  annotate(geom = "text", x = 3.2, y = 3.9, label = "High Neuroticism: b = -0.07 [-0.23, 0.10], p = .447",
           angle = -9) + 
  annotate(geom = "text", x = 4, y = 4.4, label = "Low Neuroticism: b = -0.21 [-0.36, -0.06], p = .005",
           angle = -22)

```

A linear mixed model using restricted maximum likelihood was used to
test whether the association of daily stress on daily energy is
moderated by baseline neuroticism scores. All predictors were included
as fixed effects and a random intercept by participant was included.
Visual diagnostics showed that energy was normally distributed, and no
outliers were present.

The daily stress x neuroticism interaction was not statistically
significant, which indicated that the relationship between stress 
and energy did not vary by neuroticism. Results from the analysis with
the interaction dropped revealed that both neuroticism and daily stress
were negatively associated with daily energy. 


# Continuous x Categorical Interactions in `R`

Continuous x Categorical interactions are conducted much as continuous
x continuous interactions. Typically with continuous x categorical
interactions, simple slopes for the continuous variable are calculated
at all levels of the categorical variable.

Let's illustrate this with a model examining the relationship between
daily energy levels and self-esteem with sex as a moderator.


```{r}

mconcat <- lmer(dSE ~ dEnergy*sex + (1| ID), data = dm) 

```

The model diagnostics look relatively good, albeit not perfect.

```{r}

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

```

With reasonable diagnostics, we can look at a summary.
There is one extreme residual but I'm choosing to 
leave it in the dataset.


```{r}

summary(mconcat)

```

Factor variables in interactions do not work currently with 
`modelTest()`, so if we wanted to use it, we'd need to manually dummy
code the categorical variable. The results are identical.

```{r}

dm[, female := as.integer(sex == "female")]

malt <- lmer(dSE ~ dEnergy * female + (1 | ID), data = dm)

APAStyler(modelTest(malt))

``` 

## Plotting

With continuous x categorical interactions, the easiest approach is to
plot the simple slope of the continuous variable by the categorical
one as shown in the following.

```{r}

visreg(mconcat, xvar = "dEnergy",
       by = "sex", overlay=TRUE,
       partial = FALSE, rug = FALSE)

```

## Simple Effects

When working with models that have interactions, a common aid to
interpretation is to test the simple effects / slopes from the
model. For example, previously we graphed the association between
daily energy and self-esteem at each level of the categorical
`sex` variable, i.e. for men and women.
However, we cannot tell from the graph whether daily energy is
significantly associated with self-esteem for men or women.

To answer that, we need to test the simple slope of energy at 
the two sex levels. Our default model does actually give us one simple slope:
it is the simple slope of energy when for men (i.e when female = 0), but we
might want more.

We will use the `emtrends()` function from the `emmeans` package to
test the simple slopes. 

The `emtrends()` function take a model as its first argument, then the
variable that you want to calculate a simple slope for, here `energy`,
the argument `at` requires a list of specific values of the moderator,
and then we tell it how we want degrees of freedom calculated (note
this only applies to `lmer` models). We store the results in an `R`
object, `mem` and then call `summary()` to get a summary table. The
`infer = TRUE` argument is needed in `summary()` if you want
p-values. 

```{r}

mem <- emtrends(mconcat, var = "dEnergy",
                at = list(sex = c("male", "female")),
                lmer.df = "satterthwaite")

summary(mem, infer=TRUE)

```

The relevant parts of the output, for us, are the columns for
`dEnergy.trend` which are the simple slopes, the values of
`sex` which tell us at what values of sex we have calculated 
simple slopes, the confidence intervals, `lower.CL` and 
`upper.CL`, 95% by default, and the p-value. From these results,
we can see that daily energy is significantly associated with 
self-esteem for any sex, although it is stronger for male participants than 
female participants.

# Categorical x Categorical Interactions in `R`


Categorical x Categorical interactions are conducted comparably,
although more contrasts / simple effect follow-ups are possible.

Here we will work with Int_Str again, which is a two-level 
categorical predictor (0 = no interaction with a stranger, 1 = 
interacted with a stranger that day) and energy as the outcome.
We also work with a three-level conscientiousness variable.


```{r}

## create a categorical conscientiousness variable
dm[, cons3 := cut(conscientiousness, breaks = quantile(conscientiousness, probs = c(0, 1/3, 2/3, 1), na.rm=TRUE),
                    labels = c("Low", "Mid", "High"),
                    include.lowest = TRUE)]

mcat2 <- lmer(dEnergy ~ cons3 * Int_Str + (1 | ID), data = dm)

```

The model diagnostics look relatively good.

```{r}

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

```

With reasonable diagnostics, we can look at a summary.

```{r}

summary(mcat2)

```

## Plotting

With categorical x categorical interactions, `visreg()` produces OK
but not great figures as shown in the following. We can see the means
of energy for all 6 cells (the cross of 3 level of
conscientiousness x 2 levels of Int_Str).

```{r}

dm[, Int_Str := factor(
  Int_Str,
  levels = c(0,1),
  labels = c("No interaction with stranger", "Interacted with stranger"))]

visreg(mcat2, xvar = "Int_Str",
       by = "cons3", overlay=TRUE,
       partial = FALSE, rug = FALSE)

```

## Simple Effects

When working with two categorical interactions (or with
a categorical predictor with >2 levels where you want to test various
group differences), the `emmeans()` function from the `emmeans`
package is helpful. We can get the means of interaction with stranger by
conscientiousness group and get confidence intervals and p-values. These p-values
test whether each mean is different from zero, by default.

```{r}

## re-run mcat2 now with int_str as factor
mcat2 <- lmer(dEnergy ~ cons3 * Int_Str + (1 | ID), data = dm)

em <- emmeans(mcat2, "Int_Str", by = "cons3",
              lmer.df = "satterthwaite")
summary(em, infer = TRUE)

```

A nice plot, with confidence intervals for the fixed effects, can be
obtained by using the `emmip()` function from the `emmeans`
package. It takes as input the results from `emmeans()`, not the
`lmer()` model results directly. Here is a simple plot showing the
categorical interactions. Note that with this approach, you could
basically fit the same model(s) that you would with a repeated
measures or mixed effects ANOVA model, with the advantage that LMMs do
not require balanced designs and allow both categorical and continuous
predictors (e.g., you could include continuous covariates
easily). GLMs and (G)LMMs can do everything that t-tests and various
ANOVAs can, but with greater flexibility.

```{r}

emmip(em, cons3~Int_Str, CIs = TRUE) +
  theme_pubr() +
  ylab("Predicted Energy")

```

# Comparing Models

For many statistical models, including LMMs, it can be informative to
compare different models. Comparing different models can be used in
lots of different ways. Here are some examples, although they are not
meant to be exhaustive.

* Evaluate which of two (or more) models provides the best fit to the
  data
* Evaluate / compare how the results for a particular predictor(s) of
  interest change across two (or more) models
* Calculate the significance of multiple predictors
* Calculate effect sizes for one or multiple predictors

We will look at examples of the different uses of model comparisons in
this topic. 

Just as there are many different kinds of models we can fit, even with
LMMs (e.g., with or without random slopes, etc.), so to there are many
different kinds and purposes for different model comparisons. 

To begin with, let's imagine we are just trying to compare two models:
Model A and Model B. We can broadly classify the type of comparison
based on whether Model A and B are nested or non-nested models. We
will talk about what these mean next.

## Nested Models

Models are considered nested when one model is a restricted or
constrained version of another model. For example, suppose that 
**Model A** predicts mood from stress and loneliness whereas 
**Model B** predicts mood from loneliness only. Written as a formula,
these could be:

$$
ModelA: mood_{ij} = b_{0j} + b1 * loneliness_j + b2 * stress_{ij} + \varepsilon_{ij}
$$

and

$$
ModelB: mood_{ij} = b_{0j} + b1 * loneliness_j + 0 * stress_{ij} + \varepsilon_{ij}
$$

In **Model B**, I purposely used $0 * stress_{ij}$ to highlight that
when a predictor is left out of a model, it is the same as fixing
(sometimes also called constraining) the coefficient ($b2$ in this
case) to be exactly 0. In this case, we would say that 
**Model B** is nested within **Model A**. In other words, 
**Model A** contains every predictor and parameter in **Model B** plus
more.

Restricted 

This idea is similar to the idea of nested data used in LMMs, but the
difference is that we are not talking about observations or
data, rather we are talking about the parameters of a model.

To summarize, briefly, we say that **Model B** is *nested* within
**Model A** if:

* setting one or more parameters in **Model A** to 0 yields the same
  model as **Model B**.
* both models use the same data, have the same number of observations,
  etc. Even if the same dataset is used in `R`, this does not
  guarantee the same data are used because if additional predictors
  are included in **Model A** and these extra predictors have some
  missing data, by default `R` would drop the missing data and result
  in a subset of cases being used for **Model A** than for **Model B**
  
If two models are nested, then we have the most options in terms of
comparing the two models. For example, we can evaluate whether 
**Model A** is a statistically significantly better fit than is
**Model B** using a Likelihood Ratio Test (LRT)^[https://en.wikipedia.org/wiki/Likelihood-ratio_test]. 

We can compare the fit of each model and use the difference in fit to
derive effect sizes. We also can attribute any differences in the two
models to the parameter(s) that have been constrainted to 0 in **Model
B** from **Model A**.

A simple definition of the LRT test statistic, $\lambda$ is based on
two times the difference in the log likelihoods.

$$
\lambda = -2 * (LL_B - LL_A)
$$

You may wonder why the likelihood *ratio* test can be conducted by
taking the difference in the log likelihoods. It is because the log of
a ratio is the same as the difference in the logs of the numerator and
denominator.

$$
log_{e}\left(\frac{6}{2}\right) = log_{e}(6) - log_{e}(2)
$$

which we can confirm is true in `R`:

```{r}

## log of ratio
log(6/2)

## difference in logs
log(6) - log(2)

```

If the null hypothesis of no difference is true in the population,
then $\lambda$ will follow a chi-squared distribution with degrees of
freedom equal to the number of parameters constrained to 0 in **Model
B** from **Model A**, the difference in degrees of freedom used for
each model, that is:

$$
\lambda \sim \chi^2(DF_A - DF_B)
$$

Thus we often use a chi-square distribution in the LRT to look up the
p-value. Finally, note that because LRTs are based on the log
likelihoods (LL) from a model, we need true log likelihoods for the
LRT to be valid. Therefore, we cannot use restricted maximum
likelihood, we need to use maximum likelihood estimation.

### Nested Models in `R`

To see the idea of nested models and the LRT in action, let's examine
a concrete example in `R`. Here are two LMMs corresponding to 
**Model A** and **Model B** formula we wrote previously. We can see in
the `R` code that the models are nested, the only difference is
`stress`. We set `REML = FALSE` to get maximum likelihood estimates so
that we get true log likelihoods. We also need to confirm that the two
models are based on the same number of observations. We can extract
just this information in `R` using the `nobs()` function. This lets
us confirm that the two models are fitted to the same data. For
example, if stress had some missing data that was not missing loneliness or
mood, it could be that **Model A** is based on fewer observations than 
**Model B**.

```{r}

modela <- lmer(dMood ~ loneliness + dStress + (1 | ID), data = dm, REML = FALSE)
modelb <- lmer(dMood ~ loneliness + (1 | ID), data = dm, REML = FALSE)

nobs(modela)
nobs(modelb)

``` 

In this case, we can see that the number of observations are
identical. Now we can find the log likelihoods of both models by using
the `logLik()` function.

```{r}

logLik(modela)
logLik(modelb)

``` 

Now we have the log likelihoods (LL) of each model and the degrees of
freedom from both models. From this, we can calculate $\lambda$ and
then look up the p-value for a chi-square distribution with $\lambda$
and 1 degree of freedom (1 is the difference in degrees of freedom
between Model A and B). To get the p-value from a chi-square, we use
the `pchisq()` function.

```{r}

## lambda
-2 * (-432.79 - -396.81)

## p-value from a chi-square
pchisq(71.96, df = 1, lower.tail = FALSE)

``` 

In practice, we do not need to do these steps by hand, we can get a
test of two nested models in `R` using the `anova()` function (which
in this case is not actually analyzing the variance really). We use
the `anova()` function with `test = "LRT"` to get a likelihood ratio
test.

```{r}

anova(modela, modelb, test = "LRT")

```

If you compare the chi-square value, degrees of freedom and p-value,
you'll see that they basically match what we calculated by hand. The
small differences are due to rounding (`R` will use more decimals of
precision whereas we only used two decimals).

In this simple case, we are only testing a single parameter, the fixed
regression coefficient for stress, because that is the only parameter
that differs between **Model A** and **Model B**. Thus in this case,
it would be easier to rely on the t-test we get from a summary of the
model. 

```{r}

summary(modela)

``` 

In this instance, the LRT and the t-test actually yield equivalent
p-values, 2e-16, however, this does not have to be. The LRT is based
on slightly different theory than the t-test and the t-test uses
approximate degrees of freedom for the t-distribution based on the
Satterthwaite method, whereas the LRT does not directly incorporate
the sample size in the same way. The two methods will be
assymptotically equivalent (at very large sample sizes) but can
differ, particularly for smaller samples.

In practice, we wouldn't usually use a LRT to evaluate whether a
single model parameter is statistically significant. The benefit of
(nested) model comparisons is that it allows us to compare two
models. Those models can be quite different.

Here are another two models, this time they differ by two predictors,
energy and sex. The LRT now tests whether the addition of energy
and sex results in significantly better fit for **Model A** than 
**Model B**.

```{r}

modela <- lmer(dMood ~ loneliness + stress + dEnergy + sex + (1 | ID),
               data = dm, REML = FALSE)
modelb <- lmer(dMood ~ loneliness + stress + (1 | ID), data = dm,
               REML = FALSE)

nobs(modela)
nobs(modelb)

anova(modela, modelb, test = "LRT")

``` 

In this case, we can see from the significant p-value that **Model A**
is a significantly better fit to the data than is **Model B**.
Note that LRTs are only appropriate when the models are nested. We
cannot use LRTs for non-nested models.

In nested models, the more complex model is often called the "Full"
model and the simpler model the "Reduced" or "Restricted" version of
the Full model.

## Non Nested Models

Models are considered non nested when one model is not strictly a
constrained version of a more complex model. For example, suppose that 
**Model A** predicts mood from loneliness and stress whereas 
**Model B** predicts mood from loneliness and sex. Written as a formula,
these could be:

$$
ModelA: mood_{ij} = b_{0j} + b1 * loneliness_j + b2 * stress_{ij} + \varepsilon_{ij}
$$

and

$$
ModelB: mood_{ij} = b_{0j} + b1 * loneliness_j + b2 * sex_j + 0 * stress_{ij} + \varepsilon_{ij}
$$

Although **Model B** does have loneliness which also is present in **Model
A** and it has restricted stress to be 0, it has another addition,
sex, which is not in **Model A**. In this case, the two models are
not nested. In this case, although we could fit both models and they
are both on the same number of observations, so the outcome is the
same in both models, the models are not nested. Although we can still
ask `R` to conduct a LRT, this LRT is not valid. It is shown here to
highlight that you as the analyst are responsible for determining
whether your two models are nested or not and therefore deciding
whether a LRT is an appropriate way to evaluate whether the models are
significantly different from each other, or not.

```{r}

modela <- lmer(dMood ~ loneliness + dStress + (1 | ID),
               data = dm, REML = FALSE)
modelb <- lmer(dMood ~ loneliness + sex + (1 | ID), data = dm,
               REML = FALSE)

nobs(modela)
nobs(modelb)

## this is NOT appropriate LRT
anova(modela, modelb, test = "LRT")

``` 

Although we cannot conduct a LRT on non nested models, it *is* still
useful to compare the fit of non-nested models. For example, if one is
a much better fit than another model, that may suggest one set of
predictors is superior or can be used to evaluate competing hypotheses
(e.g., Theory 1 says that stress and family history are the most
important predictors of mood and Theory 2 says that age and sleep are
the best predictors of mood --- these two theories are competing, not
nested versions of each other).

We cannot use LRTs to compare non nested models, but we can use other
measures, including performance measures such as variance explained or
model accuracy and information criterion. We will talk about
information criterion next.

## Information Criterion

Two common information criterion are:

* AIC: the Akaike Information Criterion
  (AIC)^[https://en.wikipedia.org/wiki/Akaike_information_criterion] 
* BIC: the Bayesian Information Criterion
  (BIC)^[https://en.wikipedia.org/wiki/Bayesian_information_criterion]
  also sometimes called the Schwarz Information Criterion

Both the AIC and BIC are calculated based primarily on the log
likelihood, LL, of a model. One way of thinking about these
information criterion is that you could think about models being a way
of approximating reality. Suppose you have two models that both
generate approximations (predictions) of reality. The model whose
predictions are closer to the observed data will have a higher LL.
The LL can be used as a *relative* measure of model fit. LL **is not**
an absolute measure of fit. That is, we do not interpret a specific LL
value as indicating "good" fit. Only which of a set of (all
potentially bad) models is the best.

However, there is a limitation with using LL alone. The LL will always
stay the same or increase as additional predictors / parameters are
added to the model. Thus if we use the LL alone, out of a set of
competing models, we would virtually always pick the more complex
models. To address this, we need to incorporate some penalty for the
complexity of a model, so that to choose a more complex over simpler
model as the "best" it has to improve the LL *enough*.

The AIC and BIC are very similar except that they use different
penalties for model complexity (technically they are derived from
rather different theoretical foundations, but for practical purposes
they are similar other than the complexity penalties).

A common way of defining model complexity is based on the number of
estimated model parameters, $k$. For example, consider the following
simple linear regression for $n$ different observations:

$$
\begin{align}
\eta_{i} &= b_0 + b_1 * x_i\\
y_i &\sim \mathcal{N}(\eta_i, \sigma_{\varepsilon_i})\\
\end{align}
$$

The parameters are:

$$
b_0, b_1, \sigma_{\varepsilon_i}
$$

so $k = 3$.

The equations for AIC and BIC are quite easy to follow and looking at
them helps understand where they are similar and different.

$$
\begin{align}
AIC &= 2 * k - 2 * LL \\
BIC &= log_{e}(n) * k - 2 * LL \\
\end{align}
$$

$n$ is the number of observations included in the model and $LL$ is
the log likelihood, where higher values indicate a better fitting
model. These equations highlight that the only difference between the
AIC and BIC is whether the number of parameters in the model, $k$ are
multiplied by $2$ (AIC) or $log_{e}(n)$. Thus, the AIC and BIC will be
identical if:

$$
\begin{align}
log_{e}(n) &= 2 \\
n &= e^2 \approx 7.39 \\
\end{align}
$$

If $n < e^2 \approx 7.39$ then the BIC will have a weaker penalty
based on the number of parameters, $k$, than the AIC.
If $n > e^2 \approx 7.39$ then the BIC will have a stronger penalty
based on the number of parameters, $k$, than the AIC.
Functionally, a stronger penalty on the number of parameters means
that a particular information criterion will tend to favor more
parsimonious (less complex) models. Thus, for all but the tiniest of
sample sizes ($\approx 7.39$) the BIC will have a stronger penalty on
model complexity and so will favor relatively more parsimonious models
than will AIC.

For both AIC and BIC, the *relatively* better model of those compared
is the model with the lower value (i.e., lower = better for both AIC
and BIC).

**There is no "right" or "wrong" information criterion to use.**
People use both the AIC and/or the BIC. If both AIC and BIC suggest
the same model is the "best" there is no ambiguity. Sometimes the AIC
and BIC disagree regarding which model is the best. In these cases one
must pick which information criterion to go with. Both criterion
require that the number of observations, $n$, be larger than the
number of parameters, $k$ for them to operate well.

A benefit of the AIC and BIC is that both can be used to compare non
nested models and they also can be used to compare nested models. It
is relatively common to use AIC and/or BIC to "choose" amongst a set
of possible models.

In `R` we can usually calculate AIC and BIC using the functions,
`AIC()` and `BIC()`. Here we will calculate the AIC and BIC for our
two models.

```{r}

AIC(modela, modelb)

BIC(modela, modelb)

```

In this case, both the AIC and BIC are lower for **Model A** than for
**Model B**, indicating that **Model A** is the optimal of those two
models. Again, the relative interpretation is important. We cannot
conclude that **Model A** is a "good" model, only that it is better
than **Model B**.

## Model Selection

Integrating what we have learned about LRT for nested models and
AIC/BIC for nested or non nested models, we can compare across several
models to select the best model.

First, we fit a series of models. In all cases, mood is the
outcome. First we have an intercept only model (m0), an unconditional
model, as a baseline reference. This is helpful as all the fit indices
are relative. Then we fit polynomials of stress with degree 1 (linear;
m1); degree 2 (quadratic; m2); degree 3 (cubic; m3). 

We fit different degree polynomials of stress using the `poly()`
function in `R`. Finally we fit a competing model, maybe a linear
effect of energy is a better predictor of mood than stress. 

Next, we check that all the observations are identical across
models. 
**If the observations were not the same across models, we would need
to create a new dataset that had any missing data on any variable used
in any of the models excluded.**  This is a critical step if needed,
because LRTs and AIC and BIC are only valid if based on the same
data.

```{r}

dm[, dStress := as.numeric(dStress)]
dm[, dMood := as.numeric(dMood)]

m0 <- lmer(dMood ~ 1 + (1 | ID), data = dm, REML = FALSE)
m1 <- lmer(dMood ~ poly(dStress, 1) + (1 | ID), data = dm, REML = FALSE)
m2 <- lmer(dMood ~ poly(dStress, 2) + (1 | ID), data = dm, REML = FALSE)
m3 <- lmer(dMood ~ poly(dStress, 3) + (1 | ID), data = dm, REML = FALSE)
malt <- lmer(dMood ~ dEnergy + (1 | ID), data = dm, REML = FALSE)

## check all the observations are the same
nobs(m0)
nobs(m1)
nobs(m2)
nobs(m3)
nobs(malt)

```

Now we can see which model is the best fit. For the nested models
(m0 - m3) we can use LRTs. Technically, m0 also is nested in malt, so
we can use a LRT for that too. We cannot use a LRT for, say, m1 and
malt as those models are not nested. We can use AIC and BIC for all
models, though. Note that with multiple models, the `anova()` function
is doing sequential LRT tests (e.g., m0 vs m1; m1 vs m2, etc.) not all
compared to one model. If you want two specific models compared,
specify just two.

```{r}

#### LRTs for nested models ####

## two specific comparisons
anova(m0, m3, test = "LRT")

## sequential comparisons
anova(m0, m1, m2, m3, test = "LRT")
anova(m0, malt, test = "LRT")

## AIC and BIC for nested and non nested models
AIC(m0, m1, m2, m3, malt)
BIC(m0, m1, m2, m3, malt)

```

From the LRT we can see that m3 is significantly better fit than
m0. Looking at the sequential tests, however, m2 is no better than m1
and m3 is no better than m2, which might suggest m1 is the best
model. If we look at the AIC values, for the models with stress (m1,m2,m3), AIC
selects m1 as the best model (by 1.45 compared to m2, very close). BIC
also selects m1 over m2.
However, both AIC and BIC indicate that the alternate model (malt) is
the best of all the models evaluated, suggesting that linear energy is a
better predictor of mood than is linear, quadratic, or cubic
polynomials of stress.

If we were picking a stress model only, we would pick the
linear (m1) model, because it is best on all indices (LRT, AIC, BIC).
If we were willing to pick energy over stress, we
would clearly go with energy.

In addition to testing fixed effects, these same methods can be used
to test random effects. Here we will fit a new model, m4, that includes
a random slope of stress and the correlation between the random
intercept and stress slope. Then we compare the linear stress with
random slope (m4), linear stress without random slope (m1) and
intercept only model (m0). Please note that m4 did not converge, so 
normally we would simplify our model but we don't want to remove the random
slope of stress for demonstration purposes - so in this case we will have
to remember not to "trust" our results for m4.

```{r}

m4 <- lmer(dMood ~ dStress + (1 + dStress | ID), data = dm, REML = FALSE)

## check all the observations are the same
nobs(m0)
nobs(m1)
nobs(m4)

anova(m0, m1, test = "LRT")

anova(m1, m4, test = "LRT")

anova(m0, m4, test = "LRT")

AIC(m0, m1, m4)
BIC(m0, m1, m4)

```

From these tests we can conclude several things. Firstly adding a
linear trend of stress as a fixed effect is significantly better than
the intercept only model. Secondly (if our model had converged), 
adding a random slope of stress
and correlation of the random slope and intercept is also significantly
better than the fixed linear effect of stress model. Thirdly, the
model with stress as both a fixed and random slope is significantly
better than the intercept only model (m4 vs m0 is kind of an omnibus
test of the effect of adding stress as both a fixed and random slope
into the model).  

Finally, AIC favours m4 and and BIC favours m1 over the other models, 
and since m4 did not converge, suggests that m1 is the best balance of complexity
and model fit (and certainly better than m0).

The main point of this final exercise is to show how you can use LRTs
and AIC and BIC to evaluate whether random effects "help" a model
fit.


# Summary

## Conceptual

Key points to take away conceptually are:

- How to interpret standard notation in LMM equations, including with interaction terms
- What the different types of LMM interactions are
- How to include interactions/moderation in LMMs in R
- How to understand whether there is a significant interaction 
- How to test and interpret interactions 
- How to test simple slopes / simple effects from different kinds of
  interactions 
- How comparing models is a powerful and flexible way to test many
  different hypotheses
- The differences between nested and non nested models, as well as
  what sorts of comparisons are valid for each
- What they are and how to interpret LRT, AIC, BIC
- How to select the "best" model empirically using model comparison techniques.

## Code


| Function       | What it does                                 |
|----------------|----------------------------------------------|
| `lmer()`     | estimate a LMM  |
| `confint()` | calculate confidence intervals for a LMM  | 
| `visreg()` | create marginal or conditional graphs from a LMM  | 
| `modelDiagnostics()` | evaluate model diagnostics for LMMs including of multivariate normality  | 
| `summary()` | get a summary of model results
| `modelTest()` | along with `APAStyler()` get a nicely formatted summary of a model results. | 
| `emmeans()` | test specific means from a model. | 
| `emtrends()` | test simple slopes from a model. | 
| `AIC()` | calculates AIC for LMMs  | 
| `BIC()` | calculates BIC for LMMs  | 
| `anova()` | can be used to compare two nested LMMs and calculate a LRT | 
| `modelTest()` | along with `APAStyler()` get a nicely formatted summary of a model results, including many automatic effect sizes. | 
| `poly()` | Fit a polynomial of specified degree for a predictor in a LMM. | 
 

## Extra (in case you ever want to fix left skews - optional material)

Applying transformations to left skewed data is more difficult as
generally transformations work on long right tails. A solution is to
reverse the variable, apply the transformation and then again reverse
it so that the direction is the same as it originally was. We could
try a square root transformations which is milder than a log
transformation. To reverse it, we subtract the variable from the sum
of its minimum and maximum. Next we take its square root, then we
reverse by again subtracting from the sum of its minimum and maximum,
but square root transformed.

Let's sidetrack a little and try this out with an example using a 
slightly skewed outcome variable, dMood.

```{r}

## First let's see how the model looks like with original dMood scores

mtest1 <- lmer(dMood ~ neuroticism * dStress + (1 | ID), data = dm)

mt1 <- modelDiagnostics(mtest1) 
plot(mt1, nrow = 2, ncol = 2, ask = FALSE)

## Now let's try to fix the very mild left skew just for the sake
## of demonstration

max(dm$dMood) + min(dm$dMood)

## transform
dm[, moodtrans := sqrt(8) - sqrt(8 - dMood)]

mtest <- lmer(moodtrans ~ neuroticism * dStress + (1 | ID), data = dm)

mt <- modelDiagnostics(mtest) 
plot(mt, nrow = 2, ncol = 2, ask = FALSE)

``` 

The transformation appears to have modestly helped the distribution of
residuals. Its not that clear whether it was bad enough to begin with
and whether the transformation improved it enough that it is worth the
difficulty in interpretation (dMood is now square root transformed and
that must be incorporated into its interpretation). For the lecture,
we did this for demonstration purposes, but in practice, consider
whether this is worth it or only adds difficulty in understanding
without improving / changing results much.