R code Chapter 14

Load packages

Remember to load the tidyverse:

library(tidyverse)

and the other packages we’ll use

library(broom)
library(broom.mixed)
library(datawizard)
library(lme4)
library(lmerTest)

Load the data

Load the data from the discover package:

cosmetic_tib <- discovr::cosmetic

If you want to read the file from the CSV and you have set up your project folder as I suggest in Chapter 1, then the code you would use is:

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

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:

levels(cosmetic_tib$clinic)

##  [1] "Clinic 1"  "Clinic 2"  "Clinic 3"  "Clinic 4"  "Clinic 5"  "Clinic 6" 
##  [7] "Clinic 7"  "Clinic 8"  "Clinic 9"  "Clinic 10" "Clinic 11" "Clinic 12"
## [13] "Clinic 13" "Clinic 14" "Clinic 15" "Clinic 16" "Clinic 17" "Clinic 18"
## [19] "Clinic 19" "Clinic 20"

levels(cosmetic_tib$reason)

## [1] "Change appearance" "Physical reason"

If they’re not in the correct order then:

cosmetic_tib <- cosmetic_tib |>
  dplyr::mutate(
    clinic = forcats::fct_relevel(clinic, paste("Clinic", seq(1, 20, 1))),
    reason = forcats::fct_relevel(reason, "Change appearance")
    )

Plot the data

Create the first plot in the chapter

ggplot2::ggplot(cosmetic_tib, aes(days, post_qol, colour = clinic)) +
  geom_point(alpha = 0.5, size = 1) +
  geom_smooth(method = "lm", se = FALSE, size = 0.5) +
  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") +
  guides(colour=guide_legend(ncol=2)) +
  theme_minimal()

Now the second plot

ggplot2::ggplot(cosmetic_tib, aes(days, post_qol)) +
  geom_point(position = position_jitter(width = 0.1, height = NULL), alpha = 0.5, size = 1) +
  geom_smooth(method = "lm", se = FALSE, size = 0.5) +
  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

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

qol_sum |>
  dplyr::mutate(
    Clinic = stringr::str_remove(.group, "clinic=")
    ) |>
  dplyr::select(Clinic, Variable:n) |>
  knitr::kable(caption = "Summary statistics for the surgery data",
               digits = 2)

Table 1: Summary statistics for the surgery data

Fixed effect fits

Fit the model to the pooled data

pooled_lm <- lm(post_qol ~ days*reason + base_qol, data = cosmetic_tib)

broom::tidy(pooled_lm, conf.int = TRUE) |>
  knitr::kable(caption = "Parameter estimates for the pooled data model",
               digits = 3)

Table 2: Parameter estimates for the pooled data model

Fit the model across contexts

clinic_lms <- cosmetic_tib |>
  dplyr::arrange(clinic) |>
  dplyr::group_by(clinic) |>
  tidyr::nest() |>
  dplyr::mutate(
    model = purrr::map(.x = data,
                       .f = \(clinic_tib) lm(post_qol ~ days*reason + base_qol, data = clinic_tib)),
    coefs = purrr::map(model, tidy, conf.int = TRUE)
    )

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

models |>
  ggplot(aes(estimate)) +
  geom_density() +
  facet_wrap(~term , scales = "free") +
  theme_minimal()

Random effect models

Try to fit the model

cosmetic_mod <- lmerTest::lmer(post_qol ~ days*reason + base_qol + (days|clinic), data = cosmetic_tib)

allFit(cosmetic_mod)

## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]

## original model:
## post_qol ~ days * reason + base_qol + (days | clinic) 
## data:  cosmetic_tib 
## optimizers (7): bobyqa, Nelder_Mead, nlminbwrap, nmkbw, optimx.L-BFGS-B,nloptwrap.NLOPT_LN_N...
## differences in negative log-likelihoods:
## max= 0.304 ; std dev= 0.12

Rescale days

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

Standardizing variables

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

cosmetic_tib <- cosmetic_tib |>
  dplyr::mutate(across(.cols = where(is.double),
                       .fns = \(column) (column-mean(column, na.rm = T))/sd(column, na.rm = T),
                       .names = "z_{.col}")
                )

cosmetic_modz <- lmerTest::lmer(z_post_qol ~ z_days*reason + z_base_qol + (z_days|clinic), data = cosmetic_tib)

broom.mixed::tidy(cosmetic_modz)

Fit the model using months

cosmetic_mod <- lmerTest::lmer(post_qol ~ months*reason + base_qol + (months|clinic), data = cosmetic_tib, REML = T)

allFit(cosmetic_mod)

## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]

## original model:
## post_qol ~ months * reason + base_qol + (months | clinic) 
## data:  cosmetic_tib 
## optimizers (7): bobyqa, Nelder_Mead, nlminbwrap, nmkbw, optimx.L-BFGS-B,nloptwrap.NLOPT_LN_N...
## differences in negative log-likelihoods:
## max= 5.3e-07 ; std dev= 1.98e-07

Change the optimizer

cosmetic_bob <- lmerTest::lmer(
  post_qol ~ months*reason + base_qol + (months|clinic),
  data = cosmetic_tib,
  control = lmerControl(optimizer="bobyqa")
  )

Use maximum likelihood estimation

cosmetic_ml <- lmerTest::lmer(
  post_qol ~ months*reason + base_qol + (months|clinic),
  data = cosmetic_tib,
  REML = F
  )

Interpreting the model

anova(cosmetic_bob) |> 
  knitr::kable(caption = "Fixed effects for the final model",
               digits = c(rep(2, 5), 3))

Table 3: Fixed effects for the final model

options(knitr.kable.NA = '')
broom.mixed::tidy(cosmetic_bob, conf.int = T) |> 
  knitr::kable(caption = "Parameter estimates for the final model",
               digits = c(rep(2, 7), 3, 2, 2))

Table 4: Parameter estimates for the final model

cosmetic_slopes <- modelbased::estimate_slopes(cosmetic_bob,
                                               trend = "months",
                                               at = "reason",
                                               ci = 0.95) 

cosmetic_slopes |> 
  knitr::kable(caption = "Change over time for different reasons for surgery",
        digits = 3)

Table 5: Change over time for different reasons for surgery

Some alternative ways to get the same simple slopes

cosmetic_slopes <- interactions::sim_slopes(
  cosmetic_bob,
  pred = months,
  modx = reason,
  confint = TRUE
  ) 

cosmetic_slopes <-  emmeans::emtrends(cosmetic_bob, var = "months", spec = "reason")

Plot the simple slopes

interactions::interact_plot(
  model = cosmetic_bob,
  pred = months,
  modx = reason,
  interval = TRUE,
  x.label = "Months since surgery",
  y.label = "Quality of life post-surgery (0-100)",
  legend.main = "Reason for surgery"
  )