# Logic behind writing custom planned contrasts for a triple interaction

I have a model with a triple interaction involving two factors and one continuous variable.

Replicable example:

# %%%
set.seed(1234)

make_times <- function(n, n_per_group) {
vector <- c()
for (i in 1:n) {
one_increment_sequence <- seq(0, to = n_per_group - 1)
increment_vecor <- c(0, c(runif(n_per_group - 1, 0, 1)))
vector <- c(vector, one_increment_sequence + increment_vecor)
}
return(vector)
}

# Define the dataset
dataset <- data.frame(
Outcome = rnorm(60),
Group = rep(c("A", "B", "C"), each = 20),
Sex = rep(c("M", "F"), each = 15, times = 4),
TimeFromBaseline = make_times(20, 3), # Continuous
Extra = sample(c("Foo", "Bar"), 60, replace = TRUE) # Fixed Covariate
)
dataset$$Group <- factor(dataset$$Group, levels = c("A", "B", "C", "D"))
dataset$$Sex <- factor(dataset$$Sex, levels = c("M", "F"))
dataset$$Extra <- factor(dataset$$Extra, levels = c("Foo", "Bar"))
str(dataset)

# Fit a model
model <- lm(Outcome ~ Sex * Group * TimeFromBaseline + Extra, data = dataset)
summary(model)

# Plot the model
sjPlot::plot_model(model, type = "pred", terms = c("TimeFromBaseline", "Group", "Sex"))


There are 13 coefficients in the model. My goal is to define length-13 contrast vectors that test certain hypotheses. This is because I will be conducting brain voxelwise analysis and the software requires manually specifying planned contrasts a priori

I have a basic understanding of what each coefficient represents in treatment encoding scheme R uses by default:

• Beta0 is the mean of the reference level (in this case Group A, Male, at the intercept)
• Beta1 is would be the increment from the base reference level, changing Male -> Female
• etc.

So, presumably, the contrast vector testing whether a male from Group A has a non-zero outcome at time zero (intercept) would be:

c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

Or the contrast vector describing the slope of a female from Group B would be

c(1, 1, 1, 0, 0.04, 0, 1, 0, -1.35, -0.36, 0, 1.86, 0)

Right?

If so, where it gets difficult for me is understanding how to construct hypothesis-testing contrast matrices.

For example, what contrast matrix here (presumably with 13 columns) would test the hypothesis: the triple interaction is greater than zero

^ Can anyone give a rundown of how such a contrast matrix would be constructed?

The contrasts provided to packages like emmeans or lsmeans or multcomp are short vectors that specifically target what is provided to their library-focus object like emmGrid. That's not what I'm looking for.

P.S If there is an R library that can provide these "full" or "a priori" contrast matrices from a model like how model.matrix provides a design matrix, that would make my day.

As I think you are demonstrating, trying to interpret coefficients becomes increasingly difficult as the number of factors increases. However, you can get the needed weights easily if you can get emmeans to do the contrasts you want.

Suppose, for example, you want Dunnett-style comparisons of Group means for various combinations of Sex and TimeFromBaseline. To do this in emmeans, you'd code:

> emm <- emmeans(model, ~ Group | Sex * TimeFromBaseline, at = list(TimeFromBaseline = 0:3))
> con <- contrast(emm, "dunnett")


This will create 2 x 4 = 8 sets of 2 contrasts; here are the first few:

> head(con)
contrast Sex TimeFromBaseline estimate    SE  df t.ratio p.value
B - A    M                  0    0.542 0.388 107   1.396  0.5307
C - A    M                  0    1.061 0.646 107   1.641  0.3803
B - A    F                  0   -1.031 0.536 107  -1.924  0.2362
C - A    F                  0   -0.920 0.493 107  -1.865  0.2628
B - A    M                  1   -0.166 0.261 107  -0.637  0.9313
C - A    M                  1    0.329 0.396 107   0.832  0.8589


con is an object of class emmGrid and it has everything it needs to estimate these contrasts - the regression coefficients, covariances, and the linear functions of those coefficients. The latter is what you need, and they are stored in the linfct slot. So to see them, just look at that slot:

> head(con@linfct)
(Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC
[1,]           0    0      1      0                0        0           0           0
[2,]           0    0      0      1                0        0           0           0
[3,]           0    0      1      0                0        0           1           0
[4,]           0    0      0      1                0        0           0           1
[5,]           0    0      1      0                0        0           0           0
[6,]           0    0      0      1                0        0           0           0
SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline
[1,]                     0                       0                       0
[2,]                     0                       0                       0
[3,]                     0                       0                       0
[4,]                     0                       0                       0
[5,]                     0                       1                       0
[6,]                     0                       0                       1
SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline
[1,]                            0                            0
[2,]                            0                            0
[3,]                            0                            0
[4,]                            0                            0
[5,]                            0                            0
[6,]                            0                            0


Each row of this matrix comprises the 13 weights needed to estimate the corresponding contrast.

### For the 3-way interaction itself

For this, we need an interaction contrast. This entails setting up two different values of each covariate (timeFromBaseline in this case); then the "consec" contrasts avoid any linear dependence among the contrasts:

> emm2 = emmeans(model, ~ Group * Sex * TimeFromBaseline, cov.reduce = range)

> (con2 <- contrast(emm2, interaction = "consec"))
Group_consec Sex_consec TimeFromBaseline_consec estimate   SE  df t.ratio p.value
B - A        F - M      2.99874080810696 - 0       3.377 1.18 107   2.866  0.0050
C - B        F - M      2.99874080810696 - 0      -0.278 1.24 107  -0.224  0.8233

Results are averaged over the levels of: Extra

> con2@linfct
(Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC
[1,]           0    0      0      0                0        0           0           0
[2,]           0    0      0      0                0        0           0           0
SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline
[1,]                     0                       0                       0
[2,]                     0                       0                       0
SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline
[1,]                     2.998741                     0.000000
[2,]                    -2.998741                     2.998741


These two contrasts define the 2 d.f. for the 3-way interaction. [They'd have come out as less messy if we had specified at = list(TimeAfterBaseline = c(0,3)) instead of cov.reduce = range.] Now that you see these, of course, you notice it involves only the 3-way-interaction coefficients, so in fact any two linearly independent combinations of just the 12th and 13th coefficients will do the job.

There is another way to do this, to get all the hypothesis-testing contrasts, and that is to use the joint_tests() function:

> jt <- joint_tests(model)
> jt
model term                 df1 df2 F.ratio p.value
Sex                          1 107   0.756  0.3865
Group                        2 107   2.190  0.1169
TimeFromBaseline             1 107   0.083  0.7741
Extra                        1 107   4.466  0.0369
Sex:Group                    2 107   0.908  0.4065
Sex:TimeFromBaseline         1 107   0.055  0.8142
Group:TimeFromBaseline       2 107   0.512  0.6009
Sex:Group:TimeFromBaseline   2 107   4.469  0.0137


The jt object is not en emmGrid; it's just an extension of a data frame. But the needed information is in an attribute called est.fcns for "estimable functions":

> attr(jt, "est.fcns")
$Sex (Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB [1,] 0 -0.5258791 0 0 0 0 -0.175293 SexF:GroupC SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline [1,] -0.175293 -0.7359572 0 0 SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline [1,] -0.2453191 -0.2453191$Group
(Intercept) SexF     GroupB     GroupC TimeFromBaseline ExtraBar SexF:GroupB
[1,]           0    0 -0.5200032  0.0000000                0        0  -0.2600016
[2,]           0    0  0.0000000 -0.5200032                0        0   0.0000000
SexF:GroupC SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline
[1,]   0.0000000                     0              -0.7277339               0.0000000
[2,]  -0.2600016                     0               0.0000000              -0.7277339
SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline
[1,]                    -0.363867                     0.000000
[2,]                     0.000000                    -0.363867

$TimeFromBaseline (Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC [1,] 0 0 0 0 -0.8090398 0 0 0 SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline [1,] -0.4045199 -0.2696799 -0.2696799 SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline [1,] -0.13484 -0.13484$Extra
(Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC
[1,]           0    0      0      0                0       -1           0           0
SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline
[1,]                     0                       0                       0
SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline
[1,]                            0                            0

$Sex:Group (Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC [1,] 0 0 0 0 0 0 -0.5813813 0.0000000 [2,] 0 0 0 0 0 0 0.0000000 -0.5813813 SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline [1,] 0 0 0 [2,] 0 0 0 SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline [1,] -0.8136313 0.0000000 [2,] 0.0000000 -0.8136313$Sex:TimeFromBaseline
(Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC
[1,]           0    0      0      0                0        0           0           0
SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline
[1,]             -0.904534                       0                       0
SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline
[1,]                   -0.3015113                   -0.3015113

$Group:TimeFromBaseline (Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC [1,] 0 0 0 0 0 0 0 0 [2,] 0 0 0 0 0 0 0 0 SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline [1,] 0 -0.8944272 0.0000000 [2,] 0 0.0000000 -0.8944272 SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline [1,] -0.4472136 0.0000000 [2,] 0.0000000 -0.4472136$Sex:Group:TimeFromBaseline
(Intercept) SexF GroupB GroupC TimeFromBaseline ExtraBar SexF:GroupB SexF:GroupC
[1,]           0    0      0      0                0        0           0           0
[2,]           0    0      0      0                0        0           0           0
SexF:TimeFromBaseline GroupB:TimeFromBaseline GroupC:TimeFromBaseline
[1,]                     0                       0                       0
[2,]                     0                       0                       0
SexF:GroupB:TimeFromBaseline SexF:GroupC:TimeFromBaseline
[1,]                           -1                            0
[2,]                            0                           -1