Nothing
Contrasts
In faux: Simulation for Factorial Designs
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 three-level categories (or worse).
I thought it was just me, but some recent Twitter discussion showed me I'm not alone. There's a serious jingle-jangle problem with contrasts. Here are some of the explainers that I found useful (thanks to everyone who recommended them).
- Contrasts in R by Marissa Barlaz
- Coding categorical predictor variables in factorial designs by Dale Barr
- Contrast Coding in Multiple Regression Analysis by Matthew Davis
- R Library Contrast Coding Systems for categorical variables by UCLA Statistical Consulting
- Rosetta store: Contrasts by GAMLj
- Coding Schemes for Categorical Variables by Phillip Alday
- Experimental personality designs: analyzing categorical by continuous variable interactions by West, Aiken & Krull
Terminology
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 |
Faux Contrast Functions
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.
{.tabset -}
Treatment {#treatment -}
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 {#anova -}
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_dog-cat is the mean value for dogs minus cats ($4 - 2$) and the term pet_ferret-cat 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 {#sum -}
Sum coding also sets the grand mean as the intercept. Each contrast compares one level with the grand mean. Therefore, the estimate for pet_cat-intercept 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)
Difference {#difference -}
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 re-order 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 {#helmert -}
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_ferret-cat.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 {#poly -}
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")
Add Contrasts
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()
Examples
2x2 Design
Let's use simulated data with empirical = TRUE to explore how to interpret interactions between two 2-level 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..$ |
{.tabset}
Treatment {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | cat when am | $Y_{ca}$ | $r Yca$ |
| pet.dog-cat | dog minus cat, when am | $Y_{da} - Y_{ca}$ | $r Yda - r Yca = r Yda-Yca$ |
| time.pm-am | pm minus am, for cats | $Y_{cp} - Y_{ca}$ | $r Ycp - r Yca = r Ycp-Yca$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp)-(Yda-Yca)$ |
df %>%
add_contrast("pet", "treatment") %>%
add_contrast("time", "treatment") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Anova {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
df %>%
add_contrast("pet", "anova") %>%
add_contrast("time", "anova") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Sum {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $3$ |
| pet.cat-intercept | mean cat minus grand mean | $Y_{c.} - Y_{..}$ | $r Yc. - r Y.. = r Yc. -Y..$ |
| time.am-intercept | mean am minus grand mean | $Y_{.a} - Y_{..}$ | $r Y.a - r Y.. = r Y.a - Y..$ |
| pet.cat-intercept:time.am-intercept | 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 ((Yca-Y.a)-(Ycp-Y.p))/2$ |
df %>%
add_contrast("pet", "sum") %>%
add_contrast("time", "sum") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Difference {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
df %>%
add_contrast("pet", "difference") %>%
add_contrast("time", "difference") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Helmert {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
df %>%
add_contrast("pet", "helmert") %>%
add_contrast("time", "helmert") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
{-}
N.B. In the case of 2-level 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.
2x3 Design
Let's use simulated data with empirical = TRUE to explore how to interpret interactions between a 2-level factor and a 3-level 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..$ |
{.tabset -}
Treatment {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | cat when am | $Y_{ca}$ | $r Yca$ |
| pet.dog-cat | dog minus cat, when am | $Y_{da} - Y_{ca}$ | $r Yda - r Yca = r Yda-Yca$ |
| pet.ferret-cat | ferret minus cat, when am | $Y_{fa} - Y_{ca}$ | $r Yfa - r Yca = r Yfa-Yca$ |
| time.pm-am | pm minus am, for cats | $Y_{cp} - Y_{ca}$ | $r Ycp - r Yca = r Ycp-Yca$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp)-(Yda-Yca)$ |
| pet.ferret-cat:time.pm-am | 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 (Yfp-Ycp)-(Yfa-Yca)$ |
df %>%
add_contrast("pet", "treatment") %>%
add_contrast("time", "treatment") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Anova {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| pet.ferret-cat | mean ferret minus mean cat | $Y_{f.} - Y_{c.}$ | $r Yf. - r Yc. = r Yf.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
| pet.ferret-cat:time.pm-am | 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 (Yfp-Ycp) - (Yfa-Yca)$ |
df %>%
add_contrast("pet", "anova") %>%
add_contrast("time", "anova") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Sum {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $3$ |
| pet.cat-intercept | mean cat minus grand mean | $Y_{c.} - Y_{..}$ | $r Yc. - r Y.. = r Yc. -Y..$ |
| pet.dog-intercept | mean dog minus grand mean | $Y_{d.} - Y_{..}$ | $r Yd. - r Y.. = r Yd. - Y..$ |
| time.am-intercept | mean am minus grand mean | $Y_{.a} - Y_{..}$ | $r Y.a - r Y.. = r Y.a - Y..$ |
| pet.cat-intercept:time.am-intercept | 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 ((Yca-Y.a)-(Ycp-Y.p))/2$ |
| pet.dog-intercept:time.am-intercept | 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 ((Yda-Y.a)-(Ydp-Y.p))/2$ |
df %>%
add_contrast("pet", "sum") %>%
add_contrast("time", "sum") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Difference {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| pet.ferret-dog | mean ferret minus mean dog | $Y_{f.} - Y_{d.}$ | $r Yf. - r Yd. = r Yf.-Yd.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
| pet.ferret-dog:time.pm-am | 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 (Yfp-Ydp) - (Yfa-Yda)$ |
df %>%
add_contrast("pet", "difference") %>%
add_contrast("time", "difference") %>%
lm(y ~ pet * time, .) %>%
broom::tidy() %>% kable() %>% kable_styling()
Helmert {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| pet.ferret-cat.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.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
| pet.ferret-cat.dog:time.pm-am | 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.
3x3 Design
Please just don't. You're going to break my head.
Try the faux package in your browser
Any scripts or data that you put into this service are public.
faux documentation built on April 20, 2023, 9:13 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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 three-level categories (or worse).
I thought it was just me, but some recent Twitter discussion showed me I'm not alone. There's a serious jingle-jangle problem with contrasts. Here are some of the explainers that I found useful (thanks to everyone who recommended them).
- Contrasts in R by Marissa Barlaz
- Coding categorical predictor variables in factorial designs by Dale Barr
- Contrast Coding in Multiple Regression Analysis by Matthew Davis
- R Library Contrast Coding Systems for categorical variables by UCLA Statistical Consulting
- Rosetta store: Contrasts by GAMLj
- Coding Schemes for Categorical Variables by Phillip Alday
- Experimental personality designs: analyzing categorical by continuous variable interactions by West, Aiken & Krull
Terminology
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 |
Faux Contrast Functions
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.
{.tabset -}
Treatment {#treatment -}
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 {#anova -}
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_dog-cat is the mean value for dogs minus cats ($4 - 2$) and the term pet_ferret-cat 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 {#sum -}
Sum coding also sets the grand mean as the intercept. Each contrast compares one level with the grand mean. Therefore, the estimate for pet_cat-intercept 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)
Difference {#difference -}
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 re-order 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 {#helmert -}
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_ferret-cat.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 {#poly -}
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")
Add Contrasts
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()
Examples
2x2 Design
Let's use simulated data with empirical = TRUE to explore how to interpret interactions between two 2-level 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..$ |
{.tabset}
Treatment {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | cat when am | $Y_{ca}$ | $r Yca$ |
| pet.dog-cat | dog minus cat, when am | $Y_{da} - Y_{ca}$ | $r Yda - r Yca = r Yda-Yca$ |
| time.pm-am | pm minus am, for cats | $Y_{cp} - Y_{ca}$ | $r Ycp - r Yca = r Ycp-Yca$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp)-(Yda-Yca)$ |
df %>% add_contrast("pet", "treatment") %>% add_contrast("time", "treatment") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Anova {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
df %>% add_contrast("pet", "anova") %>% add_contrast("time", "anova") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Sum {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $3$ |
| pet.cat-intercept | mean cat minus grand mean | $Y_{c.} - Y_{..}$ | $r Yc. - r Y.. = r Yc. -Y..$ |
| time.am-intercept | mean am minus grand mean | $Y_{.a} - Y_{..}$ | $r Y.a - r Y.. = r Y.a - Y..$ |
| pet.cat-intercept:time.am-intercept | 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 ((Yca-Y.a)-(Ycp-Y.p))/2$ |
df %>% add_contrast("pet", "sum") %>% add_contrast("time", "sum") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Difference {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
df %>% add_contrast("pet", "difference") %>% add_contrast("time", "difference") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Helmert {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
df %>% add_contrast("pet", "helmert") %>% add_contrast("time", "helmert") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
{-}
N.B. In the case of 2-level 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.
2x3 Design
Let's use simulated data with empirical = TRUE to explore how to interpret interactions between a 2-level factor and a 3-level 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..$ |
{.tabset -}
Treatment {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | cat when am | $Y_{ca}$ | $r Yca$ |
| pet.dog-cat | dog minus cat, when am | $Y_{da} - Y_{ca}$ | $r Yda - r Yca = r Yda-Yca$ |
| pet.ferret-cat | ferret minus cat, when am | $Y_{fa} - Y_{ca}$ | $r Yfa - r Yca = r Yfa-Yca$ |
| time.pm-am | pm minus am, for cats | $Y_{cp} - Y_{ca}$ | $r Ycp - r Yca = r Ycp-Yca$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp)-(Yda-Yca)$ |
| pet.ferret-cat:time.pm-am | 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 (Yfp-Ycp)-(Yfa-Yca)$ |
df %>% add_contrast("pet", "treatment") %>% add_contrast("time", "treatment") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Anova {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| pet.ferret-cat | mean ferret minus mean cat | $Y_{f.} - Y_{c.}$ | $r Yf. - r Yc. = r Yf.-Yc.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
| pet.ferret-cat:time.pm-am | 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 (Yfp-Ycp) - (Yfa-Yca)$ |
df %>% add_contrast("pet", "anova") %>% add_contrast("time", "anova") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Sum {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $3$ |
| pet.cat-intercept | mean cat minus grand mean | $Y_{c.} - Y_{..}$ | $r Yc. - r Y.. = r Yc. -Y..$ |
| pet.dog-intercept | mean dog minus grand mean | $Y_{d.} - Y_{..}$ | $r Yd. - r Y.. = r Yd. - Y..$ |
| time.am-intercept | mean am minus grand mean | $Y_{.a} - Y_{..}$ | $r Y.a - r Y.. = r Y.a - Y..$ |
| pet.cat-intercept:time.am-intercept | 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 ((Yca-Y.a)-(Ycp-Y.p))/2$ |
| pet.dog-intercept:time.am-intercept | 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 ((Yda-Y.a)-(Ydp-Y.p))/2$ |
df %>% add_contrast("pet", "sum") %>% add_contrast("time", "sum") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Difference {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| pet.ferret-dog | mean ferret minus mean dog | $Y_{f.} - Y_{d.}$ | $r Yf. - r Yd. = r Yf.-Yd.$ |
| time.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
| pet.ferret-dog:time.pm-am | 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 (Yfp-Ydp) - (Yfa-Yda)$ |
df %>% add_contrast("pet", "difference") %>% add_contrast("time", "difference") %>% lm(y ~ pet * time, .) %>% broom::tidy() %>% kable() %>% kable_styling()
Helmert {-}
| term | interpretation | formula | value |
|:-----|:---------------|:--------|:------|
| intercept | grand mean | $Y_{..}$ | $r Y..$ |
| pet.dog-cat | mean dog minus mean cat | $Y_{d.} - Y_{c.}$ | $r Yd. - r Yc. = r Yd.-Yc.$ |
| pet.ferret-cat.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.pm-am | mean pm minus mean am | $Y_{.p} - Y_{.a}$ | $r Y.p - r Y.a = r Y.p-Y.a$ |
| pet.dog-cat:time.pm-am | 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 (Ydp-Ycp) - (Yda-Yca)$ |
| pet.ferret-cat.dog:time.pm-am | 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.
3x3 Design
Please just don't. You're going to break my head.
Try the faux package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.