Categorical predictors

Treatment coding, sum coding and factors

Dr Massimiliano Canzi

Categorical predictors

Did you know: mantis shrimp have 12 types of photoreceptor for colour in their eyes. Humans have three.

Source: Thoen et al. (2014), Science

Today

  • What a categorical predictor is, and how R stores it: factors
  • Treatment coding: the default in R
  • Dummy coding by hand, and changing the reference level
  • Sum coding
  • Predictors with more than two levels
  • Combining categorical and continuous predictors

So far

Until now every predictor was a number: frequency, length, body mass. Each slope was a change in y per unit of x.

Many linguistic predictors are categories instead: voicing (voiced, voiceless), speaking style (casual, clear, loud), language, condition, speaker group, word class.

What does “one unit more” mean for voicing?

The data

Vowels are longer before voiced consonants (bead) than before voiceless ones (beat). A simulated production study:

set.seed(7)
vowels <- tibble(
  voicing  = rep(c("voiceless", "voiced"), each = 60),
  duration = round(c(rnorm(60, 120, 20), rnorm(60, 160, 22)))   # ms
)
vowels |> summarise(mean = mean(duration), sd = sd(duration), n = n(), .by = voicing)
# A tibble: 2 × 4
  voicing    mean    sd     n
  <chr>     <dbl> <dbl> <int>
1 voiceless  124.  19.7    60
2 voiced     162.  19.8    60

The data

White: the mean of each group. A model with voicing as a predictor should recover these two means.

Treatment coding

Did you know: mantis shrimp come in two kinds: smashers break open crabs, clams and snails with club-like limbs, and spearers ambush fish from their burrows and impale them.

Source: UC Berkeley News

lm() with a categorical predictor

m <- lm(duration ~ voicing, data = vowels)
coef(m)
     (Intercept) voicingvoiceless 
       162.06667        -37.76667

Compare with the group means:

vowels |> summarise(mean = mean(duration), .by = voicing)
# A tibble: 2 × 2
  voicing    mean
  <chr>     <dbl>
1 voiceless  124.
2 voiced     162.
  • The intercept is the mean of voiced
  • The slope voicingvoiceless is the difference: voiceless minus voiced

The reference level

Why is voiced the intercept?

R turns text into a factor, and puts its levels in alphabetical order:

levels(factor(vowels$voicing))
[1] "voiced"    "voiceless"

The first level becomes the reference level. Every other level is compared against it.

R names the coefficient voicingvoiceless: the predictor, then the level being compared with the reference.

Treatment coding, visually

Code the reference level as 0 and the other as 1: the regression line goes through the two means. A one-unit change in x is the change from one category to the other.

Dummy coding by hand

This is exactly what R does behind the scenes:

vowels <- vowels |>
  mutate(voiceless01 = if_else(voicing == "voiceless", 1, 0))

coef(lm(duration ~ voiceless01, data = vowels))
(Intercept) voiceless01 
  162.06667   -37.76667

Same numbers as before. A variable coded 0 and 1 like this is called a dummy variable.

Changing the reference level

Voiceless is the more natural baseline here. Set the order of the levels with fct_relevel() from {forcats}:

vowels <- vowels |>
  mutate(voicing = fct_relevel(voicing, "voiceless"))

levels(vowels$voicing)
[1] "voiceless" "voiced"   
coef(lm(duration ~ voicing, data = vowels))
  (Intercept) voicingvoiced 
    124.30000      37.76667

The intercept is now the voiceless mean, and the slope is +38 ms for voiced. The model and its fit are unchanged; only the comparison is different.

Choosing the reference level

  • Choose the level that makes the coefficients easiest to read, often a natural baseline: a control condition, the standard variety, the most frequent category
  • Set it explicitly, rather than relying on the alphabet

In R, the reference level of a factor is simply its first level.

Discussion

Your study compares L1 English speakers with three groups of L2 speakers of English, whose first languages are German, Italian and Mandarin.

  • Which group would you choose as the reference level, and why?
  • Would your choice change if your question were about differences among the three L2 groups?

Sum coding

Did you know: the peacock mantis shrimp strikes snail shells at up to 23 m/s (about 83 km/h), with accelerations of up to 8,000 times gravity. The strike also makes bubbles that collapse with destructive force.

Sources: Patek, Korff & Caldwell (2004), Nature; UC Berkeley News

Sum coding

Instead of 0 and 1, code the two levels as +1 and −1:

contrasts(vowels$voicing) <- contr.sum(2)
contrasts(vowels$voicing)
          [,1]
voiceless    1
voiced      -1
coef(lm(duration ~ voicing, data = vowels))
(Intercept)    voicing1 
  143.18333   -18.88333

Reading sum-coded coefficients

means <- vowels |> summarise(mean = mean(duration), .by = voicing)
mean(means$mean)                      # the mean of the group means
[1] 143.1833
diff(means$mean) / 2                  # half the difference
[1] 18.88333
  • The intercept is the grand mean: the average of the group means
  • The coefficient is how far each group lies from that grand mean: + for the first level, − for the second

Sum coding, visually

Arrows: the sum-coded coefficient, ± 18.9 ms from the grand mean.

Treatment or sum coding?

Treatment (0 / 1) Sum (+1 / −1)
Intercept mean of the reference level grand mean
Coefficient difference from the reference deviation from the grand mean
Default in R yes no: contr.sum()

The fit is identical. The choice changes what the intercept means, which matters most when you add interactions (next week).

More than two levels

Did you know: despite their 12 receptors, mantis shrimp were worse than humans, bees and butterflies at telling similar colours apart. They seem to encode colour unlike any other animal studied: perhaps by recognising colours rather than comparing them.

Sources: Thoen et al. (2014), Science; Sci.News

Speaking style

A simulated study of vowel duration in three speaking styles:

set.seed(9)
styles <- tibble(
  style    = rep(c("casual", "clear", "loud"), each = 50),
  duration = round(c(rnorm(50, 110, 18), rnorm(50, 150, 20), rnorm(50, 132, 19)))
)
styles |> summarise(mean = mean(duration), .by = style)
# A tibble: 3 × 2
  style   mean
  <chr>  <dbl>
1 casual  108.
2 clear   150.
3 loud    130.

Speaking style

Three levels, two coefficients

m_style <- lm(duration ~ style, data = styles)
coef(m_style)
(Intercept)  styleclear   styleloud 
     108.34       41.40       21.68

The intercept is the mean of casual, the reference level. styleclear is clear minus casual, and styleloud is loud minus casual.

A predictor with k levels gets k − 1 coefficients: each level against the reference.

Three levels, visually

Each arrow is one coefficient: how far clear and loud speech lie from casual speech, the dashed line.

Behind the scenes

R creates one dummy variable per non-reference level:

model.matrix(~ style, data = styles)[c(1, 51, 101), ]
    (Intercept) styleclear styleloud
1             1          0         0
51            1          1         0
101           1          0         1
style styleclear styleloud predicted duration
casual 0 0 b₀
clear 1 0 b₀ + b₁
loud 0 1 b₀ + b₂

What is missing?

The model does not directly compare clear with loud.

Change the reference level to get that comparison:

styles |>
  mutate(style = fct_relevel(style, "clear")) |>
  lm(duration ~ style, data = _) |>
  coef()
(Intercept) stylecasual   styleloud 
     149.74      -41.40      -19.72

data = _ uses the native pipe’s placeholder: the data come from the pipe.

Sum coding with three levels

styles <- styles |> mutate(style = factor(style))
contrasts(styles$style) <- contr.sum(3)
coef(lm(duration ~ style, data = styles))
(Intercept)      style1      style2 
  129.36667   -21.02667    20.37333

The intercept is now the grand mean of the three styles. style1 is casual’s deviation from it and style2 is clear’s. Loud’s deviation is not shown, because it is minus the sum of the other two.

Model fit and assumptions

glance(m_style)[, 1:2]
# A tibble: 1 × 2
  r.squared adj.r.squared
      <dbl>         <dbl>
1     0.459         0.451

The same assumptions as before. With categories, check that the residuals have a similar spread in every group.

Mixing predictors

Did you know: mantis shrimp were the first animals shown to see circularly polarised light. Some species also reflect it from body parts they display in courtship.

Source: Chiou et al. (2008), Current Biology

Categorical and continuous together

Back to the warning from lecture 4: bill length by bill depth.

bills <- drop_na(penguins, bill_length_mm, bill_depth_mm)
coef(lm(bill_length_mm ~ bill_depth_mm, data = bills))
  (Intercept) bill_depth_mm 
   55.0673698    -0.6498356

Add species as a second predictor:

m_bills <- lm(bill_length_mm ~ bill_depth_mm + species, data = bills)
coef(m_bills)
     (Intercept)    bill_depth_mm speciesChinstrap    speciesGentoo 
       13.216370         1.394011         9.938955        13.403279

With species in the model, the bill depth slope turns positive, as it is within each species.

Three parallel lines

One slope for bill depth, and a different intercept for each species (Adelie is the reference). The lines are parallel: no interaction.

Reading the mixed model

tidy(m_bills)[, 1:2]
# A tibble: 4 × 2
  term             estimate
  <chr>               <dbl>
1 (Intercept)         13.2 
2 bill_depth_mm        1.39
3 speciesChinstrap     9.94
4 speciesGentoo       13.4
Coefficient Meaning
bill_depth_mm 1.39 mm of length per mm of depth, within a species
speciesChinstrap Chinstrap bills are 9.9 mm longer than Adelie bills of the same depth
speciesGentoo Gentoo bills, 13.4 mm longer

Other coding schemes

Treatment and sum coding are the two you will meet most. There are others, for example:

  • Helmert coding: each level against the mean of the previous levels
  • Successive differences: each level against the previous one, for ordered categories

All of them give the same fit. They differ only in which comparisons the coefficients show: choose the one that matches your research question.

Your turn

Exercises

With penguins:

  1. Model body_mass_g by species. Check that the coefficients match the group means.
  2. Make Gentoo the reference level and refit. What changed?
  3. Refit with sum coding. What is the intercept now?
  4. Model body_mass_g by sex and flipper_length_mm together. Interpret both coefficients.

Solutions

pengs <- drop_na(penguins, body_mass_g, sex, flipper_length_mm)

coef(lm(body_mass_g ~ species, data = pengs))
     (Intercept) speciesChinstrap    speciesGentoo 
      3706.16438         26.92385       1386.27259 
pengs |> summarise(mean = mean(body_mass_g), .by = species)
# A tibble: 3 × 2
  species    mean
  <fct>     <dbl>
1 Adelie    3706.
2 Gentoo    5092.
3 Chinstrap 3733.
pengs |>
  mutate(species = fct_relevel(species, "Gentoo")) |>
  lm(body_mass_g ~ species, data = _) |>
  coef()
     (Intercept)    speciesAdelie speciesChinstrap 
        5092.437        -1386.273        -1359.349

Solutions

pengs_sum <- pengs
contrasts(pengs_sum$species) <- contr.sum(3)
coef(lm(body_mass_g ~ species, data = pengs_sum))
(Intercept)    species1    species2 
  4177.2299   -471.0655   -444.1416 
coef(lm(body_mass_g ~ sex + flipper_length_mm, data = pengs))
      (Intercept)           sexmale flipper_length_mm 
      -5410.30022         347.85025          46.98218

Males are about 348 g heavier than females with the same flipper length.

Summary

  • R stores categories as factors; the first level is the reference
  • Treatment coding compares each level with the reference; sum coding compares each with the grand mean
  • k levels → k − 1 coefficients; change the reference to see other comparisons
  • Categorical and continuous predictors can be combined: one slope, several intercepts

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ws26_statistics/07_categorical_predictors.html

Code, data and all lectures

github.com/mcanzi/ws26_statistics