---
title: "Chapter 14: Multilevel models"
author: "Andy Field"
date: "2025-05-14"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
    code-overflow: wrap
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(glmmTMB)
library(tidyverse)


# load data
cosmetic_tib <- here::here("data/cosmetic.csv") |>
  read_csv() |>
  mutate(
    clinic = as_factor(clinic),
    reason = as_factor(reason)
    )

# or load from the discovr package
# cosmetic_tib <- discovr::cosmetic

```

## Data preparation

This code reads the file in and converts the variables **clinic** and **reason** to factors. It's a good idea to check that the levels of **clinic** are sequential from clinic 1 to clinic 21, and that **reason** has levels ordered change appearance and physical reason. Check the factor levels by executing:

```{r}
#| label: chk_levels
#| 
levels(cosmetic_tib$clinic)
levels(cosmetic_tib$reason)
```

If they're not in the correct order then:

```{r}
#| label: relevel_factors
#| 
cosmetic_tib <- cosmetic_tib |>
  mutate(
    clinic = fct_relevel(clinic, paste("Clinic", seq(1, 20, 1))),
    reason = fct_relevel(reason, "Change appearance")
    )
```

View data for screenshot

```{r}
#| label: view_sample

head(cosmetic_tib, n = 10)
```


## Plot the data

```{r}
#| label: plot_clinics

ggplot(cosmetic_tib, aes(days, post_qol)) +
  geom_point(position = position_jitter(width = 0.1, height = NULL), alpha = 0.5, size = 1, colour = "#999933") +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.75, colour = "#CC6677") +
  coord_cartesian(xlim = c(0, 400), ylim = c(0, 100)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(x = "Days post surgery", y = "Quality of life after surgery (%)", colour = "Clinic") +
  facet_wrap(~ clinic, ncol = 4) +
  theme_minimal()
```

## Summary statistics

```{r}
#| label: summary_by_clinic

qol_sum <- cosmetic_tib |>
  group_by(clinic) |>
  describe_distribution(select = c("base_qol", "post_qol"), ci = 0.95)

qol_sum |> 
  display()
```


## Rescale predictor

```{r}
#| label: rescale_days

cosmetic_tib <- cosmetic_tib |>
  mutate(
    months = days*12/365
    )
```



## Fixed effect fits

### Fit the model to the pooled data

```{r}
#| label: pooled_model
#| 
pooled_lm <- lm(post_qol ~ months*reason + base_qol, data = cosmetic_tib)
model_parameters(pooled_lm) |> 
  display()
```

### Fit the model across contexts

```{r}
#| label: models_by_clinic
 
clinic_lms <- cosmetic_tib |>
  arrange(clinic) |>
  group_by(clinic) |>
  nest() |>
  mutate(
    model = purrr::map(.x = data,
                       .f = \(clinic_tib) lm(post_qol ~ months*reason + base_qol, data = clinic_tib)),
    coefs = purrr::map(model, model_parameters)
    )

models <- clinic_lms |>
  select(-c(data, model)) |>
  unnest(coefs)

# display the model for a particular clinic
models |> 
  filter(clinic == "Clinic 9") |> 
  display()


# plot the parameters
ggplot(data = models, aes(Coefficient)) +
  geom_density() +
  facet_wrap(~Parameter , scales = "free") +
  theme_minimal()
```

## Random effect models

### Fit the model

```{r}
#| label: build_model

# optionally build the model in steps
intcpt_mlm <- glmmTMB::glmmTMB(post_qol ~ 1 + (1|clinic), data = cosmetic_tib) # create intercept only model
months_mlm <- update(intcpt_mlm, .~. + months) # add months
monthsre_mlm <- update(months_mlm, .~  months + (1 + months|clinic)) # add random slope for months
reason_mlm <- update(monthsre_mlm, .~.  + reason) # add reason
qol_mlm <- update(reason_mlm, .~.  + base_qol) # add baseline QoL
cosmetic_mlm <- update(qol_mlm, .~.  + months:reason) # add the interaction term

# Compare the models

test_lrt(intcpt_mlm, months_mlm, monthsre_mlm, reason_mlm, qol_mlm, cosmetic_mlm) |> 
  display()
```


Alternatively, because we have a theory driven model we can fit only the final one:

```{r}
#| label: fit_model
#| eval: false

cosmetic_mlm <- glmmTMB::glmmTMB(post_qol ~ months*reason + base_qol + (months|clinic), data = cosmetic_tib)
```


```{r}
#| label: instect_mlm

# Optional F-tests of effects:

car::Anova(cosmetic_mlm, type = 3) |> 
  model_parameters() |> 
  display(use_symbols = TRUE)

# inspect the model

model_performance(cosmetic_mlm) |> 
  display()
model_parameters(cosmetic_mlm) |> 
  display()
```


### Standardizing variables

We don't use this model, but if you wanted to fit the model to standardized variables then use this code.

```{r}
#| label: z_mlm
#| eval: false

cosmetic_z <- cosmetic_tib |>
  standardize()

cosmetic_modz <- glmmTMB::glmmTMB(post_qol ~ months*reason + base_qol + (months|clinic), data = cosmetic_z)
model_parameters(cosmetic_modz) |> 
  display()
```


### Simple slopes

```{r}
#| label: simple_slopes
#| 
cosmetic_slopes <- estimate_slopes(cosmetic_mlm,
                                               trend = "months",
                                               by = "reason",
                                               ci = 0.95) 
cosmetic_slopes |> 
  display()
```

Plot the simple slopes

```{r}
#| label: plot_slopes
#| 
estimate_means(cosmetic_mlm, c("months", "reason")) |> 
  plot() +
  labs(x = " Months since surgery", y = "Quality of life post-surgery (0-100)", colour = "Reason for surgery", fill = "Reason for surgery") +
  theme_minimal()
```

## Diagnostic checks

```{r}
#| label: mlm_diag

check_model(cosmetic_mlm)

```



