Mixed Effects activity in R, based on Brown (2021)
Jamil Palacios Bhanji - Statistical Methods
2024-04-10
Today’s activity is based entirely on Violet Brown’s 2021 article, “An Introduction to Linear Mixed-Effects Modeling in R.” Advances in Methods and Practices in Psychological Science, 4(1), pp. 1–19.
Goals for today
- Learn differences between a linear mixed effects model approach and
a standard linear model approach with aggregated data for repeated
measures in a crossed design
- Learn about fixed and random effects and how to specify them in a
formula using the
lme4::lmer()function in R
- Learn how to address problems in model fitting (failure to converge,
or singular fit)
- Learn how to interpret output from a linear mixed effects
model
- Use model comparison for significance tests (for full model compared to null model, for individual terms in a model)
Step 0
- install the required packages:
install.packages("optimx")*note: some versions of R may make it tricky to install this package - it is only needed for one step so you can skip it if needed
lme4you should already have installed
afexyou should already have installed
- setup your project folder (with data and r_docs folders within) and
start a new project in this existing folder
- start a new R markdown for your code and notes and include the usual
code in the setup chunk, but add
library(lme4)to the setup chunk (and then run the chunk)- Older versions of R may give an error along the lines of
cholmod_factor_ldetA not provided by package Matrixwhen you try to load the package. If this happens, trying installing theMatrixandlme4packages directly from CRAN by running these commands in the console (I recommend that you restart RStudio first - don’t save the workspace - Thanks to Amanda Derrell for looking this up!):
oo <- options(repos = "https://cran.r-project.org/")changes the package repository setting and saves original setting inoo
install.packages("Matrix")installs from CRAN
install.packages("lme4")installs from CRAN
options(oo)restores repository setting to what you had before
- Older versions of R may give an error along the lines of
- Download these files and place in the data folder within your
project folder:
- rt_dummy_data.csv for example 1
Example 1: trial level analysis of a within-subject (2 level) effect
Step 1 - import the first data file and check it out
Data description:
This data set is adapted from a study investigating comprehension of
spoken words based on audio alone (Audio-only condition) vs audio with
video of the speaker (Audiovisual condition). On each of 553 trials,
participants heard and repeated a single word (either Audio-only or
Audiovisual) while simultaneously performing a second task (judging the
duration of a vibrating stimulus on their knee). The researchers
hypothesized that response times for the secondary task would be longer
in the Audiovisual condition than the Audio-only condition, which would
indicate greater listening effort required in the Audiovisual condition.
This finding is relevant to theories of speech processing.
Variables:
- Dependent variable: RT=response time
(ms) for the secondary task
- Independent variable: modality=listening
condition (“Audio-only” or “Audiovisual”)
- stim=word stimulus for each trial
- PID=participant identifier
Each of 53 participants complete 553 trials. The response times have been modified as described in Brown (2021) to illustrate issues involved in linear mixed effects modeling.
What to do:
- Use the
read_csv()function to import thert_dummy_data.csvfile and store it in a variable namedrt_data. There are no missing data that you need to deal with, because trials with a missed response are not included in the data file. Store themodalityvariable as a factor.
#import the data
rt_data <- readr::read_csv(
"data/rt_dummy_data.csv", show_col_types = FALSE) |>
mutate(
modality = forcats::as_factor(modality)
)
- Use the
group_by()andsummarise()functions to compute the following and print it to the screen - notice that you are ignoring thePIDwhen computing these values (not something you would normally do but we do it here to aid learning):- mean RT across all trials in each modality condition
- median RT across all trials in each modality condition
- standard devation of RT across all trials in each modality condition
- mean RT across all trials in each modality condition
#descriptives across trials (ignoring participant ID)
rt_data |> dplyr::group_by(modality) |> dplyr::summarise(
totalN = n(),
trialwiseRTmean = mean(RT, na.rm = TRUE),
trialwiseRTmedian = median(RT, na.rm = TRUE),
trialwiseRTsd = sd(RT, na.rm = TRUE)
) |> ungroup() |>
kableExtra::kbl(
caption = "descriptives by modality condition (trialwise, PID ignored)",
digits = 2) |>
kableExtra::kable_classic(full_width=FALSE)| modality | totalN | trialwiseRTmean | trialwiseRTmedian | trialwiseRTsd |
|---|---|---|---|---|
| Audio-only | 10647 | 1041.46 | 1013 | 302.42 |
| Audiovisual | 11032 | 1124.52 | 1100 | 316.56 |
- Now, first aggregate the data by participant (calculate the number
of trials each participant did of each condition, and the mean RT for
each condition for each subject) and then print out:
- mean number of trials in each condition for a participant
- min and max number of trials in each condition for a
participant
- mean RT across participants (different than before bc you are first
aggregating by participant)
- standard deviation of RT across participants
#first, calc num trials and mean RT for each participant (by condition)
rt_bysub <- rt_data |> dplyr::group_by(modality,PID) |> dplyr::summarise(
bysub_N = n(),
bysub_RT = mean(RT),
bysub_medRT = median(RT),
) |> ungroup()## `summarise()` has grouped output by 'modality'. You can override using the
## `.groups` argument.# now print descriptives, calculated by first aggregating by subject
rt_bysub |> dplyr::group_by(modality) |> dplyr::summarise(
bysub_mean_N = mean(bysub_N),
bysub_min_N = min(bysub_N),
bysub_max_N = max(bysub_N),
bysub_mean_RT = mean(bysub_RT),
bysub_sd_RT = sd(bysub_RT)
) |> ungroup() |>
kableExtra::kbl(
caption = "descriptives, first aggregated by participant",
digits = 2) |>
kableExtra::kable_classic(full_width=FALSE)| modality | bysub_mean_N | bysub_min_N | bysub_max_N | bysub_mean_RT | bysub_sd_RT |
|---|---|---|---|---|---|
| Audio-only | 200.89 | 141 | 244 | 1044.07 | 169.89 |
| Audiovisual | 208.15 | 132 | 264 | 1127.23 | 177.91 |
- Now make a histogram and box plot (split by condition) of the
trialwise RTs (ignoring
PID), then do the same for the RTs that were first aggregated by participant.
# histogram and box plot of the trialwise data
rt_data |> ggplot( aes(x=RT, fill=modality) ) +
geom_histogram(position = "identity", alpha = .5, binwidth = 100) +
theme_classic() + labs(title = "trialwise RTs - histogram")
rt_data |> ggplot( aes(x=modality, y=RT) ) +
geom_boxplot() + theme_classic() + labs(title = "trialwise RTs - box plot")
# histogram and box plot of the by participant data
rt_bysub |> ggplot( aes(x=bysub_RT, fill=modality) ) +
geom_histogram(position = "identity", alpha = .5, binwidth = 100) +
theme_classic() + labs(title = "RTs averaged by participant")
rt_bysub |> ggplot( aes(x=modality, y=bysub_RT) ) +
geom_boxplot() + theme_classic() + labs(title = "RTs averaged by participant - box plot")
- in the trialwise data averaged by condition (ignoring
PID) you should see that mean RT in the Audio-only condition is 1041 ms and in the Audiovisual condition it is 1125 ms
- in the data that is first averaged by participant and then across participants, the mean RT in the Audio-only condition is 1044 ms and in the Audiovisual condition it is 1127 ms
- the difference in means calculated in the two different ways is due
to different numbers of trials for each participant and condition in the
trialwise calculation (look at the min and max trial numbers), but in
the second calculation each participant’s average RTs are weighted
equally. The difference in the two calculations is small here but can be
larger in datasets with larger imbalances.
- You should also see that the data are positively skewed (whether we look at trialwise data or aggregated by subject, this is common with RT data) - but we have a sufficient sample size so we aren’t concerned by the non-normal distribution.
Step 2: Use a traditional method to analyze aggregated RTs (by participant or by item)
Step 2.1 - Model the data after aggregating by participant
First, let’s use each participant’s average RT per condition as our
unit of analysis. Once we average across trials to get a two values
(Audio-only mean RT, Audio-visual mean RT) for each participant, then we
can just do a paired t-test on those values like we have done
before.
Try it now, you can use the tibble of RTs averaged by subject and
modality that you created in the step above (called
rt_bysub in the solution code), and pipe it to the
t.test() function that we’ve used before.
# paired t-test with mean RT per subject in each condition
rt_bysub |> dplyr::arrange(PID) |> t.test(bysub_RT ~ modality, data = _, paired=TRUE)##
## Paired t-test
##
## data: bysub_RT by modality
## t = -6.6122, df = 52, p-value = 2.055e-08
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -108.39746 -57.92312
## sample estimates:
## mean difference
## -83.16029# equivalently, you can run the same test with the lmer function. the formula will make
# more sense after doing the whole activity
rt_bysub_lmer <- rt_bysub |> lme4::lmer(formula = bysub_RT ~ modality + 1 + (1|PID), data = _)
summary(rt_bysub_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: bysub_RT ~ modality + 1 + (1 | PID)
## Data: rt_bysub
##
## REML criterion at convergence: 1305.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.59477 -0.46495 -0.00148 0.37616 2.64062
##
## Random effects:
## Groups Name Variance Std.Dev.
## PID (Intercept) 26066 161.45
## Residual 4192 64.74
## Number of obs: 106, groups: PID, 53
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1044.07 23.89 43.697
## modalityAudiovisual 83.16 12.58 6.612
##
## Correlation of Fixed Effects:
## (Intr)
## modltyAdvsl -0.263
- The mean response time in the Audiovisual condition is 83ms slower
than in the Audio-only condition. The t-statistic (t(52) = -6.612) and
low p-value (.00000002055) show that it is unlikely to observe sample
means this different if there is no true difference.
- This is a simple and commonly used approach to analyze data with
many trials per condition, where the first step is to average over all
trials in a condition for each participant (referred to as
by-participant analysis in Judd, Westfall, and Kenny
(2012)).
- today we’re going over an alternate approach, linear mixed modeling,
where the unit of analysis can be individual trials. Some advantages of
this approach are explained in Brown (2021) and
Judd, Westfall, and Kenny
(2012):
- captures variance (e.g. trial-to-trial) that may be important for
accurately estimating effects of interest
- allows modeling of multiple random variables (e.g., participants and
stimuli), whereas ANOVA approaches allow only one random variable
- handles unbalanced designs/missed trials
- handles categorical or continuous predictors, coefficients give
magnitude and direction of effects (as opposed to ANOVA approach)
- extends to other types of outcome variables (e.g., binary)
- captures variance (e.g. trial-to-trial) that may be important for
accurately estimating effects of interest
Next, we will use a linear mixed model approach where trial RTs are the outcome variable, and the model allows for varying intercepts for each participant (e.g. some participants may tend to respond slower overall) and varying intercepts for each word stimulus (e.g., some words may tend to require longer response times) - we refer to these random effects as random intercepts by participant and by stimulus. The model also allows for a varying effect of modality for different subjects and for different stimuli (random slopes by participant and by stimulus).
Step 3: Linear mixed-effects model of the effect of modality on trial response time
Reminder: what are fixed and random effects?
modalityis systematically varied so we consider it a fixed effect. In designs (e.g., correlation) where no variable is systematically manipulated we can consider a fixed effect to be one where all levels of interest are covered (e.g.,agein the lumosity data we used long ago). We expect fixed effects to be consistent across experiments.
- If the levels of a variable are a sample of a larger set we can
consider it a random effect (in this activity: participants sampled from
a population, and a set of words sampled from a language)
- participants are a sample of a larger population so we consider
PIDa random effect
- word stimuli used in the study are a sample from a language so we
consider
stima random effect
Barr et al (2013) advise using the maximal random effects structure justified by the design. In this case we are justified to include the full random effect structure (random intercepts and slopes for modality, by participant and by stimulus).
Full random effect structure
Here is how we include random effects in our formula specification
(see
Martin Speekenbrink’s book for a more complete explanation -
screenshot at the bottom of this guide):
- by-participant varying intercepts: (1|PID)
- by-stimulus varying intercepts: (1|stim)
- by-participant varying effect of modality, aka random slope
(modality|PID)
- by-stimulus varying effect of modality, aka random slope
(modality|stim)
So let’s specify the model in a typical R formula syntax. Just the
fixed effects part would be RT ~ 1 + modality to specify
that RT is the DV and it is predicted by
modality plus an intercept term. Remember that
lm() will automatically dummy code (aka treatment code) a
factor variable. The function for linear mixed models,
lmer(), does the same. So entering modality as
a predictor is equivalent to a variable where “Audio-only” is set to 0
and “Audiovisual” is set to 1.
So we start with the formula RT ~ 1 + modality - Then we
add the by-participant varying intercepts like this:
RT ~ 1 + modality + (1 |PID)
- add the by-stimulus varying intercepts:
RT ~ 1 + modality + (1|PID) + (1|stim)
- add the by-participant varying effect of modality (random slope by
participant):
RT ~ 1 + modality + (1 + modality|PID) + (1|stim)
- add the by-stimulus varying effect of modality (random slope by
stimulus):
RT ~ 1 + modality + (1 + modality|PID) + (1 + modality|stim)
Step 3.1 - Fit the model
Okay, we have the model formula. We will use the
lme4::lmer() function to estimate the model - you call it
in basically the same way that we have used the lm()
function before. Pass the data and the formula to the function and store
the result in a variable (let’s call it rt_full_mod). Then
pass that variable to the summary() function.
rt_full_mod <- lmer(RT ~ 1 + modality + (1 + modality|PID) +
(1 + modality|stim), data = rt_data)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00220168 (tol = 0.002, component 1)summary(rt_full_mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: RT ~ 1 + modality + (1 + modality | PID) + (1 + modality | stim)
## Data: rt_data
##
## REML criterion at convergence: 302385.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3646 -0.6964 -0.0140 0.5886 5.0003
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## stim (Intercept) 304.0 17.44
## modalityAudiovisual 216.8 14.73 0.16
## PID (Intercept) 28551.1 168.97
## modalityAudiovisual 7708.7 87.80 -0.17
## Residual 65258.9 255.46
## Number of obs: 21679, groups: stim, 543; PID, 53
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1044.14 23.36 44.706
## modalityAudiovisual 83.18 12.57 6.615
##
## Correlation of Fixed Effects:
## (Intr)
## modltyAdvsl -0.178
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00220168 (tol = 0.002, component 1)
Did you get a message that the “model failed to converge”? The
message means that the algorithm used to estimate the model parameters
could not find a good fit within the allotted number of iterations.
Although the summary() function gives reasonable looking
parameters in this case, you should not report results of a
model that failed to converge - the nonconvergence means the
model has not been reliably estimated.
What to do when a model does not converge
Here are strategies covered in Brown
(2021):
1. Recheck your model and make sure it is not misspecified
2. Are there major imbalances in the data? (e.g., if one participant or
item has very few observations it can cause nonconvergence - consider
dropping the participant/item)
3. Adjust the parameters that control the model fitting process (e.g.,
change the optimizer method, increase iterations allotted)
4. Simplify the model by removing correlation between random effects
(see Brown
(2021) p. 10)
5. Last resort, consider a simpler random effects structure - but
document this decision fully, the random effects structure should be
theoretically motivated so this is a compromise.
Step 3.2 - use afex::all_fit to test different optimizers
Note: if you were unable to install the “optimx” package then
skip to Step 3.3
We will try changing the optimizer (strategy 3). First, use the
lme4::allFit() function (include only the model we fit
above as an argument). This function will try several different
optimizer algorithms and report warning messages for each. Then you can
pick one optimizer that converges and use it. *Note:
lme4::allFit() replaces the afex::all_fit()
function that is referred to in Brown (2021).
Click the button for the code to run the all_fit() function
and see the output:
allFit(rt_full_mod)## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## optimx.L-BFGS-B :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0240828 (tol = 0.002, component 1)## [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00220168 (tol = 0.002, component 1)## [OK]## original model:
## RT ~ 1 + modality + (1 + modality | PID) + (1 + modality | stim)
## data: rt_data
## optimizers (6): bobyqa, Nelder_Mead, nlminbwrap, optimx.L-BFGS-B,nloptwrap.NLOPT_LN_NELDERME...
## differences in negative log-likelihoods:
## max= 0.00015 ; std dev= 6.1e-05
Step 3.3 - Re-fit the model with a different optimizer
The output of lme4::allFit() shows us that the “bobyqa”,
“Nelder_Mead”, and a 2 other methods all converge (others do not). When
we first ran the model above, by default lmer() used the
“nloptwrap” method.
So we will use one of the converging methods in its place now - let’s
use “bobyqa”. We do this by passing an extra argument to
lmer() called control with the value set to be
control = lmerControl(optimizer = "bobyqa"). Then callsummary()`
like before
rt_full_mod <- lmer(RT ~ 1 + modality + (1 + modality|PID) + (1 + modality|stim),
data = rt_data,
control = lmerControl(optimizer = "bobyqa"))
summary(rt_full_mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: RT ~ 1 + modality + (1 + modality | PID) + (1 + modality | stim)
## Data: rt_data
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 302385.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3646 -0.6964 -0.0140 0.5886 5.0003
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## stim (Intercept) 303.9 17.43
## modalityAudiovisual 216.7 14.72 0.16
## PID (Intercept) 28552.8 168.98
## modalityAudiovisual 7709.8 87.81 -0.17
## Residual 65258.8 255.46
## Number of obs: 21679, groups: stim, 543; PID, 53
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1044.14 23.36 44.704
## modalityAudiovisual 83.18 12.58 6.615
##
## Correlation of Fixed Effects:
## (Intr)
## modltyAdvsl -0.178
Step 3.4 - hypothesis testing
The first thing you may notice in the output is that there is no
p-value for our fixed effect (if you do see a p-value then you probably
loaded the afex the lmerTest package, so you
ran a slightly different function). The recommended way to compute a
p-value to test a null hypothesis is by comparing a model with the
effect of interest to a model that is identical but with the fixed
effect of interest removed (Brown, 2021). Our
effect of interest is the fixed effect of modality, so we’ll remove it
in the reduced model (but leave the random effects as is, including
modality random slopes). Let’s do that now. First, fit a model called
rt_modalityremoved_mod, then use the anova()
function to compare it to the full model with a likelihood ratio test.
Use the same parameters to fit the reduced model that were used for the
full model.
rt_modalityremoved_mod <- lmer(RT ~ 1 + (1 + modality|PID) + (1 + modality|stim),
data = rt_data,
control = lmerControl(optimizer = "bobyqa"))
anova(rt_modalityremoved_mod,rt_full_mod)## refitting model(s) with ML (instead of REML)## Data: rt_data
## Models:
## rt_modalityremoved_mod: RT ~ 1 + (1 + modality | PID) + (1 + modality | stim)
## rt_full_mod: RT ~ 1 + modality + (1 + modality | PID) + (1 + modality | stim)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rt_modalityremoved_mod 8 302449 302513 -151217 302433
## rt_full_mod 9 302419 302491 -151200 302401 32.385 1 1.264e-08
##
## rt_modalityremoved_mod
## rt_full_mod ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(rt_full_mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: RT ~ 1 + modality + (1 + modality | PID) + (1 + modality | stim)
## Data: rt_data
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 302385.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3646 -0.6964 -0.0140 0.5886 5.0003
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## stim (Intercept) 303.9 17.43
## modalityAudiovisual 216.7 14.72 0.16
## PID (Intercept) 28552.8 168.98
## modalityAudiovisual 7709.8 87.81 -0.17
## Residual 65258.8 255.46
## Number of obs: 21679, groups: stim, 543; PID, 53
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1044.14 23.36 44.704
## modalityAudiovisual 83.18 12.58 6.615
##
## Correlation of Fixed Effects:
## (Intr)
## modltyAdvsl -0.178
Step 4 - Interpret the output of the model comparison and the full model summary
Model comparison (output of anova()):
- the “Chisq” value of 32.385 and low p-value (.00000001264) indicate that the full model (including the fixed effect of modality) fits better than the reduced model, so we can reject the null hypothesis that modality has no effect on response time.
Full model summary (output of summary() for the full
model):
- The coefficient/estimate for the effect of modality is 83.18,
meaning that a 1 unit increase in modality predicts an RT increase of
83.18 ms. In other words, participants were 83 ms slower in the
audiovisual (coded as 1) relative to the audio-only (coded as 0)
condition.
- The fixed intercept is 1044.14, meaning that when
modalityis 0 (i.e., in the Audio-only condition) the mean response time is 1044.14
- standard error estimates how much the coefficient varies across
samples
- the t-value is the coefficient divided by its standard error
We could report the results like this (from Brown
(2021)):
> “A linear mixed effects model was fit to trial response times, with
a fixed effect of modality, and random effects structure allowing
intercepts and slopes of the modality effect to vary by-participant and
by-item. A likelihood-ratio test indicated that the model including the
fixed effect of modality provided a better fit for the data than a model
without it, χ2(1) = 32.39, p < .001. Examination of the
summary output for the full model indicated that response times were on
average an estimated 83 ms slower in the audiovisual relative to the
audio-only condition (β = 83.18, SE = 12.58, t = 6.62).”
coef() function (you will get over 500 lines of
output)coef(rt_full_mod)## $stim
## (Intercept) modalityAudiovisual
## babe 1038.921 82.11513
## back 1050.914 86.52640
## bad 1041.122 81.12265
## bag 1042.896 86.40593
## bake 1039.395 81.75285
## balk 1042.560 84.17721
## ball 1035.114 80.50240
## bane 1055.740 83.15440
## bang 1048.652 85.51504
## bar 1042.990 82.57039
## base 1043.169 81.00300
## bash 1036.544 79.31502
## batch 1047.445 84.84042
## bath 1039.152 81.59483
## bathe 1043.921 77.52535
## bead 1037.029 78.75310
## beak 1042.654 84.36897
## beam 1040.067 79.28125
## bear 1040.514 83.76434
## beat 1049.312 84.73022
## beg 1044.947 85.99852
## bib 1048.505 83.33006
## bide 1046.706 84.50435
## big 1044.807 84.21705
## bin 1045.735 83.48758
## birth 1022.948 72.73989
## bit 1046.026 84.01296
## boat 1039.346 79.28316
## bog 1041.345 80.92928
## bomb 1045.924 85.51498
## booth 1039.411 80.25989
## booze 1051.868 86.61135
## botch 1040.909 78.78912
## bud 1034.365 79.13424
## budge 1044.947 84.51831
## bun 1034.580 79.54525
## burp 1033.454 75.26893
## bus 1033.021 79.51754
## bush 1035.837 74.66136
## but 1039.495 80.06999
## cab 1040.445 83.75955
## cage 1038.917 78.57026
## call 1049.019 85.14604
## calm 1047.082 87.14903
## cape 1040.686 78.54312
## cash 1042.070 86.83023
## catch 1036.478 79.55022
## caught 1046.641 86.60317
## cave 1054.375 88.71955
## cease 1045.759 86.50896
## cell 1050.967 92.27456
## chair 1029.443 73.05484
## chalk 1044.948 83.93337
## chase 1059.059 88.57720
## cheap 1052.476 87.50220
## cheat 1049.722 90.15106
## cheek 1035.924 79.15037
## chef 1056.613 91.35629
## chewed 1042.900 80.91129
## chin 1040.826 84.85155
## chirp 1033.836 76.85190
## choke 1054.485 90.25234
## choose 1037.296 79.01739
## chug 1042.780 80.73139
## chum 1047.075 86.13041
## cite 1034.088 80.05814
## coat 1053.603 91.02674
## cone 1038.833 81.21495
## cope 1057.129 90.66435
## core 1048.035 83.86187
## cot 1051.651 87.05233
## couch 1040.402 80.76850
## cough 1045.715 78.56381
## cub 1030.915 82.09348
## cuff 1046.868 86.57878
## cup 1037.053 80.95711
## curl 1033.687 75.00310
## curt 1034.330 77.60745
## curve 1051.558 89.37790
## czar 1046.183 82.42070
## dab 1046.619 84.25423
## dash 1037.596 79.10409
## date 1036.570 83.03279
## dawn 1037.905 76.16523
## days 1038.708 81.21403
## dear 1045.196 83.98793
## death 1040.987 80.77035
## deck 1049.092 87.14966
## dice 1035.631 80.64568
## did 1030.882 76.37728
## died 1042.946 80.63664
## dig 1039.994 81.75872
## dim 1040.517 85.43423
## dime 1036.267 75.99750
## dip 1041.740 80.13675
## dire 1048.370 82.92115
## ditch 1058.635 90.70649
## dive 1033.222 79.65575
## dock 1035.795 81.31098
## dodge 1048.735 84.84144
## does 1039.454 82.17599
## dog 1039.891 84.28055
## door 1033.471 77.61828
## dot 1031.332 76.88081
## doubt 1040.535 82.81439
## doug 1035.189 81.62758
## douse 1044.627 86.85940
## dove 1045.114 83.05949
## doze 1039.404 81.31138
## duke 1034.318 77.57747
## dull 1053.447 87.50619
## face 1034.552 76.66391
## fad 1053.157 89.38940
## fake 1047.266 86.25887
## fame 1042.216 76.66439
## fang 1056.606 87.39907
## fate 1042.709 84.33710
## fawn 1050.121 86.16232
## fed 1050.624 86.57675
## fern 1053.202 87.63107
## fetch 1026.848 74.19896
## fief 1038.397 80.43287
## fight 1038.246 82.56968
## fill 1048.568 83.73583
## fin 1055.797 86.71934
## fire 1037.484 78.56707
## fit 1044.405 83.27911
## five 1050.794 84.80899
## fizz 1043.524 81.92920
## foal 1046.642 83.33906
## foam 1038.740 78.93390
## fog 1032.641 76.77786
## foil 1045.297 82.55937
## folk 1041.681 80.78197
## fought 1048.857 84.01334
## foul 1044.304 82.80384
## fudge 1037.576 81.84432
## fuzz 1046.588 85.08888
## gab 1055.214 87.85660
## gag 1040.243 83.27562
## gain 1050.870 85.24620
## gal 1048.022 88.26238
## gall 1045.231 82.68150
## game 1039.056 78.09731
## gang 1048.991 83.11020
## gave 1047.072 84.53772
## gawk 1040.297 82.07796
## gear 1048.268 84.41224
## gel 1040.727 84.26028
## gem 1036.839 77.17615
## get 1038.837 80.07484
## ghoul 1047.729 85.36133
## gig 1049.111 85.43305
## give 1043.300 84.10929
## gnat 1039.071 86.00424
## gnome 1051.898 86.95742
## goal 1047.900 88.15657
## gob 1041.922 80.02941
## gone 1049.454 86.30119
## gong 1042.177 82.48138
## goof 1043.616 82.18387
## goon 1047.555 80.33854
## goose 1043.661 89.58170
## gore 1043.550 87.56203
## got 1028.070 72.92420
## gown 1032.794 77.64271
## guide 1053.576 87.74422
## guile 1049.651 86.22694
## gum 1042.929 81.95106
## hack 1033.212 72.79558
## had 1056.474 89.79908
## hair 1035.562 75.88222
## half 1055.069 90.36750
## hall 1045.000 78.93671
## hash 1045.853 89.40260
## hat 1038.921 79.46537
## hate 1029.835 76.25324
## have 1052.226 89.29447
## head 1053.048 87.21683
## hear 1040.925 79.32771
## heat 1044.058 83.03276
## heed 1044.397 81.81069
## height 1039.468 80.97627
## hem 1047.760 83.68884
## hid 1030.819 77.66540
## hike 1046.925 88.09888
## hill 1045.766 86.79680
## hip 1027.933 75.13207
## hire 1043.359 84.57237
## hiss 1042.363 84.64477
## hit 1040.295 76.04834
## hitch 1044.371 81.72511
## hive 1043.310 83.95391
## hog 1050.314 84.20909
## home 1049.890 87.47964
## hone 1049.453 84.30109
## hoof 1038.460 80.14797
## hook 1029.895 76.23789
## hose 1042.307 82.44217
## house 1041.671 85.98863
## huff 1052.831 84.78264
## hug 1031.841 80.62322
## hull 1044.747 85.12870
## hung 1050.495 89.46245
## hurt 1038.114 80.32454
## hush 1032.224 80.90897
## hut 1057.149 91.93565
## hutch 1041.535 84.58059
## jab 1048.368 82.74098
## jail 1044.345 80.30181
## jam 1045.090 87.85350
## jar 1046.065 82.50548
## jazz 1050.055 87.11876
## jeer 1043.386 90.41935
## jewel 1036.120 80.82945
## job 1051.172 90.10607
## jock 1039.001 74.91884
## jog 1046.207 86.46544
## jowl 1027.885 77.33049
## judge 1036.912 82.00576
## jug 1044.486 84.41406
## juice 1045.157 82.95526
## juke 1046.280 84.24499
## kale 1032.415 77.02185
## kid 1034.978 81.12625
## kill 1046.351 84.34470
## kin 1036.976 80.18406
## kite 1040.837 81.98890
## kneel 1040.971 78.96488
## knit 1058.487 90.19146
## known 1049.701 84.77769
## lab 1041.940 81.94606
## lace 1037.523 81.57194
## lad 1046.570 82.58021
## lamb 1050.551 84.52504
## latch 1035.295 77.17097
## late 1039.832 82.65885
## lawn 1041.522 83.11971
## leaf 1043.802 85.12625
## league 1059.531 92.30763
## leak 1038.067 77.24232
## leave 1055.274 87.47576
## ledge 1038.830 81.44114
## leg 1049.767 84.41886
## less 1046.770 80.75696
## let 1037.587 80.30923
## lied 1060.315 88.13072
## life 1047.393 87.24208
## like 1053.049 89.35596
## lime 1054.262 87.92199
## line 1042.631 80.72760
## load 1034.127 81.92004
## lob 1042.670 82.09141
## lobe 1048.628 84.05326
## lodge 1035.484 77.99526
## log 1046.951 81.64542
## lone 1040.065 82.17616
## loose 1052.169 86.28392
## loot 1042.175 81.71561
## lop 1049.378 88.01630
## lore 1038.699 78.73430
## loud 1042.115 85.30276
## louse 1055.125 90.88445
## lug 1055.800 88.50842
## luge 1057.174 91.53722
## mail 1044.915 85.50498
## maim 1042.724 82.31876
## mall 1048.419 88.54724
## man 1037.297 83.64758
## mass 1036.112 78.22457
## mat 1055.562 89.98442
## math 1049.406 85.80852
## mauve 1054.478 90.32598
## maze 1047.851 86.82600
## mead 1051.824 83.94758
## meal 1032.525 75.84278
## meek 1048.708 93.00292
## men 1031.343 81.29828
## mess 1044.152 82.58488
## might 1044.477 83.70255
## mill 1043.427 88.04708
## miss 1051.055 85.65212
## mitt 1046.520 82.73159
## mob 1046.300 81.09604
## mole 1048.742 87.84247
## moon 1043.335 84.76673
## moot 1046.066 82.33530
## mop 1058.143 89.48166
## mope 1047.581 84.08314
## moth 1043.878 83.01256
## mouse 1040.950 81.26744
## mouth 1043.925 81.07537
## mud 1043.979 82.25271
## muff 1033.502 76.85967
## mum 1051.619 84.76574
## myth 1051.944 87.61113
## nab 1052.703 83.27204
## name 1042.368 85.79213
## nape 1046.271 82.35672
## nash 1044.102 84.36209
## neat 1039.380 79.00130
## neck 1041.276 85.96619
## need 1053.045 84.24198
## nerve 1045.939 83.56613
## newt 1048.357 82.53336
## niece 1037.787 75.62437
## night 1052.631 82.61608
## nod 1063.382 88.84001
## none 1039.936 79.73739
## noon 1049.854 86.69282
## nose 1041.298 79.44303
## not 1048.859 84.43756
## notch 1035.649 81.58975
## null 1036.649 76.29606
## pair 1060.324 93.20114
## pan 1044.253 84.08786
## pass 1047.290 83.56259
## peace 1021.884 76.96530
## peach 1032.528 74.07108
## peak 1043.150 82.73248
## peat 1032.655 76.14684
## peep 1035.021 80.66507
## peeve 1049.469 86.97180
## pen 1045.176 85.58502
## perch 1046.296 87.01463
## perm 1042.549 84.15966
## pet 1035.778 79.20365
## phase 1037.388 81.19836
## phone 1040.226 76.50912
## pick 1048.105 87.14125
## pile 1039.829 84.99850
## pin 1035.327 77.83476
## ping 1038.957 81.51854
## pipe 1038.738 79.05160
## pitch 1049.535 82.15329
## pod 1049.224 82.11442
## poise 1045.467 85.12786
## pooch 1048.726 85.13146
## pool 1031.427 81.95062
## pop 1039.851 83.78195
## pose 1040.228 78.89918
## pouch 1048.508 86.99761
## pour 1020.649 75.59048
## pout 1044.446 85.35956
## puff 1046.550 84.63828
## pun 1046.204 83.41679
## pup 1042.187 82.84781
## putt 1051.320 86.30392
## rack 1052.194 88.42031
## rag 1049.513 85.34213
## raid 1045.075 85.21239
## rail 1043.475 80.06065
## rake 1048.632 85.15197
## ram 1046.203 85.17706
## ran 1033.117 77.91632
## rang 1037.375 82.99756
## rap 1036.655 74.45906
## rash 1062.557 90.51986
## rate 1025.637 75.05706
## rear 1043.327 81.29075
## red 1050.696 82.86388
## ref 1047.366 83.69489
## rev 1041.169 81.84952
## rib 1048.010 89.44283
## rich 1037.605 76.67944
## ridge 1060.859 89.01771
## right 1050.706 91.67910
## ring 1053.756 90.31850
## rise 1044.146 85.83980
## roach 1045.981 82.95105
## rogue 1038.807 80.01357
## roof 1047.225 83.69871
## rook 1041.931 84.79130
## room 1056.964 96.29308
## root 1037.153 75.93280
## rose 1048.279 87.77433
## rot 1038.003 78.05591
## rouge 1047.958 81.11532
## rove 1042.736 82.53232
## rude 1048.663 85.51996
## rule 1045.350 87.68941
## rut 1030.108 76.84224
## sack 1041.282 78.54588
## safe 1047.540 86.21768
## sage 1039.689 82.04662
## sail 1044.078 83.95963
## same 1045.230 85.12612
## sane 1043.497 79.31376
## sang 1043.811 80.94921
## sap 1047.028 84.43262
## sash 1054.596 91.60744
## sauce 1037.759 83.45014
## save 1054.777 88.22368
## seem 1033.716 77.48933
## seen 1035.783 75.03905
## sees 1049.030 85.99762
## serve 1047.696 83.73240
## set 1048.900 87.96606
## sewed 1053.234 86.18341
## shack 1042.167 84.73908
## shade 1046.022 86.01450
## sham 1055.725 89.49458
## shave 1054.705 87.20292
## shear 1031.092 77.33960
## sheath 1052.023 82.78834
## shed 1040.827 84.25482
## sheep 1051.119 85.38800
## sheet 1038.318 79.86355
## shell 1046.058 83.42337
## shim 1053.477 91.12424
## shin 1047.999 83.94257
## shine 1058.400 88.40587
## ship 1058.944 86.36568
## shook 1050.533 86.07398
## shoot 1057.926 87.96622
## shop 1040.877 80.48169
## shore 1041.212 80.23309
## shot 1053.560 85.04721
## shout 1045.631 82.74686
## shown 1047.444 86.41158
## shun 1036.757 82.00065
## shut 1036.340 80.70995
## sick 1043.321 80.00558
## side 1054.130 86.85368
## siege 1053.585 82.58664
## sieve 1049.320 85.29026
## sing 1043.518 77.59442
## sip 1037.098 76.41068
## sit 1055.054 91.27707
## soar 1045.964 83.13872
## sob 1049.354 89.09940
## sock 1041.142 79.43845
## sod 1059.770 90.60555
## soil 1048.732 80.73996
## some 1051.421 89.68941
## song 1053.748 87.10236
## soothe 1048.824 87.14040
## south 1043.370 84.77203
## sown 1055.497 87.34025
## sub 1048.541 83.47911
## such 1041.860 81.79824
## tag 1039.211 80.82300
## tail 1041.578 81.15355
## take 1047.383 83.86050
## talk 1043.455 80.17272
## tall 1045.471 86.83288
## tap 1043.258 80.37828
## tar 1035.998 78.06458
## taught 1043.637 82.89166
## teach 1032.380 76.94175
## tease 1033.791 75.30195
## tech 1039.392 81.65804
## teethe 1069.811 93.19231
## tell 1035.733 78.58382
## ten 1043.146 86.01562
## than 1034.099 81.12536
## thatch 1049.705 84.13242
## their 1039.480 82.88014
## them 1049.661 83.35430
## theme 1038.007 78.09569
## these 1048.204 84.25600
## thin 1046.338 83.67447
## thine 1053.484 91.28764
## those 1037.470 81.00035
## thug 1056.629 88.59648
## thus 1048.983 86.75589
## tide 1055.038 88.95381
## ties 1050.464 84.31076
## tight 1058.002 84.35955
## time 1042.171 79.42650
## tin 1036.582 82.24119
## tip 1046.729 87.79698
## toad 1047.643 84.62195
## toil 1050.046 89.78463
## toll 1040.978 81.53409
## tone 1027.634 76.98424
## top 1050.038 90.07484
## tote 1031.957 74.84663
## touch 1041.758 79.51840
## tough 1048.554 86.93229
## tout 1044.901 86.62052
## town 1048.702 81.11463
## tub 1038.111 77.01853
## tube 1044.586 76.69827
## use 1041.596 83.12440
## vague 1047.785 84.44383
## van 1053.601 85.54461
## vane 1051.523 90.12707
## vase 1042.641 82.29998
## veal 1039.174 79.67358
## verge 1043.832 80.22405
## verse 1046.594 84.56838
## vet 1039.852 80.92725
## vice 1038.963 82.34232
## vile 1044.623 82.69476
## vine 1038.305 80.89245
## vogue 1044.954 77.68277
## voice 1047.961 84.33096
## wade 1054.396 88.92327
## wage 1049.496 92.53737
## walk 1043.360 79.97585
## war 1043.699 83.18414
## watch 1030.853 77.33569
## wave 1043.158 81.72411
## web 1046.780 84.93750
## wed 1050.873 85.64588
## wedge 1049.558 86.09362
## weed 1048.798 84.90534
## week 1042.000 74.80555
## weep 1042.335 81.80561
## well 1041.712 81.23102
## wet 1041.312 79.81610
## whale 1035.809 78.57471
## whole 1051.191 88.42446
## wick 1052.483 85.12277
## wife 1043.994 81.43968
## will 1049.199 84.49703
## wine 1045.220 86.53611
## wipe 1038.198 79.55266
## wire 1034.326 78.78424
## woke 1045.808 85.73513
## womb 1047.488 80.21622
## won 1043.131 86.43064
## wool 1035.060 78.73339
## work 1038.318 75.58146
## wort 1052.709 84.56164
## worth 1040.409 82.85559
## wove 1052.240 87.98749
## wrath 1041.035 80.92772
## wren 1044.410 81.95183
## writhe 1042.264 81.78374
## wrong 1055.399 85.59014
## wrote 1042.828 80.63865
## yak 1044.860 81.99311
## yam 1048.245 90.16337
## yawn 1039.379 81.14480
## year 1037.308 82.97026
## yen 1044.379 83.33265
## yolk 1032.187 81.66625
## young 1032.830 78.03479
## your 1052.891 90.65923
## zeal 1042.742 77.45131
## zing 1045.024 83.55766
##
## $PID
## (Intercept) modalityAudiovisual
## 301 1024.0669 -16.936434
## 302 1044.1377 1.842595
## 303 882.8306 57.789295
## 304 1232.7545 -27.919818
## 306 1042.3420 33.886522
## 307 1111.3632 -9.939667
## 308 1250.7673 71.164806
## 309 795.2446 15.481884
## 310 1176.1359 104.151774
## 311 1012.9321 29.956256
## 312 1109.8745 153.103921
## 313 1114.2740 29.601075
## 314 1169.1517 -73.242657
## 315 877.9771 17.015014
## 316 1419.5222 -37.700573
## 318 945.4723 58.220367
## 319 1017.4086 98.276016
## 320 987.4941 101.341112
## 321 1025.1270 139.711305
## 324 1031.3305 136.247296
## 325 826.6336 34.599816
## 326 1048.6754 40.359209
## 328 1236.8998 128.530597
## 329 1042.3467 10.735128
## 331 1406.8491 160.773353
## 333 1644.5199 56.375054
## 334 943.7441 93.903857
## 335 1171.5855 64.055478
## 337 872.3788 169.983219
## 338 806.7056 121.221931
## 340 1179.8488 104.695827
## 341 867.6691 123.864759
## 342 1253.5632 30.102294
## 343 987.7640 208.128262
## 344 1027.5597 96.902626
## 346 895.7248 36.422803
## 348 755.2919 188.633628
## 349 940.6337 54.966477
## 350 1073.1133 302.721654
## 351 1120.9818 214.077525
## 352 796.4338 78.590348
## 353 1103.8172 155.151511
## 354 1074.6543 154.892053
## 355 1192.9344 89.264166
## 356 907.2052 397.764490
## 357 910.5672 81.928221
## 358 963.0962 33.338236
## 359 1087.6474 -39.273827
## 360 1070.4021 26.771297
## 361 982.1810 131.300517
## 362 953.3700 56.194587
## 363 920.7194 43.641049
## 364 1003.5099 75.844926
##
## attr(,"class")
## [1] "coef.mer"
If you are interested in seeing more examples, you can find Violet
Brown’s R markdown document for this data at https://osf.io/v6qag/.
- For a great linear mixed modeling activity (with more visualization)
check out this
activity in PsyTeachR
- For an example publication utilizing generalized mixed modeling
(dichotomous outcome variable), see: Raio, C.M., Konova, A.B. &
Otto, A.R. Trait impulsivity and acute stress interact to influence
choice and decision speed during multi-stage decision-making. Sci
Rep 10, 7754 (2020). https://doi.org/10.1038/
That’s all for this activity!
Extras
Table from Speekenbrink (2023):
.
sjPlot::plot_model(rt_full_mod,type = "eff", terms = "modality") +
coord_cartesian(ylim = c(500,1200)) +
theme_classic()