Mixed logistic regression

Binary outcomes with glmer()

Dr Massimiliano Canzi

Beyond the normal distribution

Did you know: the name axolotl comes from Nahuatl, the language of the Aztecs, and is thought to mean roughly “water monster”. Legend links it to Xolotl, a god who disguised himself as a salamander.

Source: National Geographic

Today

  • Binary outcomes, and why lm() is the wrong tool for them
  • The generalised linear model: distribution, linear predictor, link function
  • Probabilities, odds and log odds
  • Fitting a mixed logistic regression with glmer()
  • Interpreting coefficients: log odds, odds ratios and probabilities
  • Testing, convergence and reporting

So far

Every model in lectures 4 to 11 had a continuous outcome, such as reaction times, vowel durations or body mass, with roughly normal residuals.

Many linguistic outcomes are binary: one of two things happens.

  • A response is correct or incorrect
  • A speaker chooses one variant or the other: give her the book vs give the book to her
  • A listener says yes or no: is this sentence acceptable?
  • A final /t/ is pronounced or deleted

The data

A simulated lexical decision task: is the response correct (1) or not (0)?

set.seed(1)
n_part <- 30
n_word <- 60

participants <- tibble(participant = sprintf("P%02d", 1:n_part),
                       u0 = rnorm(n_part, 0, 0.8),     # own accuracy
                       u1 = rnorm(n_part, 0, 0.3))     # own frequency effect
words <- tibble(word     = sprintf("w%02d", 1:n_word),
                log_freq = round(runif(n_word, 1, 6), 2),
                w0       = rnorm(n_word, 0, 0.5))      # own difficulty

acc <- crossing(participants, words) |>
  mutate(freq_c  = log_freq - mean(words$log_freq),
         logit   = 1.5 + u0 + w0 + (0.8 + u1) * freq_c,
         correct = rbinom(n(), size = 1, prob = plogis(logit))) |>
  select(participant, word, log_freq, freq_c, correct)

rbinom(n, size = 1, prob) flips n weighted coins: 1 with probability prob, 0 otherwise. plogis() turns log odds into probabilities: more on both below.

The data

acc
# A tibble: 1,800 × 5
  participant word  log_freq freq_c correct
  <chr>       <chr>    <dbl>  <dbl>   <int>
1 P01         w01       5.96  2.38        1
2 P01         w02       3.48 -0.101       1
3 P01         w03       3.42 -0.161       1
4 P01         w04       1.87 -1.71        0
# ℹ 1,796 more rows
mean(acc$correct)
[1] 0.7461111

About 75% of responses are correct.

Accuracy and frequency

Grey: every response, 0 or 1. Orange: the proportion correct in 8 frequency bins. Frequent words are recognised more accurately.

Why not lm()?

A straight line predicts probabilities above 1 and below 0 (the pink zones). Its residuals can take only two values, so they cannot be normal, and their variance depends on the mean, peaking at p = 0.5.

Discussion

Many outcomes in linguistics are binary: correct or incorrect, one variant or another, a word recognised or not.

  • Which binary outcomes come up in your field?
  • Have you seen them analysed as percentages averaged per participant? What could go wrong with that?

The generalised linear model

Did you know: axolotls keep juvenile features such as feathery external gills all their lives (neoteny). Turning into a land salamander is a yes-or-no outcome that in the wild almost always comes out no.

Source: National Geographic

Data-generating processes

Winter (2019, section 12.2): think about the process that generates each data point.

In a normal model, each value is a mean plus some normally distributed noise. In a Bernoulli model, each response is a weighted coin that comes up 1 with probability p and 0 otherwise.

rbinom(10, size = 1, prob = 0.8)       # 10 responses, each correct with p = 0.8
 [1] 1 1 1 0 0 1 1 1 1 1

A logistic regression models how p changes with the predictors.

Three ingredients

A generalised linear model (GLM) has:

  1. A distribution for the outcome: here, the binomial (Bernoulli) distribution
  2. A linear predictor: b₀ + b₁ · x, exactly as before
  3. A link function, which connects the two: here, the logit

Why a link function? The linear predictor can take any value, from −∞ to +∞. A probability must stay between 0 and 1.

Probabilities, odds and log odds

tibble(p = c(0.1, 0.25, 0.5, 0.75, 0.9)) |>
  mutate(odds     = p / (1 - p),
         log_odds = log(odds)) |>
  print(n = 5)
# A tibble: 5 × 3
      p  odds log_odds
  <dbl> <dbl>    <dbl>
1  0.1  0.111    -2.20
2  0.25 0.333    -1.10
3  0.5  1         0   
4  0.75 3         1.10
5  0.9  9         2.20

The odds say how much more likely something is to happen than not: p = 0.75 gives odds of 3, or 3 to 1. Their log, the log odds or logit, runs from −∞ to +∞ and is 0 when p = 0.5.

The logistic function

The model is linear on the log odds scale. On the probability scale, it becomes an S-shaped curve that never goes below 0 or above 1.

In R

qlogis(0.75)      # probability -> log odds
[1] 1.098612
plogis(1.0986)    # log odds -> probability
[1] 0.7499977

qlogis() turns a probability into log odds, and plogis() turns log odds back into a probability.

Fitting with glmer()

Did you know: axolotls can regrow lost or damaged limbs, hearts and spinal cords, and even parts of their brain, usually without scarring.

Source: National Geographic

Which random effects?

The same reasoning as for lmer() (lectures 10 and 11):

  • Frequency varies within participants: (1 + freq_c | participant)
  • Frequency does not vary within words: (1 | word)
m <- glmer(correct ~ freq_c + (1 + freq_c | participant) + (1 | word),
           data = acc, family = binomial)

glmer() = lmer() + family. family = binomial tells R that the outcome is 0/1, and uses the logit link by default.

summary()

print(summary(m), correlation = FALSE, signif.stars = FALSE)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: correct ~ freq_c + (1 + freq_c | participant) + (1 | word)
   Data: acc

      AIC       BIC    logLik -2*log(L)  df.resid 
   1688.5    1721.5    -838.3    1676.5      1794 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.0720 -0.3731  0.3047  0.5124  2.7287 

Random effects:
 Groups      Name        Variance Std.Dev. Corr 
 word        (Intercept) 0.10693  0.3270        
 participant (Intercept) 0.60900  0.7804        
             freq_c      0.07684  0.2772   0.09 
Number of obs: 1800, groups:  word, 60; participant, 30

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.47695    0.16817   8.783   <2e-16
freq_c       0.85598    0.08612   9.939   <2e-16

Reading the output

Part of the output What it tells you
Laplace Approximation how the model was fitted: maximum likelihood, there is no REML
Random effects SDs on the log odds scale, and no residual variance: for 0/1 data the distribution fixes it
Fixed effects estimates on the log odds scale, with z values and p-values

The p-values are Wald tests: estimate / SE, compared with a normal distribution. They are approximate; likelihood ratio tests are usually preferred (below).

Interpreting coefficients

Did you know: wild axolotls live only in the lakes and canals of Xochimilco and Chalco, near Mexico City. They are critically endangered: one assessment found only 50 to 1,000 left in the wild.

Source: National Geographic

The intercept

fixef(m)
(Intercept)      freq_c 
  1.4769474   0.8559824

The intercept, 1.48, is the log odds of a correct response for a word of average frequency, by an average participant (truth: 1.5).

plogis(fixef(m)[1])
(Intercept) 
  0.8141111

On the probability scale: about 81% correct.

The slope: log odds

The slope, 0.86: each extra unit of log frequency adds 0.86 to the log odds of a correct response (truth: 0.8).

Log odds are hard to think in. Two ways out:

  1. Exponentiate, to get an odds ratio
  2. Compute predicted probabilities

Odds ratios

exp(fixef(m)[2])
  freq_c 
2.353686

Each extra unit of log frequency multiplies the odds of a correct response by about 2.35.

  • Odds ratio > 1: the outcome becomes more likely
  • Odds ratio < 1: less likely
  • Odds ratio = 1, that is a log odds slope of 0: no effect

Predicted probabilities

new <- tibble(freq_c = c(-2, -1, 0, 1, 2))
new |> mutate(p = predict(m, newdata = new, re.form = NA, type = "response")) |>
  print(n = 5)
# A tibble: 5 × 2
  freq_c     p
   <dbl> <dbl>
1     -2 0.442
2     -1 0.650
3      0 0.814
4      1 0.912
5      2 0.960

type = "response": probabilities, not log odds. re.form = NA: an average participant and word.

The same slope, different changes in probability: +21 points from −2 to −1, only +5 from 1 to 2.

The divide-by-4 rule

The logistic curve is steepest at p = 0.5, where its slope is exactly 1/4 of the log odds slope.

fixef(m)[2] / 4
   freq_c 
0.2139956

So one unit of log frequency changes the probability of a correct response by at most about 21 percentage points. Elsewhere on the curve, the change is smaller.

Predicted curves

Blue: each participant’s curve, with their own intercept and slope on the log odds scale. Orange: the average participant.

Testing and checking

Did you know: wild axolotls are mostly greyish-brown. The white and pink axolotls you see in pictures are usually bred as pets.

Source: National Geographic

Likelihood ratio tests

As in lecture 11: compare the model with and without the predictor.

m_null <- glmer(correct ~ 1 + (1 + freq_c | participant) + (1 | word),
                data = acc, family = binomial)
anova(m_null, m)
Data: acc
Models:
m_null: correct ~ 1 + (1 + freq_c | participant) + (1 | word)
m: correct ~ freq_c + (1 + freq_c | participant) + (1 | word)
       npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
m_null    5 1737.4 1764.9 -863.72    1727.4                         
m         6 1688.5 1721.5 -838.25    1676.5 50.936  1  9.545e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No REML = FALSE needed: glmer() always uses maximum likelihood.

Convergence

Mixed logistic models are harder to fit than lmer() models, and convergence warnings are common. The checklist from lecture 11 still applies:

  1. Is every random slope possible in the design?
  2. Centre or scale continuous predictors
  3. Try other optimisers, or more iterations: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
  4. Compare optimisers with allFit()
  5. Simplify the random effects, and report what you did

Binary data carry less information per observation than continuous data: complex random effects need a lot of data.

Reporting

Responses were more accurate for more frequent words: each unit of log frequency increased the log odds of a correct response by 0.86 (SE = 0.09; odds ratio = 2.35). A likelihood ratio test showed that frequency improved the model, χ²(1) = 50.9, p < .001. The model included by-participant random intercepts and slopes for frequency, and by-word random intercepts.

Report estimates on the log odds scale, and help your reader with odds ratios or predicted probabilities.

Counts and proportions

Sometimes data come as counts: 17 correct out of 20 trials. The binomial family handles this too, with successes and failures:

glmer(cbind(n_correct, n_incorrect) ~ predictor + (1 | participant),
      data = counts, family = binomial)

When you have the individual trials, model them directly, as we did today: you keep the information about words.

Your turn

Real data: accuracy in lexdec

library(languageR)
lexdec |> count(Correct)
    Correct    n
1   correct 1594
2 incorrect   65

Correct is a factor: correct or incorrect.

Exercises

  1. Create a column correct, coded 1 for correct responses and 0 otherwise
  2. Fit correct ~ Frequency + NativeLanguage with glmer(), by-subject and by-word random intercepts
  3. Interpret the Frequency effect as an odds ratio. What is the predicted probability of a correct response for a native and a non-native speaker, at the mean frequency?
  4. Test NativeLanguage with a likelihood ratio test

Careful: if you give glmer() the factor Correct directly, it models the probability of the second level, "incorrect".

Solutions

lex <- lexdec |> mutate(correct = as.numeric(Correct == "correct"))

m_lex <- glmer(correct ~ Frequency + NativeLanguage + (1 | Subject) + (1 | Word),
               data = lex, family = binomial)
fixef(m_lex)
        (Intercept)           Frequency NativeLanguageOther 
          1.7466051           0.5598952          -0.5725132 
exp(fixef(m_lex)["Frequency"])
Frequency 
 1.750489

Each unit of log frequency multiplies the odds of a correct response by about 1.75.

Solutions

new_lex <- tibble(Frequency = mean(lex$Frequency),
                  NativeLanguage = factor(c("English", "Other"), levels = levels(lex$NativeLanguage)))
predict(m_lex, newdata = new_lex, re.form = NA, type = "response")
        1         2 
0.9879519 0.9788392 
m_lex_null <- update(m_lex, . ~ . - NativeLanguage)
anova(m_lex_null, m_lex)
Data: lex
Models:
m_lex_null: correct ~ Frequency + (1 | Subject) + (1 | Word)
m_lex: correct ~ Frequency + NativeLanguage + (1 | Subject) + (1 | Word)
           npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
m_lex_null    4 496.65 518.31 -244.33    488.65                     
m_lex         5 497.21 524.28 -243.60    487.21 1.4443  1     0.2294

Both groups are very accurate, and native language does not reliably improve the model.

Summary

  • Binary outcomes need glmer(..., family = binomial) instead of lmer()
  • A GLM: a distribution, a linear predictor, a link function
  • Coefficients are on the log odds scale: exp() gives odds ratios, plogis() probabilities
  • The same slope means different changes in probability along the S-curve; divide by 4 for the maximum
  • Test with likelihood ratio tests; expect more convergence issues
  • The random effects logic is the same as for lmer()

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ws26_statistics/12_logistic_regression.html

Code, data and all lectures

github.com/mcanzi/ws26_statistics