Mixed Models Part 2
Overview
In this second workshop on mixed models we will look at how to build LMMs for factorial design, including ensuring we choose a contrast coding scheme for our experimental factors that allows us to interpret our parameter estimates appropriately. We’ll also look at generalised linear mixed models (GLMMs) in the context of modelling binomial data (where our DV is a 0 or a 1), and ordinal mixed models for cases where our DV is measured on an ordinal scale (as might be the case with Likert-scale data).
Slides
You can download the slides in .odp format by clicking here and in .pdf format by clicking on the image below.
Once you’ve watched the video above, run the code below on your own machines.
Mixed Models
2 x 2 Factorial Design
In this first case imagine that we have a 2 x 2 repeated measures design. The first factor is Context (Negative vs. Positive) and the second is Sentence Type (Negative vs. Positive). The DV is reading time duration to a Target Sentence (measured in ms.). We have 60 subjects, and 28 items.
Let’s first load the libraries that we’re going to use.
library(tidyverse)
library(lme4)
library(lmerTest)
library(emmeans)
library(performance)
library(pbkrtest)Now we’re going to read in our data and turn our
subject, item, context, and
sentence columns into factors.
factorial_data <- read_csv("https://raw.githubusercontent.com/ajstewartlang/15_mixed_models_pt1/master/data/2x2.csv")
tidied_factorial_data <- factorial_data %>%
mutate(subject = factor(Subject), item = factor(Item), RT = RT,
context = factor(Context), sentence = factor(Sentence))When we generate summary statistics, we encounter the issues of NAs where we expected to see mean values.
tidied_factorial_data %>%
group_by(context, sentence) %>%
summarise(mean_rt = mean(RT), sd_rt = sd(RT))## `summarise()` has grouped output by 'context'. You can override using the
## `.groups` argument.## # A tibble: 4 × 4
## # Groups: context [2]
## context sentence mean_rt sd_rt
## <fct> <fct> <dbl> <dbl>
## 1 Negative Negative 1474. 729.
## 2 Negative Positive NA NA
## 3 Positive Negative NA NA
## 4 Positive Positive 1579. 841.This indicates we likely have missing data. R never decides what it should do about missing data - you need to make that decision. Let’s use the {visdat} package to investigate…
library(visdat)vis_miss(tidied_factorial_data)## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
We can calculate how much data RT is missing another way…
tidied_factorial_data %>%
filter(is.na(RT)) %>%
count()## # A tibble: 1 × 1
## n
## <int>
## 1 12So we can see we have 12 missing data points. We can filter our cases where we have missing data and re-do our summary statistics.
tidied_factorial_data %>%
filter(!is.na(RT)) %>%
group_by(context, sentence) %>%
summarise(mean_rt = mean(RT), sd_rt = sd(RT))## `summarise()` has grouped output by 'context'. You can override using the
## `.groups` argument.## # A tibble: 4 × 4
## # Groups: context [2]
## context sentence mean_rt sd_rt
## <fct> <fct> <dbl> <dbl>
## 1 Negative Negative 1474. 729.
## 2 Negative Positive 1595. 887.
## 3 Positive Negative 1633. 877.
## 4 Positive Positive 1579. 841.Let’s visualise the data next. I’m using the
stat_summary() function to add means and bootstrapped
condidence intervals around those means for each of our conditions. I
use the syntax context:sentence in the aes()
expression to indicate that I want to plot all combinations of these two
factors (i.e., the interactions). I use a few other aesthetic tweaks to
make the graph look a little better than the default. I should have
added a title - can you modify the code below to do that?
tidied_factorial_data %>%
filter(!is.na(RT)) %>%
ggplot(aes(x = context:sentence, y = RT, colour = context:sentence)) +
geom_violin() +
geom_jitter(width = .1, alpha = .2) +
stat_summary(fun.data = "mean_cl_boot", colour = "black") +
guides(colour = FALSE) +
labs(x = "Context X Sentence",
y = "RT (ms.)") +
theme_minimal() +
coord_flip()## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.## Warning: Computation failed in `stat_summary()`:
Before we build our model, we need to set the contrast coding for our factors. By default, R using dummy (treatment) coding (think back to how we used this in the context of understanding ANOVA as a case of regression). The problem with dummy coding for factorial designs is that you can end up misinterpreting simple effects (e.g., an effect of Factor 2 at one level of Factor 1) as a main effect (e.g., an effect of Factor 2 at the average of Factor 1). To address this, we can use sum or deviation coding. This will result in the Intercept of our model correspoding to the grand mean of our conditions (i.e., the mean of means) and makes the interpretation of our fixed effects (and any interaction effects) more straightforward.
If you’re interested in reading more about this topic, I recommend this great paper by Schad and colleagues. Just click on the image below to be access it.
Let’s set the contrast coding of our factors using deviation coding as follows. This will allow us to compare each of our conditions to the average of the other conditions.
contrasts(tidied_factorial_data$context) <- matrix(c(.5, -.5))
contrasts(tidied_factorial_data$sentence) <- matrix(c(.5, -.5))Now that our contrasts are coded we can go ahead and build our mixed
model. Note that the maximal model did not converge so we dropped the
interaction term (context:sentence) from our subject random
effect. For this random effect, we are modelling just additive effects
of context and sentence.
factorial_model <- lmer(RT ~ context * sentence +
(1 + context + sentence | subject) +
(1 + context * sentence | item),
data = tidied_factorial_data)Let’s check the assumptions of our model.
check_model(factorial_model)
So it looks like we may have an issue with the normality of our residuals. Things looks ok apart from the right 20% of our residuals. We may want to try to model under a different distribution. We can plot our RT values on a Cullen and Frey plot.
library(fitdistrplus)missing_data_removed <- tidied_factorial_data %>%
filter(!is.na(RT))
descdist(missing_data_removed$RT)
## summary statistics
## ------
## min: 204 max: 7000
## median: 1387.5
## mean: 1569.854
## estimated sd: 836.5866
## estimated skewness: 1.96018
## estimated kurtosis: 9.283192On the Cullen and Frey plot we see our data is quite close to a Gamma distribution. We can try to model our data using a generalised linear model assuming sampling from the Gamma distribution as follows. One of the challenges with such models is that very often the random effects structure needs to be radically simplified. This can increase the Type I error rate and may result in us thinking we have an effect when really we don’t. This is one of the challenges of building models in general. To paraphrase George Box, all models are wrong but some are useful. Let’s try to build a Gamma model anyway.
gamma_factorial_model <- glmer(RT ~ context * sentence +
(1 + context + sentence | subject) +
(1 + sentence | item),
family = Gamma,
nAGQ = 0,
data = tidied_factorial_data)In order to fit this model, I had to simplify the random effects terms and set nAGQ to 0 (its default is 1). This means our parameter estimates are a little less exact than if we had gone with the default (but at least the model converges on a solution).
Let’s look at the summaries of our two models and see if they differ.
First we’ll look at our factorial_model which we built
using lmer().
summary(factorial_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ context * sentence + (1 + context + sentence | subject) +
## (1 + context * sentence | item)
## Data: tidied_factorial_data
##
## REML criterion at convergence: 26628.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3460 -0.5439 -0.1563 0.3202 7.1297
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 106850 326.88
## context1 21848 147.81 -0.68
## sentence1 28914 170.04 -0.08 -0.28
## item (Intercept) 105902 325.43
## context1 4919 70.14 0.27
## sentence1 163485 404.33 -0.02 -0.11
## context1:sentence1 56612 237.93 -0.71 -0.81 -0.21
## Residual 432879 657.94
## Number of obs: 1668, groups: subject, 60; item, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1568.68 76.31 50.04 20.556 <2e-16 ***
## context1 -68.87 39.77 31.62 -1.732 0.0931 .
## sentence1 -36.35 85.81 29.71 -0.424 0.6749
## context1:sentence1 -168.45 78.68 26.32 -2.141 0.0417 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cntxt1 sntnc1
## context1 -0.109
## sentence1 -0.025 -0.070
## cntxt1:snt1 -0.327 -0.147 -0.112So we have an interaction between Context and
Sentence.
What about our Gamma model built using glmer()?
summary(gamma_factorial_model)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: Gamma ( inverse )
## Formula: RT ~ context * sentence + (1 + context + sentence | subject) +
## (1 + sentence | item)
## Data: tidied_factorial_data
##
## AIC BIC logLik deviance df.resid
## 26590.3 26666.1 -13281.1 26562.3 1654
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1002 -0.6623 -0.2033 0.4071 7.4220
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 8.326e-09 9.125e-05
## context1 3.188e-10 1.785e-05 -1.00
## sentence1 8.385e-09 9.157e-05 0.08 -0.08
## item (Intercept) 1.137e-01 3.372e-01
## sentence1 1.071e-01 3.273e-01 0.01
## Residual 1.469e-01 3.833e-01
## Number of obs: 1668, groups: subject, 60; item, 28
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 6.995e-04 6.372e-02 0.011 0.99124
## context1 1.964e-05 1.176e-05 1.670 0.09484 .
## sentence1 1.035e-05 6.186e-02 0.000 0.99987
## context1:sentence1 6.323e-05 2.312e-05 2.735 0.00624 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cntxt1 sntnc1
## context1 0.000
## sentence1 0.010 0.000
## cntxt1:snt1 0.000 0.032 0.000OK, this is interesting. The interaction between Context
and Sentence is still there. You are probably fine to
report either the results of the LMM or the GLMM. Whichever you choose,
be transparent and open when you report them. Indeed, you might even
want to report both models but highlight the fact that the interaction
term is significant regardless of which model is chosen.
To interpret the interaction, you’ll need to run pairwise comparisons
using emmeans().
emmeans(factorial_model, pairwise ~ context*sentence, adjust = "none")## $emmeans
## context sentence emmean SE df lower.CL upper.CL
## Negative Negative 1474 80.9 38.9 1310 1638
## Positive Negative 1627 99.0 43.9 1428 1827
## Negative Positive 1595 96.5 37.1 1399 1790
## Positive Positive 1579 90.2 48.4 1398 1761
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Negative Negative - Positive Negative -153.1 51.7 28.4 -2.960 0.0061
## Negative Negative - Negative Positive -120.6 90.3 29.1 -1.335 0.1922
## Negative Negative - Positive Positive -105.2 92.0 29.0 -1.144 0.2621
## Positive Negative - Negative Positive 32.5 97.1 31.4 0.335 0.7399
## Positive Negative - Positive Positive 47.9 98.3 28.4 0.487 0.6301
## Negative Positive - Positive Positive 15.4 59.9 28.7 0.256 0.7996
##
## Degrees-of-freedom method: kenward-rogerWe can see the interaction is being driven by reading times to Negative sentences in Negative vs. Positive contexts.
Is this also the case for our Gamma model?
emmeans(gamma_factorial_model, pairwise ~ context*sentence, adjust = "none")## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $emmeans
## context sentence emmean SE df asymp.LCL asymp.UCL
## Negative Negative 0.000730 0.0711 Inf -0.139 0.140
## Positive Negative 0.000679 0.0711 Inf -0.139 0.140
## Negative Positive 0.000688 0.0705 Inf -0.138 0.139
## Positive Positive 0.000700 0.0705 Inf -0.138 0.139
##
## Results are given on the inverse (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Negative Negative - Positive Negative 5.13e-05 1.68e-05 Inf 3.059 0.0022
## Negative Negative - Negative Positive 4.20e-05 6.19e-02 Inf 0.001 0.9995
## Negative Negative - Positive Positive 3.00e-05 6.19e-02 Inf 0.000 0.9996
## Positive Negative - Negative Positive -9.29e-06 6.19e-02 Inf 0.000 0.9999
## Positive Negative - Positive Positive -2.13e-05 6.19e-02 Inf 0.000 0.9997
## Negative Positive - Positive Positive -1.20e-05 1.62e-05 Inf -0.738 0.4604
##
## Note: contrasts are still on the inverse scaleWe see here that the pattern of pairwise comparisons is the same as with our previous linear model - the interaction is still being driven by reading times to Negative sentences in Negative vs. Positive contexts.
Binomial Data
In the following video I will examine how you build a generalised linear mixed model where your dependent variable isn’t continuous. I’ll focus on binomial data (where your DV is either a 1 or a 0).
Once you’ve watched the video above, run the code below on your own machines.
One Factor With Three Levels
In an eye-movement study we measured readers’ regressions as they read sentences. We had sentences conveying three different types of meaning (Negative ve. Neutral vs. Positive). For each, we measures where people did or did not make a regressive eye-movement. This is our DV and is coded as 0 or 1. We had 24 subjects and 24 items.
First, let’s read in the data and ensure the appropriate columns are coded as factors.
regressions_data <- read_csv("https://raw.githubusercontent.com/ajstewartlang/16_mixed_models_pt2/master/data/regressions.csv")
tidied_regressions_data <- regressions_data %>%
mutate(Subject = factor(Subject), Item = factor(Item),
Condition = factor(Condition), DV = DV)
str(tidied_regressions_data)## tibble [553 × 5] (S3: tbl_df/tbl/data.frame)
## $ Subject : Factor w/ 24 levels "S1","S10","S11",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Item : Factor w/ 24 levels "I1","I10","I11",..: 12 18 19 20 21 22 23 24 2 3 ...
## $ Condition: Factor w/ 3 levels "Negative","Neutral",..: 3 1 2 3 1 2 3 1 2 3 ...
## $ Region : chr [1:553] "R3" "R3" "R3" "R3" ...
## $ DV : num [1:553] 0 0 0 0 0 1 0 0 0 0 ...Let’s work out some descriptives - for each of the three conditions we can calculate the average number of eye-movement regressions (and the associated standard deviation).
tidied_regressions_data %>%
group_by(Condition) %>%
summarise(mean_DV = mean(DV), sd_DV = sd(DV))## # A tibble: 3 × 3
## Condition mean_DV sd_DV
## <fct> <dbl> <dbl>
## 1 Negative 0.247 0.433
## 2 Neutral 0.265 0.443
## 3 Positive 0.269 0.445So, in terms of the average number of eye-movement regressions, things look pretty similar from condition to condition. Let’s build a binomial model to check.
The maximal model (i.e., with random intercepts and slopes for our
subjects and for our items) doesn’t converge so we need to simplify the
random effects structure. The following is the most complex one we can
find to fit our data. Note, we are using the glmer()
function and specifying the distribution family as
binomial.
binomial_model <- glmer(DV ~ Condition + (1 | Subject),
data = tidied_regressions_data,
family = binomial)summary(binomial_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: DV ~ Condition + (1 | Subject)
## Data: tidied_regressions_data
##
## AIC BIC logLik deviance df.resid
## 603.7 621.0 -297.9 595.7 549
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3096 -0.6118 -0.4044 0.7636 3.1872
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 0.8363 0.9145
## Number of obs: 553, groups: Subject, 24
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3096 0.2647 -4.947 7.55e-07 ***
## ConditionNeutral 0.1093 0.2539 0.431 0.667
## ConditionPositive 0.1162 0.2515 0.462 0.644
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CndtnN
## ConditnNtrl -0.485
## ConditnPstv -0.488 0.507So, it doesn’t look like there’s much going on in our data. We can also compare the binomial model with the fixed effect (above) to a model with only the random effect.
null_binomial_model <- glmer(DV ~ (1 | Subject),
data = tidied_regressions_data,
family = binomial)We can then use the Likilhood Ratio Test to see if they differ (they don’t).
anova(binomial_model, null_binomial_model)## Data: tidied_regressions_data
## Models:
## null_binomial_model: DV ~ (1 | Subject)
## binomial_model: DV ~ Condition + (1 | Subject)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null_binomial_model 2 599.97 608.60 -297.98 595.97
## binomial_model 4 603.70 620.97 -297.85 595.70 0.2617 2 0.8773If our binomial model (or the comparison between our binomial and
null model) suggested an effect, we would use the emmeans()
package to run pairwise comparisons between the different levels of our
factor (as we did previously). We would also check our model assumptions
by creating a binned residuals plots as I have done in the video
above.
Ordinal Data
One Factor With Three Levels
In this final video in this session we will examine how to model ordinal data. Such data are quite common, and people often (incorrectly) model them using ANOVA. However, as the data are ordinal (where differences between adjacent points may not be equivalent), rather than continuous, this is not appropriate.
Imagine we had 42 participants rate images of sports on a scale of
0-10 corresponding to how much they liked each one. Before each rating
measure, they saw a video of a sport - called SportType - that matched
or mismatched the one they then had to rate (with a neutral video as
baseline). Our DV is Rating and our fixed effect is
Condition (Match vs. Mismatch vs. Neutral).
We’re going to use the {ordinal} package so let’s load
that into our library.
library(ordinal)Let’s read in our data an ensure our columns are coded appropriately. Our VideoCondition column is currently coded as 2, 3, or 4. 2 corresponds to the Match condition, 3 to Mismatch, and 4 as Neutral. Let’s recode the levels of that factor with the labels instead of the numbers.
Note, I am using the select() function to select the 4
key columns we need for our analysis. I have had to put
dplyr:: beforehand so that R know that it’s the
select() function that’s part of the {dplyr}
package from the Tidyverse that we want to use. There’s
also a function called select() in the {MASS}
package - oftentimes you will have both the Tidyverse
packages and {MASS} loaded so it’s important to remember
that you may need to make explicit which version of a particular
function you want to use.
ordinal_data <- read_csv("https://raw.githubusercontent.com/ajstewartlang/16_mixed_models_pt2/master/data/ordinal_data.csv")
ordinal_data_tidied <- ordinal_data %>%
mutate(Subject = factor(Subject), SportsType = factor(SportType)) %>%
mutate(Ratings = ratings) %>%
mutate(VideoCondition = as.character(VideoCondition)) %>%
mutate(VideoCondition = factor(recode(VideoCondition, "2" = "Match",
"3" = "Mismatch", "4" = "Neutral"))) %>%
dplyr::select(Subject, SportType, VideoCondition, Ratings)We need to do one more thing before our data set is ready for
modelling. We need to make sure that our DV is coded as an ordered
factor. We can do that simply with the following. Remember, the name
after the $ symbol refers to the column name in the data
frame that preceded the $.
ordinal_data_tidied$Ratings <- as.ordered(ordinal_data_tidied$Ratings)Let’s have a look at our data.
ordinal_data_tidied %>%
ggplot(aes(x = VideoCondition, y = Ratings, group = VideoCondition)) +
geom_jitter(aes(colour = VideoCondition), width = .1, alpha = .25, size = 1) +
theme_minimal() +
guides(colour = FALSE) +
theme(text = element_text(size = 14)) +
stat_summary(fun = "median", size = 2, alpha = .5)
It does look as if the Mismatch condition has a lower average rating
than the other two. Let’s build our ordinal mixed model to investigate.
Note, we use the clmm() function from the
{ordinal} package to do this. The syntax is pretty similar
to that used in the {lme4} models.
ordinal_model <- clmm(Ratings ~ VideoCondition +
(1 + VideoCondition | Subject) +
(1 + VideoCondition | SportType),
data = ordinal_data_tidied) ## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.Let’s also build a null model (i.e., with no fixed effects) to see how it compares to the one above.
null_ordinal_model <- clmm(Ratings ~ 1 +
(1 + VideoCondition | Subject) +
(1 + VideoCondition | SportType),
data = ordinal_data_tidied) ## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.anova(null_ordinal_model, ordinal_model)## Likelihood ratio tests of cumulative link models:
##
## formula:
## null_ordinal_model Ratings ~ 1 + (1 + VideoCondition | Subject) + (1 + VideoCondition | SportType)
## ordinal_model Ratings ~ VideoCondition + (1 + VideoCondition | Subject) + (1 + VideoCondition | SportType)
## link: threshold:
## null_ordinal_model logit flexible
## ordinal_model logit flexible
##
## no.par AIC logLik LR.stat df Pr(>Chisq)
## null_ordinal_model 22 10829 -5392.6
## ordinal_model 24 10826 -5389.0 7.3433 2 0.02543 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(ordinal_model)## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: Ratings ~ VideoCondition + (1 + VideoCondition | Subject) + (1 +
## VideoCondition | SportType)
## data: ordinal_data_tidied
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 2520 -5388.96 10825.93 6267(25072) 2.23e-03 1.8e+03
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 0.6673344 0.81691
## VideoConditionMismatch 0.0048430 0.06959 -1.000
## VideoConditionNeutral 0.0007284 0.02699 -1.000 1.000
## SportType (Intercept) 0.4860165 0.69715
## VideoConditionMismatch 0.0292862 0.17113 -0.407
## VideoConditionNeutral 0.0309516 0.17593 -0.429 1.000
## Number of groups: Subject 42, SportType 10
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## VideoConditionMismatch -0.3148 0.1019 -3.088 0.00202 **
## VideoConditionNeutral -0.2902 0.1025 -2.832 0.00463 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -4.8838 0.3090 -15.804
## 1|2 -3.7131 0.2791 -13.303
## 2|3 -2.4619 0.2676 -9.198
## 3|4 -1.5642 0.2643 -5.917
## 4|5 -0.8457 0.2630 -3.216
## 5|6 -0.2706 0.2625 -1.031
## 6|7 0.5937 0.2627 2.260
## 7|8 1.3795 0.2640 5.226
## 8|9 2.2840 0.2676 8.535
## 9|10 3.3942 0.2789 12.169Let’s examine which condition differs from each other condition using
{emmeans}.
emmeans(ordinal_model, pairwise ~ VideoCondition)## $emmeans
## VideoCondition emmean SE df asymp.LCL asymp.UCL
## Match 0.609 0.262 Inf 0.0948 1.123
## Mismatch 0.294 0.243 Inf -0.1832 0.771
## Neutral 0.319 0.245 Inf -0.1624 0.800
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Match - Mismatch 0.3148 0.1019 Inf 3.088 0.0057
## Match - Neutral 0.2902 0.1025 Inf 2.832 0.0129
## Mismatch - Neutral -0.0245 0.0855 Inf -0.287 0.9556
##
## P value adjustment: tukey method for comparing a family of 3 estimatesThe pairwise comparisons indicate that the Match vs Mismatch
conditions differ from each other, as do the Match vs. Neutral
conditions. You can change the adjustment for multiple comparisons using
the adjust = parameter in emmeans(). Does the
same story emerge when you use the Bonferroni adjustment?
emmeans(ordinal_model, pairwise ~ VideoCondition, adjust = "Bonferroni")Your Challenge
There are two parts to your challenge.
Part One
Read in the dataset “factorial_data.csv” from here:
https://raw.githubusercontent.com/ajstewartlang/16_mixed_models_pt2/master/data/factorial_data.csv.
These data are from a repeated measures experiment where participants
had to respond to a target word (measured by our DV which is labelled
Time). The target word always followed a prime word.
Prime and Target are our two factors – each
with two levels – Positive vs. Negative. We are interested in whether
there is a priming effect (i.e., Positive target words responded to more
quickly after Positive than after Negative Primes, and Negative target
words responded to more quickly after Negative than after Positive
Primes). We need to build the appropriate LMM to determine whether this
is indeed correct.
Build the appropriate linear mixed model to model the outcome
variable (Time) and determine whether it was predicted by
the fixed effects of Prime and Target.
Critically, we are interested in whether the interaction of these fixed
effects predicts our outcome variable.
Remember to set your contrast coding. If you find a singular fit
warning, consider simplyfing your random effects structure. When you
have a model you are happy with, use the
performance::check_model() function to check the model
assumptions. In your model, you should find a significant interaction
effect. Exploring this with emmeans() should reveal a
priming effect: Positive targets were responded to more quickly
following Positive vs. Negative primes, and Negative targets were
responded to more quickly following Negative vs. Positive primes. In
both cases, this priming effect was around 40ms.
Part Two
Read in the data file from here:
https://raw.githubusercontent.com/ajstewartlang/16_mixed_models_pt2/master/data/accuracy_data.csv
These data are from a study where 32 partcipants saw images of 32
faces and had to classify images as whether each was happy or sad. We
want to determine whether people were more accurate (indicatad by a 1 in
the Acc column) for Sad or Happy faces (FaceExpression column). Our
random effects are Subject and Face (which
corresponds to the trial ID).
Build the appropriate generalised linear mixed model to model the
binomial outcome variable (Acc) and determine whether it
was predicted by the fixed effect of FaceExpression. Model
both Subject and Face as random effects.
Your probably found that you needed to simplify your random effects
structure due to over-parameterisation - indicated by a
singular fit warning. Once you were able to build a model
that didn’t produce this error, you should have found that accuracy was
significantly better for Sad than for Happy
faces.
Improve this Workshop
If you spot any issues/errors in this workshop, you can raise an issue or create a pull request for this repo.