The model

A multigroup model has one coefficient matrix per group. Write

• \(n_{\text{groups}}\) for the number of groups,
• \(q\) for the number of repeated measures,
• \(p + 1\) for the number of predictors, the intercept included,
• \(\mathbf{B}_g \in \mathbb{R}^{(p + 1) \times q}\) for the coefficients of group \(g\).

The \(\mathbf{B}_g\) are stacked into a three-dimensional array of dimension \((p + 1) \times q \times n_{\text{groups}}\). In semnova that array is the contrast (pi) layer’s coefficients, returned by regression_coefficients().

The hypothesis

The hypotheses of interest have the form

\[ \mathbf{H}_0: \mathbf{G} \cdot (\mathbf{L} \cdot \mathbf{B}_g \cdot \mathbf{M}) = \mathbf{K} \]

where each matrix selects one aspect of the design:

A main effect, an interaction and a covariate effect differ only in which of these is non-trivial, which is why one implementation covers them all.

Vectorisation

The hypothesis is stated in matrices but tested as a linear restriction on a parameter vector, so it has to be rewritten in terms of \(\text{vec}(\mathbf{B})\).

Dimensions

Implementation

contrast_matrix() builds \(\mathbf{R}\) in the three steps derived above:

contrast_matrix
#> function (G, L, M, B_array = NULL) 
#> {
#>     if (!is.matrix(G) || !is.matrix(L) || !is.matrix(M)) {
#>         stop("semnova: G, L, and M must all be matrices", call. = FALSE)
#>     }
#>     h <- nrow(L)
#>     r <- ncol(M)
#>     d <- nrow(G)
#>     n_groups <- ncol(G)
#>     Rg <- kronecker(t(M), L)
#>     R_within <- kronecker(diag(n_groups), Rg)
#>     R <- kronecker(G, diag(h * r)) %*% R_within
#>     dim_out <- c(rows = d * h * r, cols = n_groups * ncol(L) * 
#>         nrow(M))
#>     names(dim_out) <- c("rows", "cols")
#>     if (!is.null(B_array)) {
#>         if (dim(B_array)[1] != ncol(L) || dim(B_array)[2] != 
#>             nrow(M) || dim(B_array)[3] != n_groups) {
#>             stop("semnova: B_array dimensions (", paste(dim(B_array), 
#>                 collapse = " x "), ") are not compatible with L/M/G (expected ", 
#>                 ncol(L), " x ", nrow(M), " x ", n_groups, ")", 
#>                 call. = FALSE)
#>         }
#>         return(list(R = R, dim = dim_out, vec_B = as.vector(B_array)))
#>     }
#>     list(R = R, dim = dim_out)
#> }
#> <bytecode: 0x5d5de133c0c8>
#> <environment: namespace:semnova>

Supplying B_array returns \(\text{vec}(\mathbf{B})\) alongside it, so the two halves of the restriction are produced by one call and cannot disagree about their ordering:

set.seed(1)

# 3 predictors, 4 repeated measures, 2 groups
B_array <- array(rnorm(3 * 4 * 2), dim = c(3, 4, 2))

L <- matrix(c(0, 1, 0), nrow = 1)          # one predictor contrast
M <- cbind(c(1, -1, 0, 0), c(0, 0, 1, -1)) # two within-subject contrasts
G <- matrix(c(1, -1), nrow = 1)            # one group contrast

out <- contrast_matrix(G = G, L = L, M = M, B_array = B_array)
out$dim
#> rows cols 
#>    2   24

These are the dimensions tabulated above: \(d \cdot h \cdot r = 2\) rows and \(n_{\text{groups}} \cdot (p + 1) \cdot q = 24\) columns.

The construction can be checked against its own definition. Computing \(\mathbf{G} \cdot (\mathbf{L} \mathbf{B}_g \mathbf{M})\) directly, group by group, must reproduce \(\mathbf{R} \cdot \text{vec}(\mathbf{B})\):

direct <- Reduce(`+`, lapply(1:2, function(g) {
    G[1, g] * as.vector(L %*% B_array[, , g] %*% M)
}))

all.equal(as.vector(out$R %*% out$vec_B), direct)
#> [1] TRUE

This vignette is rendered when the package is built, so that equality is verified on every build rather than asserted here once.

Where each part comes from

The restriction \(\mathbf{R} \cdot \text{vec}(\mathbf{B}) - \mathbf{k}\) is what those functions test, and the same quantity underlies the Wald test reported in the ANOVA table, the per-contrast estimates in summary(), and the intervals drawn by plot().