1 Preparation

library(tidyverse)
library(MASS)

Simulate a dataset:

# a 2*2 between-subjects design
N_subj <- 50   # each group
A_labels <- c("a1", "a2")
B_labels <- c("b1", "b2")

# true values (used for simulating data)
a1_b1 <- 30
a2_b1 <- 20
a1_b2 <- 25
a2_b2 <- 50

# sd for the residuals
e_sd <- 10

set.seed(20201111)
A <- rep(gl(n = length(A_labels), k = N_subj, labels = A_labels), times = 2)
B <- gl(n = length(B_labels), k = N_subj *2, labels = B_labels)
y_ <- rep(c(a1_b1, a2_b1, a1_b2, a2_b2), each = N_subj)
Y <- y_ + rnorm(length(y_), 0, e_sd)

(df <- tibble(A, B, Y))

## # A tibble: 200 x 3
##    A     B         Y
##    <fct> <fct> <dbl>
##  1 a1    b1    35.9 
##  2 a1    b1    41.7 
##  3 a1    b1    30.6 
##  4 a1    b1    35.6 
##  5 a1    b1    36.2 
##  6 a1    b1     9.82
##  7 a1    b1    29.8 
##  8 a1    b1    21.0 
##  9 a1    b1    41.3 
## 10 a1    b1    32.6 
## # … with 190 more rows

2 Linear models with the default contrast (treatment contrast)

lm_default <- lm(Y ~ A * B, data = df)
summary(lm_default)

## 
## Call:
## lm(formula = Y ~ A * B, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.003  -5.772   1.103   6.386  24.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   31.087      1.386  22.425  < 2e-16 ***
## Aa2          -11.797      1.960  -6.018 8.55e-09 ***
## Bb2           -7.713      1.960  -3.934 0.000116 ***
## Aa2:Bb2       37.535      2.773  13.538  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.802 on 196 degrees of freedom
## Multiple R-squared:  0.5813, Adjusted R-squared:  0.5749 
## F-statistic: 90.72 on 3 and 196 DF,  p-value: < 2.2e-16

With the default contrast, the true values of Intercept and other coefficients are (compared with the fitted values):

From the table above, we can see the fitted parameters of lm_default match the values used for simulating data.

3 Linear model with custom contrasts

With default contrasts, we obtained the parameters estimating the mean of a1_b1, the differences between a2_b1 and a1_b1 etc. What if we would like to test the differences between the two levels of A under b1 and b2 separately (i.e., the simple effects)?

First, we need to create a new independent variable, which is the combinations of A and B, and write down the hypothesis matrix:

# level names
contrast_names <- c("a1_b1", "a2_b1", "a1_b2", "a2_b2")

# create the new independnet variable Custom
df <- df %>% 
  mutate(C_ = paste(A, B, sep="_"),
         C_ = factor(C_, levels = contrast_names),
         CH_ = C_)

hypothesis_matrix <- matrix(c(  -1,    1,    0,    0, # simple effect 1: a2_b1 - a1_b1
                                 0,    0,   -1,    1, # simple effect 2: a2_b2 - a1_b2
                              -1/2, -1/2,  1/2,  1/2), nrow = 4) # the main effects of B*
fractions(hypothesis_matrix)

##      [,1] [,2] [,3]
## [1,]   -1    0 -1/2
## [2,]    1    0 -1/2
## [3,]    0   -1  1/2
## [4,]    0    1  1/2

*Why the third contrast is the main effects of B? The short answer is that then this model can explain as most of the variances as possible. For a long answer, please refer to Schad et al (2020, P26).

Second, we convert this hypothesis matrix to a contrast matrix, and assign it to the newly added independent variable C_:

# apply the generalized matrix inverse 
contrast_matrix <- fractions(t(ginv(hypothesis_matrix)))
rownames(contrast_matrix) <- levels(df$C_)
colnames(contrast_matrix) <- c("a2_b1-a1_b1", "a2_b2-a1_b2", "B2-1")

# apply the contrasts
(contrasts(df$C_) <- contrast_matrix)

##       a2_b1-a1_b1 a2_b2-a1_b2 B2-1
## a1_b1 -1/2           0        -1/2
## a2_b1  1/2           0        -1/2
## a1_b2    0        -1/2         1/2
## a2_b2    0         1/2         1/2

Then we can fit the model with the addedly independent variable as usual:

# fit lm with the custom contrasts
lm_custom <- lm(Y ~ 1+ C_, data = df)
summary(lm_custom)

## 
## Call:
## lm(formula = Y ~ 1 + C_, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.003  -5.772   1.103   6.386  24.202 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    30.7153     0.6931  44.314  < 2e-16 ***
## C_a2_b1-a1_b1 -11.7974     1.9605  -6.018 8.55e-09 ***
## C_a2_b2-a1_b2  25.7376     1.9605  13.128  < 2e-16 ***
## C_B2-1         11.0542     1.3863   7.974 1.25e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.802 on 196 degrees of freedom
## Multiple R-squared:  0.5813, Adjusted R-squared:  0.5749 
## F-statistic: 90.72 on 3 and 196 DF,  p-value: < 2.2e-16

From the table above, we can see the fitted parameters of lm_custom match the values used for simulating data. Specifically, the three “beta” parameters estimate the simple effects of A under B and the main effects of B.

4 Custom contrasts with library(hypr)

To create the above custom contrasts, we do not have to write down the hypothesis matrix and apply the generalized matrix inverse by ourselves. We can use library(hypr) to obtain the contrasts directly from the hypothesis.

library(hypr) is used to create the same contrasts as the previous section. First, we write down the hypothesis (but not the hypothesis matrix) and then we obtain the hypothesis matrix and contrast matrix directly, which are the same as the previous section.

library(hypr)
hypr_custom <- hypr(a2_b1-a1_b1 ~ 0, # simple effect 1: a2_b1 - a1_b1
                    a2_b2-a1_b2 ~ 0, # simple effect 2: a2_b2 - a1_b2
                    (a1_b1+a2_b1)/2-(a1_b2+a2_b2)/2 ~ 0, # the main effects of B
                    levels = contrast_names)
hypr_custom

## hypr object containing 3 null hypotheses:
## H0.1: 0 = a2_b1 - a1_b1
## H0.2: 0 = a2_b2 - a1_b2
## H0.3: 0 = 1/2*a1_b1 + 1/2*a2_b1 - 1/2*a1_b2 - 1/2*a2_b2
## 
## Hypothesis matrix (transposed):
##       [,1] [,2] [,3]
## a1_b1   -1    0  1/2
## a2_b1    1    0  1/2
## a1_b2    0   -1 -1/2
## a2_b2    0    1 -1/2
## 
## Contrast matrix:
##       [,1] [,2] [,3]
## a1_b1 -1/2    0  1/2
## a2_b1  1/2    0  1/2
## a1_b2    0 -1/2 -1/2
## a2_b2    0  1/2 -1/2

Next, we apply the contrast matrix to the newly added independent variable:

# CH_ was copied from C_ in last section
(contrasts(df$CH_) <- contr.hypothesis(hypr_custom))

##       [,1] [,2] [,3]
## a1_b1 -0.5  0.0  0.5
## a2_b1  0.5  0.0  0.5
## a1_b2  0.0 -0.5 -0.5
## a2_b2  0.0  0.5 -0.5
## attr(,"class")
## [1] "hypr_cmat" "matrix"    "array"

Then we fit the model:

lm_hypr <- lm(Y ~ 1+ CH_, data = df)
summary(lm_hypr)

## 
## Call:
## lm(formula = Y ~ 1 + CH_, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.003  -5.772   1.103   6.386  24.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.7153     0.6931  44.314  < 2e-16 ***
## CH_1        -11.7974     1.9605  -6.018 8.55e-09 ***
## CH_2         25.7376     1.9605  13.128  < 2e-16 ***
## CH_3        -11.0542     1.3863  -7.974 1.25e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.802 on 196 degrees of freedom
## Multiple R-squared:  0.5813, Adjusted R-squared:  0.5749 
## F-statistic: 90.72 on 3 and 196 DF,  p-value: < 2.2e-16

As in the last two sections, the fitted values match the values used for simulating data. Someone may notice that there is discrepancy between the results of lm_hypr and lm_custom. The main effects of B in lm_hypr and lm_custom have opposite signs, though their absolute values are similar. This is because in lm_custom, the main effect of B (i.e., the third contrast) is defined as b2 - b1, while in lm_hypr the same effect is defined as b1 - b2.

5 Simple effect analysss with R code (nested effects)

In the previous two sections, we applied custom contrasts to test the simple effects of A on the two levels of B separately. Instead of using custom contrasts, we can use R code in linear models to test these effects direclty:

lm_R <- lm(Y ~ B/A, df, 
           contrasts = list(A = contr.sdif(2),
                            B = contr.sdif(2)))
summary(lm_R)

## 
## Call:
## lm(formula = Y ~ B/A, data = df, contrasts = list(A = contr.sdif(2), 
##     B = contr.sdif(2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.003  -5.772   1.103   6.386  24.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.7153     0.6931  44.314  < 2e-16 ***
## B2-1         11.0542     1.3863   7.974 1.25e-13 ***
## Bb1:A2-1    -11.7974     1.9605  -6.018 8.55e-09 ***
## Bb2:A2-1     25.7376     1.9605  13.128  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.802 on 196 degrees of freedom
## Multiple R-squared:  0.5813, Adjusted R-squared:  0.5749 
## F-statistic: 90.72 on 3 and 196 DF,  p-value: < 2.2e-16

As always, the fitted value match the values used for simulating data.

6 Summary

This post is kind of a super short summary of Schad et al (2020). For anyone who is interested in using linear model (or lmm), I highly recommend reading this paper. Probably applying priori contrasts in linear models (or lmm) is not that straightforward (at least at the beginning), even it may scare some users off. I do think it is a more beneficial approach, which at least (should) generates more accurate estimations.