1 Preparation
2 Hypotheses
3 Hypothesis 1: simple effects
- 3.1 Custom contrast matrix
- 3.2 With R built-in formula
4 Hypothesis 2: nested interaction
- 4.1 Custom contrast matrix
- 4.2 With R built-in formula
5 Hypothesis 3: three-way interation
6 Summary

In the previous post, I applied priori contrasts in linear models with hypothesis matrix and library(hypr) to test the simple effects, where a 2 \(\times\) 2 (between-subjects) design was used. It is relatively straightforward to understand the 2 \(\times\) 2 experimental design. By contrast, things get more complicated when more independent variables (IV) are included in the models. In this post, I will try to apply the prior contrasts for different hypotheses in linear models with three IVs.

1 Preparation

library(tidyverse)
library(MASS)
library(hypr)

Simulate a dataset with three independent variables (A, B, and C) and each of them has two levels (i.e., a1 vs. a2, b1 vs. b2, and c1 vs. c2).

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

# true values of 8 conditions (used for simulating data)
a1_b1_c1 <- 40
a2_b1_c1 <- 15
a1_b2_c1 <- 30
a2_b2_c1 <- 30
a1_b1_c2 <- 25
a2_b1_c2 <- 15
a1_b2_c2 <- 20
a2_b2_c2 <- 15

# sd for the residuals
e_sd <- 10

set.seed(20201113)
A <- rep(gl(n = length(A_labels), k = N_subj, labels = A_labels), times = 4)
B <- rep(gl(n = length(B_labels), k = N_subj*2, labels = B_labels), times = 2)
C <- gl(n = length(C_labels), k = N_subj*4, labels = C_labels)
y_ <- rep(c(a1_b1_c1, a2_b1_c1, a1_b2_c1, a2_b2_c1,
            a1_b1_c2, a2_b1_c2, a1_b2_c2, a2_b2_c2), each = N_subj)
Y <- y_ + rnorm(length(y_), 0, e_sd)

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

## # A tibble: 400 x 4
##    A     B     C         Y
##    <fct> <fct> <fct> <dbl>
##  1 a1    b1    c1     23.9
##  2 a1    b1    c1     36.1
##  3 a1    b1    c1     36.0
##  4 a1    b1    c1     45.6
##  5 a1    b1    c1     25.1
##  6 a1    b1    c1     36.0
##  7 a1    b1    c1     24.1
##  8 a1    b1    c1     41.7
##  9 a1    b1    c1     36.0
## 10 a1    b1    c1     41.4
## # … with 390 more rows

Create a new variable by combining A, B, and C:

level_names <- paste(rep(A_labels, length.out = 8),
                     rep(B_labels, length.out = 8, each = 2), 
                     rep(C_labels, each = 4),
                     sep = "_")
# further create a new variable by combining A, B, and C
df <- df %>% 
  mutate(I_ = paste(A, B, C, sep = "_"),
         I_ = factor(I_, levels = level_names))
df

## # A tibble: 400 x 5
##    A     B     C         Y I_      
##    <fct> <fct> <fct> <dbl> <fct>   
##  1 a1    b1    c1     23.9 a1_b1_c1
##  2 a1    b1    c1     36.1 a1_b1_c1
##  3 a1    b1    c1     36.0 a1_b1_c1
##  4 a1    b1    c1     45.6 a1_b1_c1
##  5 a1    b1    c1     25.1 a1_b1_c1
##  6 a1    b1    c1     36.0 a1_b1_c1
##  7 a1    b1    c1     24.1 a1_b1_c1
##  8 a1    b1    c1     41.7 a1_b1_c1
##  9 a1    b1    c1     36.0 a1_b1_c1
## 10 a1    b1    c1     41.4 a1_b1_c1
## # … with 390 more rows

2 Hypotheses

With this dataset, suppose we would like to test the following hypotheses:
H1: simple effects (or nested effects) of A nested B and C, e.g., a1_b1_c1 - a2_b1_c1, a1_b2_c1 - a2_b2_c1, etc.

H2: the interaction between A and B for the two levels of C separately, i.e., (a1_b1_c1 - a2_b1_c1) - (a1_b2_c1 - a2_b2_c1) and (a1_b1_c2 - a2_b1_c2) - (a1_b2_c2 - a2_b2_c2).

H3: the three-way interaction among A, B and C, i.e., [(a1_b1_c1 - a2_b1_c1) - (a1_b2_c1 - a2_b2_c1)] - [(a1_b1_c2 - a2_b1_c2) - (a1_b2_c2 - a2_b2_c2)].

3 Hypothesis 1: simple effects

3.1 Custom contrast matrix

The first hypothesis is similar to that in the previous post, where mainly the simple effects were tested.
First, let’s create the corresponding hypothesis matrix:

# "a1_b1_c1" "a2_b1_c1" "a1_b2_c1" "a2_b2_c1" "a1_b1_c2" "a2_b1_c2" "a1_b2_c2" "a2_b2_c2"
hypo_matrix1 <- matrix(c(1, -1,  0,  0,  0,  0,  0,  0,  # H1
                         0,  0,  1, -1,  0,  0,  0,  0,  # H1
                         0,  0,  0,  0,  1, -1,  0,  0,  # H1
                         0,  0,  0,  0,  0,  0,  1, -1,  # H1
                         1/4,  1/4, -1/4, -1/4,  1/4,  1/4, -1/4, -1/4,  # main effect of B
                         1/4,  1/4,  1/4,  1/4, -1/4, -1/4, -1/4, -1/4,  # main effect of C 
                         1/2,  1/2, -1/2, -1/2, -1/2, -1/2,  1/2,  1/2), # interaction between B and C
                       nrow = 8)

Convert the hypothesis matrix to contrast matrix:

con_matrix1 <- fractions(t(ginv(hypo_matrix1)))
rownames(con_matrix1) <- levels(df$I_)
colnames(con_matrix1) <- c("A_b1_c1", "A_b2_c1", "A_b1_c2", "A_b2_c2", "B", "C", "BxC")
# these are orthogonal contrasts (off-diagonal correlations are 0)
con_matrix1

##          A_b1_c1 A_b2_c1 A_b1_c2 A_b2_c2 B    C    BxC 
## a1_b1_c1  1/2       0       0       0     1/2  1/2  1/4
## a2_b1_c1 -1/2       0       0       0     1/2  1/2  1/4
## a1_b2_c1    0     1/2       0       0    -1/2  1/2 -1/4
## a2_b2_c1    0    -1/2       0       0    -1/2  1/2 -1/4
## a1_b1_c2    0       0     1/2       0     1/2 -1/2 -1/4
## a2_b1_c2    0       0    -1/2       0     1/2 -1/2 -1/4
## a1_b2_c2    0       0       0     1/2    -1/2 -1/2  1/4
## a2_b2_c2    0       0       0    -1/2    -1/2 -1/2  1/4

Similarly, we can create the contrast matrix with library(hypr).

# "a1_b1_c1" "a2_b1_c1" "a1_b2_c1" "a2_b2_c1" "a1_b1_c2" "a2_b1_c2" "a1_b2_c2" "a2_b2_c2"
hypr_h1 <- hypr(a1_b1_c1 - a2_b1_c1 ~ 0, # H1
                a1_b2_c1 - a2_b2_c1 ~ 0, # H1
                a1_b1_c2 - a2_b1_c2 ~ 0, # H1
                a1_b2_c2 - a2_b2_c2 ~ 0, # H1
                (a1_b1_c1 + a2_b1_c1 + a1_b1_c2 + a2_b1_c2)/4 - (a1_b2_c1 + a2_b2_c1 + a1_b2_c2 + a2_b2_c2)/4 ~ 0,
                (a1_b1_c1 + a2_b1_c1 + a1_b2_c1 + a2_b2_c1)/4 - (a1_b1_c2 + a2_b1_c2 + a1_b2_c2 + a2_b2_c2)/4 ~ 0,
                ((a1_b1_c1 + a2_b1_c1)/2 - (a1_b2_c1 + a2_b2_c1)/2) - ((a1_b1_c2 + a2_b1_c2)/2 - (a1_b2_c2 + a2_b2_c2)/2) ~ 0,
                levels = level_names)
hypr_h1

## hypr object containing 7 null hypotheses:
## H0.1: 0 = a1_b1_c1 - a2_b1_c1
## H0.2: 0 = a1_b2_c1 - a2_b2_c1
## H0.3: 0 = a1_b1_c2 - a2_b1_c2
## H0.4: 0 = a1_b2_c2 - a2_b2_c2
## H0.5: 0 = 1/4*a1_b1_c1 + 1/4*a2_b1_c1 + 1/4*a1_b1_c2 + 1/4*a2_b1_c2 - 1/4*a1_b2_c1 - 1/4*a2_b2_c1 - 1/4*a1_b2_c2 - 1/4*a2_b2_c2
## H0.6: 0 = 1/4*a1_b1_c1 + 1/4*a2_b1_c1 + 1/4*a1_b2_c1 + 1/4*a2_b2_c1 - 1/4*a1_b1_c2 - 1/4*a2_b1_c2 - 1/4*a1_b2_c2 - 1/4*a2_b2_c2
## H0.7: 0 = 1/2*a1_b1_c1 + 1/2*a2_b1_c1 - 1/2*a1_b2_c1 - 1/2*a2_b2_c1 - 1/2*a1_b1_c2 - 1/2*a2_b1_c2 + 1/2*a1_b2_c2 + 1/2*a2_b2_c2
## 
## Hypothesis matrix (transposed):
##          [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## a1_b1_c1    1    0    0    0  1/4  1/4  1/2
## a2_b1_c1   -1    0    0    0  1/4  1/4  1/2
## a1_b2_c1    0    1    0    0 -1/4  1/4 -1/2
## a2_b2_c1    0   -1    0    0 -1/4  1/4 -1/2
## a1_b1_c2    0    0    1    0  1/4 -1/4 -1/2
## a2_b1_c2    0    0   -1    0  1/4 -1/4 -1/2
## a1_b2_c2    0    0    0    1 -1/4 -1/4  1/2
## a2_b2_c2    0    0    0   -1 -1/4 -1/4  1/2
## 
## Contrast matrix:
##          [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## a1_b1_c1  1/2    0    0    0  1/2  1/2  1/4
## a2_b1_c1 -1/2    0    0    0  1/2  1/2  1/4
## a1_b2_c1    0  1/2    0    0 -1/2  1/2 -1/4
## a2_b2_c1    0 -1/2    0    0 -1/2  1/2 -1/4
## a1_b1_c2    0    0  1/2    0  1/2 -1/2 -1/4
## a2_b1_c2    0    0 -1/2    0  1/2 -1/2 -1/4
## a1_b2_c2    0    0    0  1/2 -1/2 -1/2  1/4
## a2_b2_c2    0    0    0 -1/2 -1/2 -1/2  1/4

The contrast matrix obtained with hypr is the same as what we calculated with generalized matrix inverse from hypothesis matrix. Then we apply this contrast matrix to I_ and fit the model:

# apply the contrast matrix
contrasts(df$I_) <- con_matrix1
# contrasts(df$I_) <- contr.hypothesis(hypr_h1) # or use the contrast from hypr
# fit lm with the custom contrasts
lm_h1 <- lm(Y ~ 1+ I_, data = df)
summary(lm_h1)

## 
## Call:
## lm(formula = Y ~ 1 + I_, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.001  -7.431  -0.280   6.616  37.293 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 22.43184    0.50775  44.179  < 2e-16 ***
## I_A_b1_c1   23.72357    2.03101  11.681  < 2e-16 ***
## I_A_b2_c1   -0.02668    2.03101  -0.013  0.98953    
## I_A_b1_c2   12.31250    2.03101   6.062 3.16e-09 ***
## I_A_b2_c2    4.87371    2.03101   2.400  0.01688 *  
## I_B         -0.31708    1.01551  -0.312  0.75502    
## I_C          9.23209    1.01551   9.091  < 2e-16 ***
## I_BxC       -5.51518    2.03101  -2.715  0.00691 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.16 on 392 degrees of freedom
## Multiple R-squared:  0.407,  Adjusted R-squared:  0.3964 
## F-statistic: 38.44 on 7 and 392 DF,  p-value: < 2.2e-16

coefficient	calculation	true value	estimate
`(Intercept)`	`(a1_b1_c1 + a2_b1_c1 + a1_b2_c1 + a2_b2_c1 + a1_b1_c2 + a2_b1_c2 + a1_b2_c2 + a2_b2_c2)/8`	23.75	22.43
`I_A_b1_c1`	`a1_b1_c1 - a2_b1_c1`	25	23.72
`I_A_b2_c1`	`a1_b2_c1 - a2_b2_c1`	0	-0.03
`I_A_b1_c2`	`a1_b1_c2 - a2_b1_c2`	10	12.31
`I_A_b2_c2`	`a1_b2_c2 - a2_b2_c2`	5	4.87
`I_B`	`(a1_b1_c1 + a2_b1_c1 + a1_b1_c2 + a2_b1_c2)/4 - (a1_b2_c1 + a2_b2_c1 + a1_b2_c2 + a2_b2_c2)/4`	0	-0.32
`I_C`	`(a1_b1_c1 + a2_b1_c1 + a1_b2_c1 + a2_b2_c1)/4 - (a1_b1_c2 + a2_b1_c2 + a1_b2_c2 + a2_b2_c2)/4`	10	9.23
`I_BxC`	`[(a1_b1_c1 + a2_b1_c1)/2 -` `(a1_b2_c1 + a2_b2_c1)/2] -` `[(a1_b1_c2 + a2_b1_c2)/2 -` `(a1_b2_c2 + a2_b2_c2)/2]`	-5	-5.52

We can see from the above table, the estimated parameters matched the values used for simulating data.

3.2 With R built-in formula

R built-in formula can also be applied to test the same hypothesis:

lm_h1R <- lm(Y ~ 1+ B*C/A, data = df, 
             contrasts = list(A = contr.sdif(2),
                              B = contr.sdif(2),
                              C = contr.sdif(2)))
summary(lm_h1R)

## 
## Call:
## lm(formula = Y ~ 1 + B * C/A, data = df, contrasts = list(A = contr.sdif(2), 
##     B = contr.sdif(2), C = contr.sdif(2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.001  -7.431  -0.280   6.616  37.293 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.43184    0.50775  44.179  < 2e-16 ***
## B2-1           0.31708    1.01551   0.312  0.75502    
## C2-1          -9.23209    1.01551  -9.091  < 2e-16 ***
## B2-1:C2-1     -5.51518    2.03101  -2.715  0.00691 ** 
## Bb1:Cc1:A2-1 -23.72357    2.03101 -11.681  < 2e-16 ***
## Bb2:Cc1:A2-1   0.02668    2.03101   0.013  0.98953    
## Bb1:Cc2:A2-1 -12.31250    2.03101  -6.062 3.16e-09 ***
## Bb2:Cc2:A2-1  -4.87371    2.03101  -2.400  0.01688 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.16 on 392 degrees of freedom
## Multiple R-squared:  0.407,  Adjusted R-squared:  0.3964 
## F-statistic: 38.44 on 7 and 392 DF,  p-value: < 2.2e-16

Not surprisingly, the parameters estimated by lm_h1R are the same as those in lm_h1. They are all similar to the values used for simulating the data.

4 Hypothesis 2: nested interaction

Possibly Hypothesis 2 is not that common in the literature. Usually, only when a three-way interaction is significant, researchers may choose to divide one of the variables, and then for each of the levels, explore the interaction of the remaining two variables. I just want to mention here that I do not think we should do this. Instead, we should test the two (or more) two-way interactions only when we have corresponding hypotheses regardless of the three-way interation results.

4.1 Custom contrast matrix

Let’s create the hypothesis matrix at first:

# "a1_b1_c1" "a2_b1_c1" "a1_b2_c1" "a2_b2_c1" "a1_b1_c2" "a2_b1_c2" "a1_b2_c2" "a2_b2_c2"
hypo_matrix2 <- matrix(c(  1,   -1,   -1,    1,    0,    0,    0,    0,  # H2: c1 AxB
                           0,    0,    0,    0,    1,   -1,   -1,    1,  # H2: c2 AxB
                         1/2, -1/2,  1/2, -1/2,    0,    0,    0,    0,  # c1 main effect of A
                           0,    0,    0,    0,  1/2, -1/2,  1/2, -1/2,  # c2 main effect of A
                         1/2,  1/2, -1/2, -1/2,    0,    0,    0,    0,  # c1 main effect of B
                           0,    0,    0,    0,  1/2,  1/2, -1/2, -1/2,  # c2 main effect of B
                         1/4,  1/4,  1/4,  1/4, -1/4, -1/4, -1/4, -1/4),  # main effect of C 
                       nrow = 8)

Convert the hypothesis matrix to contrast matrix:

con_matrix2 <- fractions(t(ginv(hypo_matrix2)))
rownames(con_matrix2) <- levels(df$I_)
colnames(con_matrix2) <- c("AxB_c1", "AxB_c2", "A_c1", "A_c2", "B_c1", "B_c2", "C")
# these are orthogonal contrasts (off-diagonal correlations are 0)
con_matrix2

##          AxB_c1 AxB_c2 A_c1 A_c2 B_c1 B_c2 C   
## a1_b1_c1  1/4      0    1/2    0  1/2    0  1/2
## a2_b1_c1 -1/4      0   -1/2    0  1/2    0  1/2
## a1_b2_c1 -1/4      0    1/2    0 -1/2    0  1/2
## a2_b2_c1  1/4      0   -1/2    0 -1/2    0  1/2
## a1_b1_c2    0    1/4      0  1/2    0  1/2 -1/2
## a2_b1_c2    0   -1/4      0 -1/2    0  1/2 -1/2
## a1_b2_c2    0   -1/4      0  1/2    0 -1/2 -1/2
## a2_b2_c2    0    1/4      0 -1/2    0 -1/2 -1/2

Same contrast matrix created with library(hypr) is shown below.

# "a1_b1_c1" "a2_b1_c1" "a1_b2_c1" "a2_b2_c1" "a1_b1_c2" "a2_b1_c2" "a1_b2_c2" "a2_b2_c2"
hypr_h2 <- hypr((a1_b1_c1 - a2_b1_c1)-(a1_b2_c1 - a2_b2_c1) ~ 0, # H2: c1 AxB
                (a1_b1_c2 - a2_b1_c2)-(a1_b2_c2 - a2_b2_c2) ~ 0, # H2: c2 AxB
                (a1_b1_c1 + a1_b2_c1)/2-(a2_b1_c1 + a2_b2_c1)/2 ~ 0, # H2: c1 main effect of A
                (a1_b1_c2 + a1_b2_c2)/2-(a2_b1_c2 + a2_b2_c2)/2 ~ 0, # H2: c2 main effect of A
                (a1_b1_c1 + a2_b1_c1)/2-(a1_b2_c1 + a2_b2_c1)/2 ~ 0, # H2: c1 main effect of B
                (a1_b1_c2 + a2_b1_c2)/2-(a1_b2_c2 + a2_b2_c2)/2 ~ 0, # H2: c2 main effect of B
                (a1_b1_c1 + a2_b1_c1 + a1_b2_c1 + a2_b2_c1)/4 - (a1_b1_c2 + a2_b1_c2 + a1_b2_c2 + a2_b2_c2)/4 ~ 0, 
                levels = level_names)
hypr_h2

## hypr object containing 7 null hypotheses:
## H0.1: 0 = a1_b1_c1 - a2_b1_c1 - a1_b2_c1 + a2_b2_c1
## H0.2: 0 = a1_b1_c2 - a2_b1_c2 - a1_b2_c2 + a2_b2_c2
## H0.3: 0 = 1/2*a1_b1_c1 + 1/2*a1_b2_c1 - 1/2*a2_b1_c1 - 1/2*a2_b2_c1
## H0.4: 0 = 1/2*a1_b1_c2 + 1/2*a1_b2_c2 - 1/2*a2_b1_c2 - 1/2*a2_b2_c2
## H0.5: 0 = 1/2*a1_b1_c1 + 1/2*a2_b1_c1 - 1/2*a1_b2_c1 - 1/2*a2_b2_c1
## H0.6: 0 = 1/2*a1_b1_c2 + 1/2*a2_b1_c2 - 1/2*a1_b2_c2 - 1/2*a2_b2_c2
## H0.7: 0 = 1/4*a1_b1_c1 + 1/4*a2_b1_c1 + 1/4*a1_b2_c1 + 1/4*a2_b2_c1 - 1/4*a1_b1_c2 - 1/4*a2_b1_c2 - 1/4*a1_b2_c2 - 1/4*a2_b2_c2
## 
## Hypothesis matrix (transposed):
##          [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## a1_b1_c1    1    0  1/2    0  1/2    0  1/4
## a2_b1_c1   -1    0 -1/2    0  1/2    0  1/4
## a1_b2_c1   -1    0  1/2    0 -1/2    0  1/4
## a2_b2_c1    1    0 -1/2    0 -1/2    0  1/4
## a1_b1_c2    0    1    0  1/2    0  1/2 -1/4
## a2_b1_c2    0   -1    0 -1/2    0  1/2 -1/4
## a1_b2_c2    0   -1    0  1/2    0 -1/2 -1/4
## a2_b2_c2    0    1    0 -1/2    0 -1/2 -1/4
## 
## Contrast matrix:
##          [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## a1_b1_c1  1/4    0  1/2    0  1/2    0  1/2
## a2_b1_c1 -1/4    0 -1/2    0  1/2    0  1/2
## a1_b2_c1 -1/4    0  1/2    0 -1/2    0  1/2
## a2_b2_c1  1/4    0 -1/2    0 -1/2    0  1/2
## a1_b1_c2    0  1/4    0  1/2    0  1/2 -1/2
## a2_b1_c2    0 -1/4    0 -1/2    0  1/2 -1/2
## a1_b2_c2    0 -1/4    0  1/2    0 -1/2 -1/2
## a2_b2_c2    0  1/4    0 -1/2    0 -1/2 -1/2

# apply the contrast matrix
contrasts(df$I_) <- con_matrix2
# contrasts(df$I_) <- contr.hypothesis(hypr_h2) # or use the contrast from hypr
# fit lm with the custom contrasts
lm_h2 <- lm(Y ~ 1+ I_, data = df)
summary(lm_h2)

## 
## Call:
## lm(formula = Y ~ 1 + I_, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.001  -7.431  -0.280   6.616  37.293 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.4318     0.5078  44.179  < 2e-16 ***
## I_AxB_c1     23.7503     2.8723   8.269 2.13e-15 ***
## I_AxB_c2      7.4388     2.8723   2.590  0.00996 ** 
## I_A_c1       11.8484     1.4361   8.250 2.43e-15 ***
## I_A_c2        8.5931     1.4361   5.983 4.93e-09 ***
## I_B_c1       -3.0747     1.4361  -2.141  0.03290 *  
## I_B_c2        2.4405     1.4361   1.699  0.09005 .  
## I_C           9.2321     1.0155   9.091  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.16 on 392 degrees of freedom
## Multiple R-squared:  0.407,  Adjusted R-squared:  0.3964 
## F-statistic: 38.44 on 7 and 392 DF,  p-value: < 2.2e-16

coefficient	calculation	true value	estimate
`(Intercept)`	`(a1_b1_c1 + a2_b1_c1 + a1_b2_c1 + a2_b2_c1 + a1_b1_c2 + a2_b1_c2 + a1_b2_c2 + a2_b2_c2)/8`	23.75	22.43
`I_AxB_c1`	`(a1_b1_c1 - a2_b1_c1) - (a1_b2_c1 - a2_b2_c1)`	25	23.75
`I_AxB_c2`	`(a1_b1_c2 - a2_b1_c2) - (a1_b2_c2 - a2_b2_c2)`	5	7.44
`I_A_c1`	`(a1_b1_c1 + a1_b2_c1)/2 - (a2_b1_c1 + a2_b2_c1)/2`	12.5	11.85
`I_A_c2`	`(a1_b1_c2 + a1_b2_c2)/2 - (a2_b1_c2 + a2_b2_c2)/2`	7.5	8.59
`I_B_c1`	`(a1_b1_c1 + a2_b1_c1)/2 - (a1_b2_c1 + a2_b2_c1)/2`	-2.5	-3.07
`I_B_c2`	`(a1_b1_c2 + a2_b1_c2)/2 - (a1_b2_c2 + a2_b2_c2)/2`	2.5	2.44
`I_C`	`(a1_b1_c1 + a2_b1_c1 + a1_b2_c1 + a2_b2_c1)/4 - (a1_b1_c2 + a2_b1_c2 + a1_b2_c2 + a2_b2_c2)/4`	10	9.23

As always, the estimates are quite close to the values used for simulating data.

4.2 With R built-in formula

lm_h2R <- lm(Y ~ 1+ C/(A*B), data = df, 
             contrasts = list(A = contr.sdif(2),
                              B = contr.sdif(2),
                              C = contr.sdif(2)))
summary(lm_h2R)

## 
## Call:
## lm(formula = Y ~ 1 + C/(A * B), data = df, contrasts = list(A = contr.sdif(2), 
##     B = contr.sdif(2), C = contr.sdif(2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.001  -7.431  -0.280   6.616  37.293 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    22.4318     0.5078  44.179  < 2e-16 ***
## C2-1           -9.2321     1.0155  -9.091  < 2e-16 ***
## Cc1:A2-1      -11.8484     1.4361  -8.250 2.43e-15 ***
## Cc2:A2-1       -8.5931     1.4361  -5.983 4.93e-09 ***
## Cc1:B2-1        3.0747     1.4361   2.141  0.03290 *  
## Cc2:B2-1       -2.4405     1.4361  -1.699  0.09005 .  
## Cc1:A2-1:B2-1  23.7503     2.8723   8.269 2.13e-15 ***
## Cc2:A2-1:B2-1   7.4388     2.8723   2.590  0.00996 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.16 on 392 degrees of freedom
## Multiple R-squared:  0.407,  Adjusted R-squared:  0.3964 
## F-statistic: 38.44 on 7 and 392 DF,  p-value: < 2.2e-16

The parameters estimated by lm_h2R are the same as those in lm_h2. They are similar to the values used for simulating data.

5 Hypothesis 3: three-way interation

Among these three hypohtese, probably we are most familiar with the third hypothesis. Results of H3 can be obtained from ANOVA or linear models with most of the contrast coding provided by R, e.g., contr.treatment(), contr.sum(), contr.sdif(). (Note this is not necesasrily true when there are more than two levels in the independent variables.)

Instead of using the custom contrast matrix, I will fit the model with differnet contrast coding available in R.

5.1 Treatment contrasts

# the default contrast coding is treatment contrasts.
lm_h3_treat <- lm(Y ~ 1+ A*B*C, data = df)
summary(lm_h3_treat)

## 
## Call:
## lm(formula = Y ~ 1 + A * B * C, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.001  -7.431  -0.280   6.616  37.293 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   37.372      1.436  26.023  < 2e-16 ***
## Aa2          -23.724      2.031 -11.681  < 2e-16 ***
## Bb2           -8.800      2.031  -4.333 1.87e-05 ***
## Cc2          -12.180      2.031  -5.997 4.57e-09 ***
## Aa2:Bb2       23.750      2.872   8.269 2.13e-15 ***
## Aa2:Cc2       11.411      2.872   3.973 8.45e-05 ***
## Bb2:Cc2        2.641      2.872   0.919    0.358    
## Aa2:Bb2:Cc2  -16.311      4.062  -4.016 7.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.16 on 392 degrees of freedom
## Multiple R-squared:  0.407,  Adjusted R-squared:  0.3964 
## F-statistic: 38.44 on 7 and 392 DF,  p-value: < 2.2e-16

5.2 Sum contrasts

lm_h3_sum <- lm(Y ~ 1+ A*B*C, data = df, 
             contrasts = list(A = contr.sum(2),
                              B = contr.sum(2),
                              C = contr.sum(2)))
summary(lm_h3_sum)

## 
## Call:
## lm(formula = Y ~ 1 + A * B * C, data = df, contrasts = list(A = contr.sum(2), 
##     B = contr.sum(2), C = contr.sum(2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.001  -7.431  -0.280   6.616  37.293 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.4318     0.5078  44.179  < 2e-16 ***
## A1            5.1104     0.5078  10.065  < 2e-16 ***
## B1           -0.1585     0.5078  -0.312  0.75502    
## C1            4.6160     0.5078   9.091  < 2e-16 ***
## A1:B1         3.8986     0.5078   7.678 1.30e-13 ***
## A1:C1         0.8138     0.5078   1.603  0.10978    
## B1:C1        -1.3788     0.5078  -2.715  0.00691 ** 
## A1:B1:C1      2.0389     0.5078   4.016 7.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.16 on 392 degrees of freedom
## Multiple R-squared:  0.407,  Adjusted R-squared:  0.3964 
## F-statistic: 38.44 on 7 and 392 DF,  p-value: < 2.2e-16

5.3 Successive difference contrast

lm_h3_sdif <- lm(Y ~ 1+ A*B*C, data = df, 
             contrasts = list(A = contr.sdif(2),
                              B = contr.sdif(2),
                              C = contr.sdif(2)))
summary(lm_h3_sdif)

## 
## Call:
## lm(formula = Y ~ 1 + A * B * C, data = df, contrasts = list(A = contr.sdif(2), 
##     B = contr.sdif(2), C = contr.sdif(2)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.001  -7.431  -0.280   6.616  37.293 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     22.4318     0.5078  44.179  < 2e-16 ***
## A2-1           -10.2208     1.0155 -10.065  < 2e-16 ***
## B2-1             0.3171     1.0155   0.312  0.75502    
## C2-1            -9.2321     1.0155  -9.091  < 2e-16 ***
## A2-1:B2-1       15.5945     2.0310   7.678 1.30e-13 ***
## A2-1:C2-1        3.2553     2.0310   1.603  0.10978    
## B2-1:C2-1       -5.5152     2.0310  -2.715  0.00691 ** 
## A2-1:B2-1:C2-1 -16.3115     4.0620  -4.016 7.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.16 on 392 degrees of freedom
## Multiple R-squared:  0.407,  Adjusted R-squared:  0.3964 
## F-statistic: 38.44 on 7 and 392 DF,  p-value: < 2.2e-16

5.4 Interim summary

Let’s compare the interaction results of these three models with the true value:

contrasts	t value	estimate	true value
`treatment`	-4.02	-16.31
`sum`	4.02	2.04 (2.04\(\times\)`8` = 16.32)	-20
`sdif`	-4.02	-16.31

The interaction is calcuated as [(a2_b2_c2 - a1_b2_c2) - (a2_b1_c2 - a1_b1_c2)] - [(a2_b2_c1 - a1_b2_c1) - (a2_b1_c1 - a1_b1_c1)].

The interaction results of these models fitted with different contrast coding are similar, which are close to the values used for simulating data.

5.5 What happened to the interaction estimate in `lm_h3_sum`?

It is noticeable that the interaction estimate in lm_h3_sum does not really match the true value. We may convert the corresponding design matrix to hypothesis (via contrast matrix) to see what it is testing?

The following codes are “copied” from Schad et al (2020, P29)..
From design matrix to contrast matrix:

# apply sum contrast to these variables
contrasts(df$A) <- contr.sum(2)
contrasts(df$B) <- contr.sum(2)
contrasts(df$C) <- contr.sum(2)

contr_matrix <- df %>% 
  group_by(A, B, C) %>% 
  summarize() %>% 
  model.matrix(~ A*B*C, .) %>% 
  as.data.frame() %>% 
  as.matrix()

## `summarise()` has grouped output by 'A', 'B'. You can override using the `.groups` argument.

contr_matrix

##   (Intercept) A1 B1 C1 A1:B1 A1:C1 B1:C1 A1:B1:C1
## 1           1  1  1  1     1     1     1        1
## 2           1  1  1 -1     1    -1    -1       -1
## 3           1  1 -1  1    -1     1    -1       -1
## 4           1  1 -1 -1    -1    -1     1        1
## 5           1 -1  1  1    -1    -1     1       -1
## 6           1 -1  1 -1    -1     1    -1        1
## 7           1 -1 -1  1     1    -1    -1        1
## 8           1 -1 -1 -1     1     1     1       -1

From conteast matrix to hypothesis matrix:

(hypo_matrix <- fractions(t(ginv(contr_matrix))))

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]  1/8  1/8  1/8  1/8  1/8  1/8  1/8  1/8
## [2,]  1/8  1/8  1/8 -1/8  1/8 -1/8 -1/8 -1/8
## [3,]  1/8  1/8 -1/8  1/8 -1/8  1/8 -1/8 -1/8
## [4,]  1/8  1/8 -1/8 -1/8 -1/8 -1/8  1/8  1/8
## [5,]  1/8 -1/8  1/8  1/8 -1/8 -1/8  1/8 -1/8
## [6,]  1/8 -1/8  1/8 -1/8 -1/8  1/8 -1/8  1/8
## [7,]  1/8 -1/8 -1/8  1/8  1/8 -1/8 -1/8  1/8
## [8,]  1/8 -1/8 -1/8 -1/8  1/8  1/8  1/8 -1/8

These eight columns correspond to the hypotheses for the eight parameters estimated in lm_h3_sum. The last column in hypo_matrix is the hypothesis for the last contrast, which is testing one eighth of the interaction among A, B, and C. That’s why I multiplied it by 8 when comparing to the true value.

6 Summary

In this post, I applied custom hypothesis matrix to test three kinds of popular hypotheses in a \(2 \times2\times2\) (between-subjects) design. I think now I have a better understanding of applying priori contrasts in linear models (possibly also linear mixed models). Thanks for reading this far; hope this post is helpful for you as well.

References

Schad, D. J., Vasishth, S., Hohenstein, S., & Kliegl, R. (2020). How to capitalize on a priori contrasts in linear (mixed) models: A tutorial. Journal of Memory and Language, 110, 104038. https://doi.org/10.1016/j.jml.2019.104038

How to apply priori contrasts in linear models (with three IVs)?