Hi all.
Is it possible to get the model matrices for both fixed and random effects from a fitted brms model?
Ideally, something like the output of the following, for lmer models:
mm.fixed <- getME(model, "X")
mm.random <- getME(model, "mmList")
I’ve also tried model.matrix(model), which produces an error: “no terms component nor attribute”.
Thanks!
Since brms syntax often mimics lme4 syntax, one workaround might be to fit the model in lme4 and then take the model matrix from that model. This works only if the brms model doesn’t have any brms-specific syntax.
You can use make_standata() and then get all model matrices. the fixed effect matrix is called X, the random effect ones are individual vectors beginning with Z:
library("brms")
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.18.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#>
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#>
#> ar
data("Machines", package = "MEMSS")
#m1 <- brm(score ~ Machine + (Machine|Worker), data=Machines)
d1 <- make_standata(score ~ Machine + (Machine|Worker), data=Machines)
str(d1)
#> List of 12
#> $ N : int 54
#> $ Y : num [1:54(1d)] 52 52.8 53.1 51.8 52.8 53.1 60 60.2 58.4 51.1 ...
#> $ K : int 3
#> $ X : num [1:54, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> .. ..$ : chr [1:3] "Intercept" "MachineB" "MachineC"
#> ..- attr(*, "assign")= int [1:3] 0 1 1
#> ..- attr(*, "contrasts")=List of 1
#> .. ..$ Machine: num [1:3, 1:2] 0 1 0 0 0 1
#> .. .. ..- attr(*, "dimnames")=List of 2
#> .. .. .. ..$ : chr [1:3] "A" "B" "C"
#> .. .. .. ..$ : chr [1:2] "B" "C"
#> $ Z_1_1 : num [1:54(1d)] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> $ Z_1_2 : num [1:54(1d)] 0 0 0 0 0 0 0 0 0 0 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> $ Z_1_3 : num [1:54(1d)] 0 0 0 0 0 0 0 0 0 0 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> $ J_1 : int [1:54(1d)] 1 1 1 2 2 2 3 3 3 4 ...
#> $ N_1 : int 6
#> $ M_1 : int 3
#> $ NC_1 : int 3
#> $ prior_only: int 0
#> - attr(*, "class")= chr [1:2] "standata" "list"
head(d1$X)
#> Intercept MachineB MachineC
#> 1 1 0 0
#> 2 1 0 0
#> 3 1 0 0
#> 4 1 0 0
#> 5 1 0 0
#> 6 1 0 0
options(contrasts=c('contr.sum', 'contr.poly'))
d2 <- make_standata(score ~ Machine + (Machine|Worker), data=Machines)
str(d2)
#> List of 12
#> $ N : int 54
#> $ Y : num [1:54(1d)] 52 52.8 53.1 51.8 52.8 53.1 60 60.2 58.4 51.1 ...
#> $ K : int 3
#> $ X : num [1:54, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> .. ..$ : chr [1:3] "Intercept" "Machine1" "Machine2"
#> ..- attr(*, "assign")= int [1:3] 0 1 1
#> ..- attr(*, "contrasts")=List of 1
#> .. ..$ Machine: num [1:3, 1:2] 1 0 -1 0 1 -1
#> .. .. ..- attr(*, "dimnames")=List of 2
#> .. .. .. ..$ : chr [1:3] "A" "B" "C"
#> .. .. .. ..$ : NULL
#> $ Z_1_1 : num [1:54(1d)] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> $ Z_1_2 : num [1:54(1d)] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> $ Z_1_3 : num [1:54(1d)] 0 0 0 0 0 0 0 0 0 0 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:54] "1" "2" "3" "4" ...
#> $ J_1 : int [1:54(1d)] 1 1 1 2 2 2 3 3 3 4 ...
#> $ N_1 : int 6
#> $ M_1 : int 3
#> $ NC_1 : int 3
#> $ prior_only: int 0
#> - attr(*, "class")= chr [1:2] "standata" "list"
head(d2$X)
#> Intercept Machine1 Machine2
#> 1 1 1 0
#> 2 1 1 0
#> 3 1 1 0
#> 4 1 1 0
#> 5 1 1 0
#> 6 1 1 0
Created on 2022-11-09 with reprex v2.0.2
To clarify, this will go inside a function that only takes a brms model as input (so I’d rather not have to refit the model, Jeremy).
But going from make_standata() should work fine! E.g.,
m1 <- brm(score ~ Machine + (Machine|Worker), data=Machines)
d1 <- make_standata(formula(m1), m1$data)
Thanks!