knitr::opts_chunk$set( collapse = TRUE, comment = "#>", out.width = "100%", fig.width = 8, fig.height = 4 ) # run a simple analysis ctable < function(df, predictor = "pet") { a < paste("y ~", predictor) %>% formula() %>% lm(df) contrast_table < kable(contrasts(df[[predictor]]), format = "html") %>% kable_styling() analysis_table < kable(broom::tidy(a), format = "html") %>% kable_styling() paste( "<table class='ctable table'>", "<thead><tr><th>Contrasts</th><th>lm(y ~ pet, df)</th></thead>", "<tbody><tr><td>", contrast_table, "</td><td>", analysis_table, "</td></tr></tbody></table>", sep = "\n" ) } options(digits = 3) set.seed(8675309) library(faux) library(kableExtra)
Contrasts other than treatment coding for factors with more than two levels have always confused me. I try to design all my own studies to never have more than two levels per categorical factor, but students in my data skills class are always trying to analyse data with threelevel categories (or worse).
I thought it was just me, but some recent Twitter discussion showed me I'm not alone. There's a serious jinglejangle problem with contrasts. Here are some of the explainers that I found useful (thanks to everyone who recommended them).
Contrasts often have multiple names. The names I'm using try to maintain the relationship with the base R function, apart from anova coding, which was suggested by Dale Barr after I got frustrated that people use so many different labels for that (extremely useful) coding and each of the terms used is also used to refer to totally different codings by others.
 My name  Other names  R function  faux function 
::::
 Treatment  Treatment (2), Dummy (1, 4, 6), Simple (5)  contr.treatment
 contr_code_treatment

 Anova  Deviation (2), Contrast (1), Simple (4)  contr.treatment  1/k
 contr_code_anova

 Sum  Sum (1, 2, 6), Effects (3), Deviation (4, 5), Unweighted Effects (7)  contr.sum
 contr_code_sum

 Difference  Contrast (3), Forward/Backward (4), Repeated (5)  MASS::contr.sdif
 contr_code_difference

 Helmert  Reverse Helmert (1, 4), Difference (5), Contrast (7)  contr.helmert / (column_i+1)
 contr_code_helmert

 Polynomial  Polynomial (5), Orthogonal Polynomial (4), Trend (3)  contr.poly
 contr_code_poly

These functions are under development
First, we'll set up a simple experimental design and analyse it with lm()
. I set empirical = TRUE
to make interpreting the estimates easier.
df < sim_design(between = list(pet = c("cat", "dog", "ferret")), n = c(50), mu = c(2, 4, 9), empirical = TRUE)
Notice that the default contrast is treatment coding, with "cat" as the baseline condition. The estimate for the intercept is the mean value for cats ($2$), while the term petdog
is the mean value for dogs minus cats ($4  2$), and the term petferret
is the mean value for ferrets minus cats ($9  2$).
r ctable(df)
In all the subsequent examples, try putting into words what the estimates for the intercept and terms mean in relation to the mean values for each group in the data.
Treatment coding (also called dummy coding) sets the mean of the reference group as the intercept. It is straightforward to interpret, especially when your factors have an obvious default or baseline value. This is the same type of coding as the default for factors shown above (unless you change your default using options()
), but with clearer term labels.
df$pet < contr_code_treatment(df$pet)
r ctable(df)
Each contrast compares one level with the reference level, which defaults to the first level, but you can set with the base
argument. Now the intercept estimates the mean value for dogs ($4$).
df$pet < contr_code_treatment(df$pet, base = "dog")
r ctable(df)
Set the reference level to the third level ("ferret").
df$pet < contr_code_treatment(df$pet, base = 3)
r ctable(df)
Anova coding is identical to treatment coding, but sets the grand mean as the intercept. Each contrast compares one level with a reference level. This gives us values that are relatively easy to interpret and map onto ANOVA values.
Below is anova coding with the first level ("cat") as the default base. Now the intercept is the grand mean, which is the mean of the three group means ($(2 + 4 + 9)/3$). Notice that this is different from the mean value of y in our dataset (r mean(df$y)
, since the number of pets in each group is unbalanced. The term pet_dogcat
is the mean value for dogs minus cats ($4  2$) and the term pet_ferretcat
is the mean value for ferrets minus cats ($9  2$).
df$pet < contr_code_anova(df$pet)
r ctable(df)
Anova coding with "dog" as the base. How does the interpretation of the terms change?
df$pet < contr_code_anova(df$pet, base = "dog")
r ctable(df)
Anova coding with the third level ("ferret") as the base.
df$pet < contr_code_anova(df$pet, base = 3)
r ctable(df)
Sum coding also sets the grand mean as the intercept. Each contrast compares one level with the grand mean. Therefore, the estimate for pet_catintercept
is the difference between the mean value for cats and the grand mean ($2  5$).
df$pet < contr_code_sum(df$pet)
r ctable(df)
You can't compare all levels with the grand mean, and have to omit one level. This is the last level by default, but you can change it with the omit
argument.
df$pet < contr_code_sum(df$pet, omit = "dog")
r ctable(df)
Omit the first level ("cat").
df$pet < contr_code_sum(df$pet, omit = 1)
r ctable(df)
A slightly different form of contrast coding is difference coding, also called forward, backward, or successive differences coding. It compares each level to the previous one, rather than to a baseline level.
df$pet < contr_code_difference(df$pet)
r ctable(df)
If you want to change which levels are compared, you can reorder the factor levels.
df$pet < contr_code_difference(df$pet, levels = c("ferret", "cat", "dog"))
r ctable(df)
# get levels back to normal df$pet < contr_code_difference(df$pet, levels = c("cat", "dog", "ferret"))
Helmert coding sets the grand mean as the intercept. Each contrast compares one level with the mean of previous levels. This coding is somewhat different than the results from stats::contr.helmert()
to make it easier to interpret the estimates. Thus, pet_ferretcat.dog
is the mean value for ferrets minus the mean value for cats and dogs averaged together ($9(2+4)/2$).
df$pet < contr_code_helmert(df$pet)
r ctable(df)
You can change the comparisons by reordering the levels.
df$pet < contr_code_helmert(df$pet, levels = c("ferret", "dog", "cat"))
r ctable(df)
Polynomial coding is the default for ordered factors in R. We'll set up a new data simulation with five ordered times.
df < sim_design(list(time = 1:5), mu = 1:5 * 0.25 + (1:5  3)^2 * 0.5, sd = 5, long = TRUE)
The function contr_code_poly()
uses contr.poly()
to set up the polynomial contrasts for the linear (^1
), quadratic (^2
), cubic (^3
), and quartic (^4
) components.
df$time < contr_code_poly(df$time)
r ctable(df, "time")
The function add_contrasts()
lets you add contrasts to a column in a data frame and also adds new columns for each contrast (unless add_cols = FALSE
). This is especially helpful if you want to test only a subset of the contrasts.
df < sim_design(list(time = 1:5), mu = 1:5 * 0.25 + (1:5  3)^2 * 0.5, sd = 5, long = TRUE, plot = FALSE) %>% add_contrast("time", "poly")
# test only the linear and quadratic contrasts lm(y ~ `time^1` + `time^2`, df) %>% broom::tidy()
lm(y ~ `time^1` + `time^2`, df) %>% broom::tidy() %>% knitr::kable() %>% kable_styling()
You can set colnames
to change the default column names for the contrasts. This can be useful if you want to add different codings for the same factor or if the default names are too long.
btwn < list(condition = c("control", "experimental")) df < sim_design(between = btwn, n = 1, plot = FALSE) %>% add_contrast("condition", "treatment", colnames = "cond.tr") %>% add_contrast("condition", "anova", colnames = "cond.aov") %>% add_contrast("condition", "difference", colnames = "cond.dif") %>% add_contrast("condition", "sum", colnames = "cond.sum") %>% add_contrast("condition", "helmert", colnames = "cond.hmt") %>% add_contrast("condition", "poly", colnames = "cond.poly")
r kable(df) %>% kable_styling()
However, if a new column has a duplicate name to an existing column, add_contrast()
will automatically add a contrast suffix to the new column.
btwn < list(pet = c("cat", "dog", "ferret")) df < sim_design(between = btwn, n = 1, plot = FALSE) %>% add_contrast("pet", "treatment") %>% add_contrast("pet", "anova") %>% add_contrast("pet", "sum") %>% add_contrast("pet", "difference") %>% add_contrast("pet", "helmert") %>% add_contrast("pet", "poly")
r kable(df, digits = 2) %>% kable_styling()
Let's use simulated data with empirical = TRUE
to explore how to interpret interactions between two 2level factors coded in different ways.
mu < c(0, 4, 6, 10) df < sim_design(between = list(time = c("am", "pm"), pet = c("cat", "dog")), n = c(50, 60, 70, 80), mu = mu, empirical = TRUE)
Yca = mu[[1]] Yda = mu[[2]] Ycp = mu[[3]] Ydp = mu[[4]] Yc. = ((Yca + Ycp)/2) %>% round(2) Yd. = ((Yda + Ydp)/2) %>% round(2) Y.a = ((Yca + Yda)/2) %>% round(2) Y.p = ((Ycp + Ydp)/2) %>% round(2) Y.. = ((Yca + Yda + Ycp + Ydp)/4) %>% round(2)
The table below shows the cell and marginal means. The notation $Y_{..}$ is used to denote a mean for a specific grouping. The .
is used to indicate the mean over all groups of that factor; Y..
is the grand mean. While you'll usually see the subscripts written as numbers to indicate the factor levels, we're using letters here so you don't have to keep referring to the order of factors and levels.
  cat  dog  mean 

 am  $Y_{ca} = r Yca
$  $Y_{da} = r Yda
$  $Y_{.a} = r Y.a
$ 
 pm  $Y_{cp} = r Ycp
$  $Y_{dp} = r Ydp
$  $Y_{.p} = r Y.p
$ 
 mean  $Y_{c.} = r Yc.
$  $Y_{d.} = r Yd.
$  $Y_{..} = r Y..
$ 
 term  interpretation  formula  value 
::::
 intercept  cat when am  $Y_{ca}$  $r Yca
$ 
 pet.dogcat  dog minus cat, when am  $Y_{da}  Y_{ca}$  $r Yda
 r Yca
= r YdaYca
$ 
 time.pmam  pm minus am, for cats  $Y_{cp}  Y_{ca}$  $r Ycp
 r Yca
= r YcpYca
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)(YdaYca)
$ 
df %>% add_contrast("pet", "treatment") %>% add_contrast("time", "treatment") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $r Y..
$ 
 pet.dogcat  mean dog minus mean cat  $Y_{d.}  Y_{c.}$  $r Yd.
 r Yc.
= r Yd.Yc.
$ 
 time.pmam  mean pm minus mean am  $Y_{.p}  Y_{.a}$  $r Y.p
 r Y.a
= r Y.pY.a
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)  (YdaYca)
$ 
df %>% add_contrast("pet", "anova") %>% add_contrast("time", "anova") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $3$ 
 pet.catintercept  mean cat minus grand mean  $Y_{c.}  Y_{..}$  $r Yc.
 r Y..
= r Yc. Y..
$ 
 time.amintercept  mean am minus grand mean  $Y_{.a}  Y_{..}$  $r Y.a
 r Y..
= r Y.a  Y..
$ 
 pet.catintercept:time.amintercept  cat minus mean when am, minus cat minus mean when pm, divided by 2  $\frac{(Y_{ca}  Y_{.a})  (Y_{cp}  Y_{.p})}{2}$  $\frac{(r Yca
 r Y.a
)  (r Ycp
 r Y.p
)}{2} = r ((YcaY.a)(YcpY.p))/2
$ 
df %>% add_contrast("pet", "sum") %>% add_contrast("time", "sum") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $r Y..
$ 
 pet.dogcat  mean dog minus mean cat  $Y_{d.}  Y_{c.}$  $r Yd.
 r Yc.
= r Yd.Yc.
$ 
 time.pmam  mean pm minus mean am  $Y_{.p}  Y_{.a}$  $r Y.p
 r Y.a
= r Y.pY.a
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)  (YdaYca)
$ 
df %>% add_contrast("pet", "difference") %>% add_contrast("time", "difference") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $r Y..
$ 
 pet.dogcat  mean dog minus mean cat  $Y_{d.}  Y_{c.}$  $r Yd.
 r Yc.
= r Yd.Yc.
$ 
 time.pmam  mean pm minus mean am  $Y_{.p}  Y_{.a}$  $r Y.p
 r Y.a
= r Y.pY.a
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)  (YdaYca)
$ 
df %>% add_contrast("pet", "helmert") %>% add_contrast("time", "helmert") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
N.B. In the case of 2level factors, anova, difference, and Helmert coding are identical. Treatment coding differs only in the intercept.
Remember that interactions can always be described in two ways, since (A1  A2)  (B1  B2) == (A1  B1)  (A2  B2)
. Therefore, "dog minus cat when pm, minus dog minus cat when am" is the same as "pm minus am for dogs, minus pm minus am for cats". The way you describe it in a paper depends on which version maps onto your hypothesis more straightforwardly. The examples above might be written as "the difference between dogs and cats was bigger in the evening than the morning" or "the difference between evening and morning was bigger for dogs than for cats". Make sure you check the plots to make sure you are describing the relationships in the right direction.
Let's use simulated data with empirical = TRUE
to explore how to interpret interactions between a 2level factor and a 3level factor coded in different ways.
mu < c(0, 5, 7, 6, 2, 1) df < sim_design(between = list(time = c("am", "pm"), pet = c("cat", "dog", "ferret")), n = c(50, 60, 70, 80, 90, 100), mu = mu, empirical = TRUE)
Yca = mu[[1]] Yda = mu[[2]] Yfa = mu[[3]] Ycp = mu[[4]] Ydp = mu[[5]] Yfp = mu[[6]] Yc. = ((Yca + Ycp)/2) %>% round(2) Yd. = ((Yda + Ydp)/2) %>% round(2) Yf. = ((Yfa + Yfp)/2) %>% round(2) Y.a = ((Yca + Yda + Yfa)/3) %>% round(2) Y.p = ((Ycp + Ydp + Yfp)/3) %>% round(2) Y.. = ((Yca + Yda + Yfa + Ycp + Ydp + Yfp)/6) %>% round(2)
  cat  dog  ferret  mean 

 am  $Y_{ca} = r Yca
$  $Y_{da} = r Yda
$  $Y_{fa} = r Yfa
$  $Y_{.a} = r Y.a
$ 
 pm  $Y_{cp} = r Ycp
$  $Y_{dp} = r Ydp
$  $Y_{fp} = r Yfp
$  $Y_{.p} = r Y.p
$ 
 mean  $Y_{c.} = r Yc.
$  $Y_{d.} = r Yd.
$  $Y_{f.} = r Yf.
$  $Y_{..} = r Y..
$ 
 term  interpretation  formula  value 
::::
 intercept  cat when am  $Y_{ca}$  $r Yca
$ 
 pet.dogcat  dog minus cat, when am  $Y_{da}  Y_{ca}$  $r Yda
 r Yca
= r YdaYca
$ 
 pet.ferretcat  ferret minus cat, when am  $Y_{fa}  Y_{ca}$  $r Yfa
 r Yca
= r YfaYca
$ 
 time.pmam  pm minus am, for cats  $Y_{cp}  Y_{ca}$  $r Ycp
 r Yca
= r YcpYca
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)(YdaYca)
$ 
 pet.ferretcat:time.pmam  ferret minus cat when pm, minus ferret minus cat when am  $(Y_{fp}  Y_{cp})  (Y_{fa}  Y_{ca})$  $(r Yfp
 r Ycp
)  (r Yfa
 r Yca
) = r (YfpYcp)(YfaYca)
$ 
df %>% add_contrast("pet", "treatment") %>% add_contrast("time", "treatment") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $r Y..
$ 
 pet.dogcat  mean dog minus mean cat  $Y_{d.}  Y_{c.}$  $r Yd.
 r Yc.
= r Yd.Yc.
$ 
 pet.ferretcat  mean ferret minus mean cat  $Y_{f.}  Y_{c.}$  $r Yf.
 r Yc.
= r Yf.Yc.
$ 
 time.pmam  mean pm minus mean am  $Y_{.p}  Y_{.a}$  $r Y.p
 r Y.a
= r Y.pY.a
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)  (YdaYca)
$ 
 pet.ferretcat:time.pmam  ferret minus cat when pm, minus ferret minus cat when am  $(Y_{fp}  Y_{cp})  (Y_{fa}  Y_{ca})$  $(r Yfp
 r Ycp
)  (r Yfa
 r Yca
) = r (YfpYcp)  (YfaYca)
$ 
df %>% add_contrast("pet", "anova") %>% add_contrast("time", "anova") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $3$ 
 pet.catintercept  mean cat minus grand mean  $Y_{c.}  Y_{..}$  $r Yc.
 r Y..
= r Yc. Y..
$ 
 pet.dogintercept  mean dog minus grand mean  $Y_{d.}  Y_{..}$  $r Yd.
 r Y..
= r Yd.  Y..
$ 
 time.amintercept  mean am minus grand mean  $Y_{.a}  Y_{..}$  $r Y.a
 r Y..
= r Y.a  Y..
$ 
 pet.catintercept:time.amintercept  cat minus mean when am, minus cat minus mean when pm, divided by 2  $\frac{(Y_{ca}  Y_{.a})  (Y_{cp}  Y_{.p})}{2}$  $\frac{(r Yca
 r Y.a
)  (r Ycp
 r Y.p
)}{2} = r ((YcaY.a)(YcpY.p))/2
$ 
 pet.dogintercept:time.amintercept  dog minus mean when am, minus dog minus mean when pm, divided by 2  $\frac{(Y_{da}  Y_{.a})  (Y_{dp}  Y_{.p})}{2}$  $\frac{(r Yda
 r Y.a
)  (r Ydp
 r Y.p
)}{2} = r ((YdaY.a)(YdpY.p))/2
$ 
df %>% add_contrast("pet", "sum") %>% add_contrast("time", "sum") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $r Y..
$ 
 pet.dogcat  mean dog minus mean cat  $Y_{d.}  Y_{c.}$  $r Yd.
 r Yc.
= r Yd.Yc.
$ 
 pet.ferretdog  mean ferret minus mean dog  $Y_{f.}  Y_{d.}$  $r Yf.
 r Yd.
= r Yf.Yd.
$ 
 time.pmam  mean pm minus mean am  $Y_{.p}  Y_{.a}$  $r Y.p
 r Y.a
= r Y.pY.a
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)  (YdaYca)
$ 
 pet.ferretdog:time.pmam  ferret minus dog when pm, minus ferret minus dog when am  $(Y_{fp}  Y_{dp})  (Y_{fa}  Y_{da})$  $(r Yfp
 r Ydp
)  (r Yfa
 r Yda
) = r (YfpYdp)  (YfaYda)
$ 
df %>% add_contrast("pet", "difference") %>% add_contrast("time", "difference") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
 term  interpretation  formula  value 
::::
 intercept  grand mean  $Y_{..}$  $r Y..
$ 
 pet.dogcat  mean dog minus mean cat  $Y_{d.}  Y_{c.}$  $r Yd.
 r Yc.
= r Yd.Yc.
$ 
 pet.ferretcat.dog  mean ferret minus mean of cat and dog  $Y_{f.}  \frac{Y_{c.} + Y_{d.}}{2}$  $r Yf.
 \frac{r Yc.
+ r Yd.
}{2} = r Yf.(Yc.+Yd.)/2
$ 
 time.pmam  mean pm minus mean am  $Y_{.p}  Y_{.a}$  $r Y.p
 r Y.a
= r Y.pY.a
$ 
 pet.dogcat:time.pmam  dog minus cat when pm, minus dog minus cat when am  $(Y_{dp}  Y_{cp})  (Y_{da}  Y_{ca})$  $(r Ydp
 r Ycp
)  (r Yda
 r Yca
) = r (YdpYcp)  (YdaYca)
$ 
 pet.ferretcat.dog:time.pmam  ferret minus mean of cat and dog when pm, minus ferret minus mean of cat and dog when am  $(Y_{fp}  \frac{Y_{cp} + Y_{dp}}{2})  (Y_{fa}  \frac{Y_{ca} + Y_{da}}{2})$  $(r Yfp
 \frac{r Ycp
+ r Ydp
}{2})  (r Yfa
 \frac{r Yca
+ r Yda
}{2}) = r (Yfp(Ycp+Ydp)/2)  (Yfa(Yca+Yda)/2)
$ 
df %>% add_contrast("pet", "helmert") %>% add_contrast("time", "helmert") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
N.B. In this case, difference coding is identical to anova coding except that the second pet contrast is ferret versus dog instead of ferret versus cat.
Please just don't. You're going to break my head.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.