|
Linear mixed models (with Dr. Elizabeth Pankratz) |
Regression refresher, intro to group-structured data |
| Modelling group-structured data using random effects | |
| Interpreting LMMs and building maximal models | |
| Troubleshooting model fit, checking assumptions + diagnostics | |
| LMMs: Practice analysis | |
|
factor analysis working with multi-item measures (with Dr. Josiah King) |
measurement and dimensionality |
| exploring underlying constructs (EFA) | |
| testing theoretical models (CFA) | |
| reliability and validity | |
| recap & exam prep |
Modelling group-structured data using random effects
Data Analysis for Psychology in R 3
Example data: Reaction times in the
Implicit Association Test (IAT)
The IAT asks you to categorise words from different topics together as quickly and accurately as you can.
For example, to test if you implicitly associate anxiety with yourself:
In theory, if you implicitly associate anxiety and yourself, then your reaction times should be faster for the category “anxious and me” (associated) compared to the category, e.g., “anxious and others” (unassociated).
Note: The IAT doesn’t actually work, but the data is structured perfectly to teach you what I want to teach you!
Example data: Log reaction times in the
Implicit Association Test (IAT)
Code
set.seed(1)
implicit_data |>
ggplot(aes(x = pairing, y = logRT)) +
geom_violin() +
geom_jitter(alpha = 0.05, size = 3) +
NULL
implicit_data |>
head(16)# A tibble: 16 x 4
ppt_id item_id pairing logRT
<chr> <chr> <fct> <dbl>
1 ppt1 item1 Associated 4.33
2 ppt1 item1 Unassociated 4.70
3 ppt2 item1 Associated 4.17
4 ppt2 item1 Unassociated 4.91
5 ppt3 item1 Associated 4.45
6 ppt3 item1 Unassociated 5.13
7 ppt4 item1 Associated 4.41
8 ppt4 item1 Unassociated 6.54
9 ppt5 item1 Associated 5.19
10 ppt5 item1 Unassociated 6.64
11 ppt6 item1 Associated 4.74
12 ppt6 item1 Unassociated 5.96
13 ppt7 item1 Associated 5.23
14 ppt7 item1 Unassociated 5.70
15 ppt8 item1 Associated 4.23
16 ppt8 item1 Unassociated 6.19Side note: Why are we log-transforming
reaction times?
Raw reaction times:
Code
implicit_data |>
mutate(RT = exp(logRT)) |>
ggplot(aes(x = pairing, y = RT)) +
geom_violin() +
geom_jitter(alpha = 0.05, size = 3) +
labs(
y = 'RT in milliseconds',
) +
NULL
Log-transformed reaction times:
Code
set.seed(1)
implicit_data |>
ggplot(aes(x = pairing, y = logRT)) +
geom_violin() +
geom_jitter(alpha = 0.05, size = 3) +
NULL
Reaction times in milliseconds have a very skewed distribution, so the assumption of normal distribution of errors is not met. To get around this, we take their logarithm.
Plot the data for a few individual participants
Code
ppts_to_plot <- c(56, 14, 33, 24, 44, 31, 67, 97, 99)
p_implicit_ppts <- implicit_data |>
mutate(
pairing_num = ifelse(pairing == 'Unassociated', 0, 1)
) |>
filter(ppt_id %in% paste0('ppt', ppts_to_plot)) |>
ggplot(aes(x = pairing_num, y = logRT)) +
geom_jitter(alpha = 0.2, width = 0.1, size = 3) +
stat_summary(fun = mean, geom = 'point', size = 5) +
stat_summary(aes(group = ppt_id), fun = mean, geom = 'line') +
facet_wrap(~ ppt_id) +
scale_x_continuous(breaks = c(0, 1), limits = c(-0.3, 1.3), labels = c('Unassoc.', 'Assoc.')) +
labs(x = 'Topic pairing') +
theme(
strip.background = element_blank(),
panel.border = element_rect(linewidth = 1),
panel.grid = element_blank()
)
p_implicit_ppts
Plot the data for a few individual items
Code
items_to_plot <- c(1, 2, 4, 13, 6, 8, 14, 27, 32)
p_implicit_items <- implicit_data |>
mutate(
pairing_num = ifelse(pairing == 'Unassociated', 0, 1)
) |>
filter(item_id %in% paste0('item', items_to_plot)) |>
ggplot(aes(x = pairing_num, y = logRT)) +
geom_jitter(alpha = 0.1, width = 0.1, size = 3) +
stat_summary(fun = mean, geom = 'point', size = 5) +
stat_summary(aes(group = item_id), fun = mean, geom = 'line') +
facet_wrap(~ item_id) +
scale_x_continuous(breaks = c(0, 1), limits = c(-0.3, 1.3), labels = c('Unassoc.', 'Assoc.')) +
labs(x = 'Topic pairing') +
theme(
strip.background = element_blank(),
panel.border = element_rect(linewidth = 1),
panel.grid = element_blank()
)
p_implicit_items
The simple LM has no way of knowing that different participants behave differently
Code
p_implicit_ppts +
geom_abline(
intercept = coef(implicit_lm)[['(Intercept)']],
slope = coef(implicit_lm)[['pairingAssociated']],
colour = 'red',
linewidth = 1
)
The red line is the line estimated by the simple LM: intercept = 4.77, slope = –0.76.
Now we have a different line for each participant
Code
indiv_ppt_coefs <- coef(implicit_lmlist) |>
rownames_to_column(var = 'ppt_id')
p_implicit_ppts +
geom_abline(
data = filter(indiv_ppt_coefs, ppt_id %in% paste0('ppt', ppts_to_plot)),
aes(intercept = `(Intercept)`, slope = pairingAssociated, colour = ppt_id),
linewidth = 1
) +
theme(legend.position = 'none')
… but these models are still averaging over all the different items and ignoring their variability.
Intercept adjustments by ppt_id
(1 | ppt_id) tells the model to estimate how much the average intercept should be nudged up or down, in order to best fit the data from each participant.
We can use ranef() to extract the list of these adjustments:
ranef(implicit_lmm_int)$ppt_id |> head() (Intercept)
ppt1 -0.1190
ppt10 0.2882
ppt100 -0.0651
ppt11 0.9000
ppt12 0.1716
ppt13 -0.4369For example, for ppt1:
For ppt10:
For ppt100:
Plot each participant’s intercept adjustments
dotplot.ranef.mer(
ranef(implicit_lmm_int)
)$ppt_id
How well does the random-intercept model fit each participant’s data?
Code
random_int_coefs <- coef(implicit_lmm_int)$ppt_id |>
rownames_to_column(var = 'ppt_id')
p_implicit_ppts +
geom_abline(
data = filter(random_int_coefs, ppt_id %in% paste0('ppt', ppts_to_plot)),
aes(intercept = `(Intercept)`, slope = pairingAssociated, colour = ppt_id),
linewidth = 1
) +
theme(legend.position = 'none')
Intercept and slope adjustments
When we include (1 + pairing | ppt_id) in the model formula, we tell the model to estimate how much the fixed intercept AND the fixed slope over pairing should be nudged up or down by in order to fit the data from each participant.
Let’s look at the adjustments that the model has estimated:
ranef(implicit_lmm_int_slp)$ppt_id |> head() (Intercept) pairingAssociated
ppt1 -0.171 0.0692
ppt10 0.307 -0.0769
ppt100 -0.120 0.0599
ppt11 1.117 -0.3740
ppt12 0.242 -0.0957
ppt13 -0.523 0.1651For example, for ppt1:
For ppt10:
Plot each participant’s intercept and slope adjustments
dotplot.ranef.mer(
ranef(implicit_lmm_int_slp)
)$ppt_id
How well does the random-intercept and random-slope model fit each participant’s data?
Code
random_int_slp_coefs <- coef(implicit_lmm_int_slp)$ppt_id |>
rownames_to_column(var = 'ppt_id')
p_implicit_ppts +
geom_abline(
data = filter(random_int_slp_coefs, ppt_id %in% paste0('ppt', ppts_to_plot)),
aes(intercept = `(Intercept)`, slope = pairingAssociated, colour = ppt_id),
linewidth = 1
) +
theme(legend.position = 'none')
It’s a lot better than the random-intercept-only model!
But still not awesome, because there’s still another source of variability that we haven’t yet modelled…
No adjustments for each item :(
Code
p_implicit_items +
geom_abline(
intercept = fixef(implicit_lmm_int_slp)[['(Intercept)']],
slope = fixef(implicit_lmm_int_slp)[['pairingAssociated']],
colour = 'red',
linewidth = 1
) 
Currently, the model currently averages over all items and imagines that they all pattern exactly the same, because we haven’t told it that each item might also vary.
The red line is the line with the fixed effect parameters: intercept = 4.77, slope = –0.76.
Look at all the adjustments
ranef(implicit_full_lmm)$ppt_id |>
head(12) (Intercept) pairingAssociated
ppt1 -0.229 0.2028
ppt10 0.250 0.0758
ppt100 -0.199 0.2488
ppt11 1.229 -0.5998
ppt12 0.316 -0.2672
ppt13 -0.542 0.1880
ppt14 -2.317 0.8716
ppt15 0.625 -0.1650
ppt16 0.578 -0.0507
ppt17 -0.454 0.1080
ppt18 0.254 -0.2137
ppt19 0.187 -0.2124For example, for ppt1:
ranef(implicit_full_lmm)$item_id |>
head(12) (Intercept) pairingAssociated
item1 0.6479 -0.392
item10 0.2301 0.189
item11 0.3806 0.111
item12 0.3166 -0.147
item13 1.1422 -0.274
item14 -0.8104 0.310
item15 0.0427 0.131
item16 0.4579 -0.178
item17 0.4274 -0.343
item18 -1.2246 0.124
item19 -1.2392 0.477
item2 1.6616 -0.689For example, for item1:
Plot the by-participant intercept and slope adjustments
dotplot.ranef.mer(ranef(implicit_full_lmm))$ppt_id
Plot the by-item intercept and slope adjustments
dotplot.ranef.mer(ranef(implicit_full_lmm))$item_id
This model captures by-participant variability
Code
full_llm_ppt_coefs <- coef(implicit_full_lmm)$ppt_id |>
rownames_to_column(var = 'ppt_id')
p_implicit_ppts +
geom_abline(
data = filter(full_llm_ppt_coefs, ppt_id %in% paste0('ppt', ppts_to_plot)),
aes(intercept = `(Intercept)`, slope = pairingAssociated, colour = ppt_id),
linewidth = 1
) +
theme(legend.position = 'none')
This model also captures by-item variability!
Code
full_llm_item_coefs <- coef(implicit_full_lmm)$item_id |>
rownames_to_column(var = 'item_id')
p_implicit_items +
geom_abline(
data = filter(full_llm_item_coefs, item_id %in% paste0('item', items_to_plot)),
aes(intercept = `(Intercept)`, slope = pairingAssociated, colour = item_id),
linewidth = 1
) +
theme(legend.position = 'none')
To do this week
Tasks:
Work on exercises in labs
Complete the weekly quiz
Get support:
Consult the flash cards
Ask questions anonymously on Piazza
We really like seeing you in office hours!
Wooclap images (1)
Is the participant’s intercept more positive or more negative than the average intercept?
Wooclap images (2)
Is the participant’s slope more positive or more negative than the average slope?
Wooclap images (3)
Are the participant’s intercept and slope each more positive or more negative than the average intercept and slope?
