Reporting Model Results

modelsummary, broom, and inline reporting

Daniela Palleschi

March 25, 2025

Overview

  • workflow for fitting a linear mixed effect model
  • extracting estimates for inline reporting
  • producing model summary table

Set-up

  • load packages: here, tidyverse, lme4, lmerTest, kableExtra
  • load data: data/processed/data_clean.csv; save it as df_clean
    • or load some data you want to work with
  • load helper function for formatting p-values: source(here("R", "helpers.R"))

Code

Model fitting

  • even with pre-processed data you typically want to do some data wrangling and exploration
  • usually I’d suggest having an Exploratory Data Analysis script before your analysis script
    • here you’d plot the data in all different ways to get a feel for it
    • you’d also want to look at the observations across participants etc.
  • we’ll skip this step

Data prep

  • we need to filter our data
  • apply contrast coding for factorial predictors
  • centre continuous predictors
  • inspect distribution of dependent variable

Fit model

  • we’ll do frequentist, and just run the published model structure
```{r}
#| eval: false
lmer_gaze <-
  lmer(log(gaze) ~ lifetime * tense + trial_c +
         (1 + lifetime | px) +
         (1 + tense | item), 
       data    = df_verb,
       control = lmerControl(optimizer = "bobyqa",        # more stable than default Nelder_Mead
                             optCtrl   = list(maxfun = 1e5)))  # max function evaluations
```

Save and load model

```{r}
#| eval: false
saveRDS(lmer_gaze, here("output", "models", "lmer_gaze.rds"))
```
```{r}
#| eval: true
lmer_gaze <- readRDS(here("output", "models", "lmer_gaze.rds"))
```
  • eval: false: fitting and saving the model
  • eval: true: loading in the saved model

✏️ Exercises — 5 minutes

Exercises 1: Analysis script

You will render a script that already follows these steps.

Reporting your model

  • report model results inline using broom.mixed::tidy()
  • extract the term of interest with filter(), then pull and format each parameter
lmer_tidy        <- broom.mixed::tidy(lmer_gaze, effects = "fixed", conf.int = TRUE)
b_lifetime_lmer  <- lmer_tidy |> filter(term == "lifetime1") |> pull(estimate)  |> round(2)
se_lifetime_lmer <- lmer_tidy |> filter(term == "lifetime1") |> pull(std.error) |> round(2)
t_lifetime_lmer  <- lmer_tidy |> filter(term == "lifetime1") |> pull(statistic) |> round(2)
p_lifetime_lmer  <- lmer_tidy |> filter(term == "lifetime1") |> pull(p.value)   |> fmt_p(3)

Inline reporting

  • embed extracted values in prose using inline `r`
  • results update automatically when the model changes
There was a main effect of lifetime, with longer first-pass reading times for dead 
versus living referents ($\beta$ = 0.05, SE = 0.02, 
*t* = 2.12, *p* < .05).

There was a main effect of lifetime, with longer first-pass reading times for dead versus living referents (\(\beta\) = 0.05, SE = 0.02, t = 2.12, p < .05).