A common alternative to treatment coding is called sum coding.
It’s easiest to see how it works using examples.
Research Question
Does the sepal length of particular iris species differ significantly from the average sepal length?
Example: Binary predictor
You have been provided with the in-built dataset, iris, which contains information concerning the sepal length (in cm), sepal width (in cm), petal length (in cm), and petal width (in cm) from three different species of iris (setosa, versicolor, and virginica). There are measurements for 50 flowers from each of the iris species (i.e., total \(n\) = 150). For this example, we’ll just look at two species, setosa and versicolor.
iris_binary <- iris |># filter out virginicafilter(Species !='virginica') |># the only factor levels we want are the two species we're keepingmutate(Species =factor(Species, levels =c('setosa', 'versicolor')))head(iris_binary)
(If this code doesn’t run for your variable, then you’ll probably need to convert your variable to a factor using factor()).
This contrast matrix shows us that setosa is represented as 1, and versicolor is represented as –1. These two values are centered around 0, so 0 will represent the mean, the average, between these two species.
When we fit a linear model to this data, the model’s intercept will represent the grand mean sepal length: that is, the mean of both species’ mean sepal lengths. And the model’s slope over Species will represent the estimated difference between the mean sepal length for setosa (the level coded as 1) and the grand mean. Or, in other words, how sepal length changes when we move from the grand mean to the mean of setosa.
We can give the column of this matrix a more useful name with the following code:
Read this contrast matrix as a table with three rows and two columns.
The first column shows the first variable, which will estimate the difference between the mean sepal length of setosa (coded as 1 in that column) and the grand mean sepal length of all three species. The second column shows the second variable, which will estimate the difference between the mean sepal length of versicolor (coded as 1 in that column) and the grand mean sepal length of all three species.
When we use Species as a predictor in a model, each of these variables will receive its own coefficient in the model summary. That coefficient represents the comparisons described for each variable.
We can give the columns of this matrix more useful names with the following code:
# A tibble: 3 x 2
Species m
<fct> <dbl>
1 setosa 5.01
2 versicolor 5.94
3 virginica 6.59
mean(c(5.01, 5.94, 6.59))
[1] 5.85
The slope coefficient SpeciesSetosa is –0.84, which corresponds to the difference between the setosa mean of 5.01 and the grand mean of 5.85.
5.01-5.85
[1] -0.84
And the slope coefficient SpeciesVersicolor is 0.09, which corresponds to the difference between the versicolor mean of 5.94 and the grand mean of 5.85.
5.94-5.85
[1] 0.09
Plotting model-fitted values
There are a couple of ways that we can visualise a model of categorical predictors, either using the sjPlot or effects packages:
library(sjPlot)plot_model(m2,type ="eff",terms ="Species") +labs(title ="Sepal Length by Species",x ="Species", y ="Sepal Length")
library(effects)effect(term =c("Species"), mod = m2) |>as.data.frame() |>ggplot(aes(x = Species, y = fit, col = Species)) +geom_pointrange(aes(ymin = lower, ymax = upper))
Reporting sum coding (three-level example)
Variables
$$
\text{Species1} = \begin{cases}
1 & \text{if Species is Setosa} \\
0 & \text{if Species is Versicolor} \\
-1 & \text{if Species is Virginica}
\end{cases}
$$
\[
\text{Species1} = \begin{cases}
1 & \text{if Species is Setosa} \\
0 & \text{if Species is Versicolor} \\
-1 & \text{if Species is Virginica}
\end{cases}
\]
$$
\text{Species2} = \begin{cases}
0 & \text{if Species is Setosa} \\
1 & \text{if Species is Versicolor} \\
-1 & \text{if Species is Virginica}
\end{cases}
$$
\[
\text{Species2} = \begin{cases}
0 & \text{if Species is Setosa} \\
1 & \text{if Species is Versicolor} \\
-1 & \text{if Species is Virginica}
\end{cases}
\]
For all coefficients, linear models always test the null hypothesis that the coefficient is equal to zero.
More specifically, for the mathematical model specified above:
For the intercept: the null hypothesis H0 is that the grand mean outcome is equal to zero. The alternative hypothesis H1 is that the grand mean outcome is different from zero. (This isn’t usually interesting to test, because the intercept will nearly always be different from zero.)
For a single slope coefficient (let’s say \(\beta_1\)): the null hypothesis is that the difference between the \(\beta_1\)’s level coded as 1 and the grand mean is equal to zero. And the alternative hypothesis is thatthe difference between the \(\beta_1\)’s level coded as 1 and the grand mean is different from zero.
For all slope coefficients: the null hypothesis is that the differences between all levels coded as 1 and the grand mean are equal to zero. And the alternative hypothesis is that any difference between any level coded as 1 and the grand mean is different from zero.
Choose your contrast coding scheme based on the hypotheses you want to test.
For example: If your RQ asks whether particular groups differ from the grand mean, this question can be tested by looking at a slope coefficient of a sum-coded categorical predictor where the particular group you’re interested in is coded as 1.
Changing which levels are coded as 1
Imagine if we wanted the Species1 coefficient to compare the grand mean to virginica instead of setosa.
To do this, reorder the factor levels using factor(). The level which appears last will not be compared to the grand mean by any of the regression coefficients.
Beware that some people will also use these labels to refer to a similar contrast coding scheme which uses –0.5 and 0.5 instead of –1 and 1. For that reason, it’s best to explicitly report what numbers you’re using for which levels. That way, the reader knows exactly what you’re doing.
