In the ANOVA and the GLM chapter we let R turn a factor into numbers for us and moved on. That was a simplification worth revisiting. R always codes a categorical predictor into numbers before it fits anything, and the code you choose decides what every coefficient means. Change the code and the intercept, the slopes, and the t-tests all change their meaning, even though the model’s fit (its \(R^2\), its \(F\), its predicted values) does not change at all.
That is the whole chapter in one idea. The coding scheme is a lens, not a model. Pick the lens that answers your question, and the coefficients read like plain English. Pick a careless lens and you get numbers that are technically correct but not meaningful, what we will call nonsense coding.
lm() output is ever a mystery.Let’s build a categorical variable from scratch so we can see the seams. Take a continuous “dogness” score, chop it into four bins, and call them cat, fish, pig, dog:
set.seed(725)
n <- 200
x <- rnorm(n) # a continuous "dogness" score
pet <- cut(x, 4, labels = c("cat", "fish", "pig", "dog"))
happy <- x + rnorm(n) # more dogness -> happier, on average
d <- data.frame(pet, happy)
table(pet)pet
cat fish pig dog
9 60 94 37 pet is now a factor: R stores it as integer codes with labels attached. The moment we put it in a model, R replaces those four labels with three columns of numbers (four groups need \(4-1=3\) columns, plus the intercept). Which three columns is the coding scheme. Here are the group means we will spend the rest of the chapter reconstructing:
means <- aggregate(happy ~ pet, data = d, FUN = mean)
means| pet | happy |
|---|---|
| cat | -2.539008 |
| fish | -1.159632 |
| pig | 0.371500 |
| dog | 1.196379 |
* Group means of happy by pet.
MEANS TABLES=happy BY pet
/CELLS=MEAN STDDEV COUNT.using Statistics, DataFrames
combine(groupby(df, :pet), :happy => mean)Dummy coding - R’s default, also called reference or indicator coding - picks one group as the reference and compares every other group to it. Let’s make dog the reference (the 4th level):
contrasts(d$pet) <- contr.treatment(4, base = 4) # reference = level 4 = dog
contrasts(d$pet) 1 2 3
cat 1 0 0
fish 0 1 0
pig 0 0 1
dog 0 0 0m_dummy <- lm(happy ~ pet, data = d)
round(coef(m_dummy), 3)(Intercept) pet1 pet2 pet3
1.196 -3.735 -2.356 -0.825 * Dummy / reference coding, reference = dog (4th category).
UNIANOVA happy BY pet
/CONTRAST(pet)=SIMPLE(4)
/PRINT=PARAMETER.using GLM, StatsModels
m_dummy = lm(@formula(happy ~ pet), df,
contrasts = Dict(:pet => DummyCoding(base = "dog")))Read the coefficients straight off the means:
(Intercept) is the reference-group mean - the mean happiness of dog.pet1 is \(\bar{y}_{cat} - \bar{y}_{dog}\), pet2 is \(\bar{y}_{fish} - \bar{y}_{dog}\), and so on.Don’t take my word for it - check:
mu <- setNames(means$happy, means$pet) # named group means
c(intercept_is_dog_mean = all.equal(unname(coef(m_dummy)[1]), unname(mu["dog"])),
cat_slope_is_cat_minus_dog = all.equal(unname(coef(m_dummy)["pet1"]),
unname(mu["cat"] - mu["dog"]))) intercept_is_dog_mean cat_slope_is_cat_minus_dog
TRUE TRUE Both TRUE. Every dummy coefficient is a simple difference from the reference. That is exactly why dummy coding is the right lens when you have a control group and want each treatment compared to it.
Add a second predictor with an interaction and a dummy slope is no longer “group mean minus reference mean” overall - it is that difference only at the reference level of the other predictor (i.e., where the other predictor’s dummy is 0). This is the single most common misreading of regression output. When an interaction is in the model, stop reading coefficients off marginal means and solve the full equation instead. (This bites hardest in the factorial models of ANOVA & the GLM.)
Dummy coding asks “how does each group differ from one group?” Effects coding (a.k.a. deviation coding, contr.sum) asks a different question: “how does each group differ from the grand mean?” The zero point moves from a reference group to the center of everything.
contrasts(d$pet) <- contr.sum(4)
contrasts(d$pet) [,1] [,2] [,3]
cat 1 0 0
fish 0 1 0
pig 0 0 1
dog -1 -1 -1m_eff <- lm(happy ~ pet, data = d)
round(coef(m_eff), 3)(Intercept) pet1 pet2 pet3
-0.533 -2.006 -0.627 0.904 * Effects / deviation coding (each group vs the grand mean).
UNIANOVA happy BY pet
/CONTRAST(pet)=DEVIATION
/PRINT=PARAMETER.using GLM, StatsModels
m_eff = lm(@formula(happy ~ pet), df,
contrasts = Dict(:pet => EffectsCoding()))Now the meanings shift:
(Intercept) is the grand mean - but note which grand mean: the mean of the four group means, not the mean of the raw data (those differ when groups are unequal in size).dog) has no coefficient of its own. Because effects sum to zero, its effect is minus the sum of the others.Watch the grand-mean subtlety, because it trips everyone:
GM_groups <- mean(means$happy) # mean of the GROUP means
GM_raw <- mean(d$happy) # mean of the raw data
c(intercept = unname(coef(m_eff)[1]),
grand_mean_of_groups = GM_groups, # this one matches
grand_mean_of_raw = GM_raw) # this one (slightly) does not intercept grand_mean_of_groups grand_mean_of_raw
-0.53269001 -0.53269001 -0.06620966 And recover the invisible fourth effect - dog - as minus the sum of the rest:
eff <- coef(m_eff)[-1] # pet1, pet2, pet3 = cat, fish, pig
dog_effect <- -sum(eff)
all.equal(unname(mu["dog"] - GM_groups), unname(dog_effect))[1] TRUETRUE. That is the whole trick of effects coding: the coefficients are deviations that must sum to zero, so the omitted group is always recoverable.
Dummy codes make the zero explicit (it is the reference group); effects codes make the zero implicit (it is the grand mean, sitting between the deviations). That is the mnemonic to carry out of this chapter.
contr.sum on a two-level factor assigns +1 to the first level and -1 to the second, which is often the opposite of what you expect, and it silently flips the sign of your effect. R’s default level ordering (alphabetical) can catch you the same way. Print contrasts(yourfactor) before you interpret anything.
Plain effects coding centers on the mean of the group means, treating a group of 10 and a group of 200 as equals. When your cells are unbalanced and you want the population grand mean (weighted by how many people are actually in each group), you need weighted effects coding.
Let’s make an unbalanced world - many pigs, few dogs - and build the weighted codes by hand, pushing the cell-size ratios onto the comparison group’s row:
set.seed(9)
cnt <- c(cat = 20, fish = 30, pig = 40, dog = 10)
petw <- factor(rep(names(cnt), cnt), levels = names(cnt))
# group-specific means so the demo has structure:
mu_w <- c(cat = 0, fish = 1, pig = 2, dog = 3)
happyw <- rnorm(length(petw), mean = mu_w[as.character(petw)], sd = 1)
contrasts(petw) <- contr.sum(4) # start from unweighted effects codes
M <- contrasts(petw)
N <- table(petw)
# reweight the comparison (last) group's row by the cell-size ratios:
M[4, ] <- c(M[4, 1] * N[1] / N[4],
M[4, 2] * N[2] / N[4],
M[4, 3] * N[3] / N[4])
contrasts(petw) <- M
round(M, 2) [,1] [,2] [,3]
cat 1 0 0
fish 0 1 0
pig 0 0 1
dog -2 -3 -4m_weff <- lm(happyw ~ petw)The payoff: the intercept is now the sample grand mean - the ordinary mean of the raw data, weighted by cell size - not the mean of the group means:
c(intercept = unname(coef(m_weff)[1]),
raw_mean = mean(happyw), # weighted effects intercept matches THIS
groups_mean = mean(tapply(happyw, petw, mean))) intercept raw_mean groups_mean
1.346485 1.346485 1.390092 The intercept tracks the raw mean. The weights get pushed onto the comparison-group means, so unbalanced data no longer distorts your grand mean.
Dummy and effects codes answer fixed questions. Contrast codes let you pose the question - “is the average of these two groups different from the average of those two?” - by writing the comparison directly into the code columns. Done right (orthogonally), each coefficient is one clean, independent comparison.
R ships one useful family, Helmert codes, which compare each level to the pooled mean of the levels before it:
contrasts(d$pet) <- contr.helmert(4)
m_helm <- lm(happy ~ pet, data = d)
round(coef(m_helm), 3)(Intercept) pet1 pet2 pet3
-0.533 0.690 0.740 0.576 * Helmert contrasts.
UNIANOVA happy BY pet
/CONTRAST(pet)=HELMERT
/PRINT=PARAMETER.using GLM, StatsModels
m_helm = lm(@formula(happy ~ pet), df,
contrasts = Dict(:pet => HelmertCoding()))The intercept is again the (unweighted) grand mean. The first Helmert slope is half the gap between the first two group means - why half? Because the code uses \(-1\) and \(+1\), a two-unit spread, so the slope is the difference divided by 2:
c(helmert_b1 = unname(coef(m_helm)[2]),
half_the_gap = unname((mu["fish"] - mu["cat"]) / 2)) helmert_b1 half_the_gap
0.689688 0.689688 The real power is writing your own. Suppose the scientific questions are: (1) are the “usual” pets (cat, fish) different from the “unusual” ones (pig, dog)? (2) cat vs. fish? (3) pig vs. dog? Encode exactly those, and name the columns so the output documents itself:
sgc <- matrix(c( 0.5, 0.5, -0.5, -0.5, # usual (cat,fish) vs unusual (pig,dog)
0.5, -0.5, 0.0, 0.0, # cat vs fish
0.0, 0.0, 0.5, -0.5), # pig vs dog
nrow = 4,
dimnames = list(levels(d$pet),
c("usual_vs_unusual", "cat_vs_fish", "pig_vs_dog")))
sgc usual_vs_unusual cat_vs_fish pig_vs_dog
cat 0.5 0.5 0.0
fish 0.5 -0.5 0.0
pig -0.5 0.0 0.5
dog -0.5 0.0 -0.5contrasts(d$pet) <- sgc
m_con <- lm(happy ~ pet, data = d)
round(coef(m_con), 3) (Intercept) petusual_vs_unusual petcat_vs_fish petpig_vs_dog
-0.533 -2.633 -1.379 -0.825 * Custom orthogonal contrasts:
* usual (cat,fish) vs unusual (pig,dog); cat vs fish; pig vs dog.
UNIANOVA happy BY pet
/CONTRAST(pet)=SPECIAL(0.5 0.5 -0.5 -0.5 0.5 -0.5 0 0 0 0 0.5 -0.5)
/PRINT=PARAMETER.using GLM, StatsModels
m_con = lm(@formula(happy ~ pet), df,
contrasts = Dict(:pet => HypothesisCoding(
[0.5 0.5 -0.5 -0.5; 0.5 -0.5 0 0; 0 0 0.5 -0.5])))These three contrasts are orthogonal (each pair of columns is uncorrelated), so the three comparisons are independent - each coefficient answers its question without contamination from the others. Confirm the first one really is the usual-vs-unusual gap:
usual <- mean(mu[c("cat", "fish")])
unusual <- mean(mu[c("pig", "dog")])
c(contrast_b1 = unname(coef(m_con)["petusual_vs_unusual"]),
usual_minus_unusual = usual - unusual) contrast_b1 usual_minus_unusual
-2.633259 -2.633259 Effects and contrast codes must be centered to move between means and coefficients - so you divide by the right denominator (the spread of the codes) when you compute a \(b\) from marginal means, and multiply back when you go the other way. Miss the denominator and your hand-check will be “off by a factor” and you’ll think you broke something. You didn’t; you forgot to center.
Now the capstone, and the point of the whole chapter. What if we assign arbitrary numbers to the groups, with no reference, no deviations, and no meaningful comparison behind them?
nonsense <- matrix(c( 3.0, -7.0, 2.5, 0.1,
1.0, 1.0, -2.0, 4.0,
-5.0, 2.0, 8.0, -1.0),
nrow = 4,
dimnames = list(levels(d$pet), c("huh", "wat", "lol")))
contrasts(d$pet) <- nonsense
m_nonsense <- lm(happy ~ pet, data = d)
round(coef(m_nonsense), 3)(Intercept) pethuh petwat petlol
-1.734 0.158 0.834 0.422 * Nonsense coding: arbitrary contrast coefficients with no meaning.
UNIANOVA happy BY pet
/CONTRAST(pet)=SPECIAL(3 -7 2.5 0.1 1 1 -2 4 -5 2 8 -1)
/PRINT=PARAMETER.using GLM, StatsModels
m_nonsense = lm(@formula(happy ~ pet), df,
contrasts = Dict(:pet => HypothesisCoding(
[3 -7 2.5 0.1; 1 1 -2 4; -5 2 8 -1])))Those coefficients are not interpretable, which is exactly why the columns are named huh, wat, and lol. You cannot read a single one off the group means. And yet look at what has not changed:
# The model fit is byte-for-byte identical to the sensible codings:
c(same_fitted = isTRUE(all.equal(fitted(m_nonsense), fitted(m_dummy))),
same_r2 = all.equal(summary(m_nonsense)$r.squared, summary(m_dummy)$r.squared),
same_F = all.equal(summary(m_nonsense)$fstatistic[1],
summary(m_dummy)$fstatistic[1]))same_fitted same_r2 same_F
TRUE TRUE TRUE All TRUE. Every coding scheme in this chapter fits the exact same model - same predicted values, same \(R^2\), same omnibus \(F\). Dummy, effects, weighted, contrast, and nonsense codes are five different descriptions of one fitted surface. The omnibus test (“does pet matter at all?”) is blind to coding. Only the individual coefficients - the intercept and the \(b\)’s - carry the fingerprint of the code you chose.
That is the lesson worth holding onto:
The data decide the fit. You decide, through the coding scheme, which questions the coefficients answer. Choose codes whose contrasts match your scientific question and the output reads like a sentence. Choose thoughtlessly and you get nonsense coding: correct arithmetic, meaningless numbers.
| Scheme | R helper | Intercept means | A slope means | Use it when |
|---|---|---|---|---|
| Dummy / reference | contr.treatment |
reference-group mean | group − reference | you have a control/baseline group |
| Effects / deviation | contr.sum |
grand mean (of group means) | group − grand mean | groups are balanced; you want deviations |
| Weighted effects | hand-built from contr.sum |
sample grand mean (weighted) | group − weighted grand mean | groups are unbalanced |
| Contrast (Helmert/custom) | contr.helmert / your matrix |
grand mean | one planned comparison | you have specific, orthogonal questions |
| Nonsense | arbitrary matrix | nothing interpretable | nothing interpretable | never (it’s the cautionary tale) |
Across every row the \(R^2\), \(F\), and fitted values are identical. The columns “Intercept means” and “A slope means” are the only things that move.
Simulated pets are friendly, but let’s prove the machinery on real data - the Module 3 dataset from the graduate GLM course (data/mod3data.csv). It has an outcome y1 and a genuine four-level grouping factor, f3, with unequal cells (real data is always lumpy):
d3 <- read.csv("data/mod3data.csv")
d3$grp <- factor(d3$f3, labels = c("A", "B", "C", "D"))
mu3 <- tapply(d3$y1, d3$grp, mean)
round(mu3, 3) # the four real group means A B C D
0.506 -0.334 -0.132 -0.019 table(d3$grp) # unequal cell sizes
A B C D
5 7 11 7 * Group means and cell sizes of y1 by grp.
MEANS TABLES=y1 BY grp
/CELLS=MEAN STDDEV COUNT.using Statistics, DataFrames
combine(groupby(df, :grp), :y1 => mean)Dummy-code it (reference = group A) and read the coefficients straight off those means, exactly as we did for the pets:
contrasts(d3$grp) <- contr.treatment(4, base = 1)
m_real <- lm(y1 ~ grp, data = d3)
round(coef(m_real), 3)(Intercept) grp2 grp3 grp4
0.506 -0.840 -0.638 -0.525 c(intercept_is_A = all.equal(unname(coef(m_real)[1]), unname(mu3["A"])),
grp2_is_B_minus_A = all.equal(unname(coef(m_real)["grp2"]), unname(mu3["B"] - mu3["A"]))) intercept_is_A grp2_is_B_minus_A
TRUE TRUE * Dummy / reference coding, reference = group A (1st category).
UNIANOVA y1 BY grp
/CONTRAST(grp)=SIMPLE(1)
/PRINT=PARAMETER.using GLM, StatsModels
m_real = lm(@formula(y1 ~ grp), df,
contrasts = Dict(:grp => DummyCoding(base = "A")))Both TRUE - the intercept is group A’s mean, each slope a group’s mean minus A’s. Switch to effects coding and the intercept becomes the grand mean of the four group means, unequal cells and all:
contrasts(d3$grp) <- contr.sum(4)
m_real_eff <- lm(y1 ~ grp, data = d3)
all.equal(unname(coef(m_real_eff)[1]), mean(mu3))[1] TRUE* Effects / deviation coding (each group vs the grand mean).
UNIANOVA y1 BY grp
/CONTRAST(grp)=DEVIATION
/PRINT=PARAMETER.using GLM, StatsModels
m_real_eff = lm(@formula(y1 ~ grp), df,
contrasts = Dict(:grp => EffectsCoding()))TRUE. Real data, real (unbalanced) cells, and every rule from the pet example still holds, because the arithmetic works the same whether the groups are house pets or the conditions of an actual study.
contr.treatment using cat as the reference instead of dog. Which coefficients change sign, and why is \(R^2\) untouched?crossprod(yourmatrix) - the off-diagonals should be 0). Fit it, and reconstruct one coefficient from the group means, remembering to divide by the code spread.Coding is where categorical data enter the linear model, and where careful thought pays off in coefficients you can actually explain. It is also the last piece of the observed-variable linear model. But the constructs we most want to study, such as depression, curiosity, and intelligence, are never observed directly; we see only their noisy indicators. The next part turns to latent variables and the factor and structural models built for them, and after that, to what you do when the normal curve no longer describes your data well.