Methods vignette

[wlandau]

I propose a codeless vignette that writes out the model in full detail and describes the methodology behind the emmeans-based post processing. I have not done this yet because both brms and emmeans are black boxes, but I hope to understand more. Fully spelled-out descriptions of the methods, plus empirical comparisons with mmrm and SAS, will be really important for helping users in the life sciences decide whether to trust the package.

[wlandau]

One confusing part is the variance components. It might be good to post a brms-only reprex to the brms repo and ask why the sigmas are negative.

[wlandau]

I think I figured it out: for us, we generate the data and formula as follows (with some variations possible on each):

library(brms)
library(tidyverse)
n_group <- 2L
n_patient <- 100L
n_time <- 4L
patients <- tibble(
  group = factor(rep(seq_len(n_group), each = n_patient)),
  patient = factor(seq_len(n_group * n_patient))
)
data <- expand_grid(patients, time = factor(seq_len(n_time))) %>%
  mutate(response = rnorm(n = n()))

data
#> # A tibble: 800 × 4
#>    group patient time  response
#>    <fct> <fct>   <fct>    <dbl>
#>  1 1     1       1       -0.137
#>  2 1     1       2        1.28 
#>  3 1     1       3        0.143
#>  4 1     1       4       -0.834
#>  5 1     2       1       -0.513
#>  6 1     2       2       -1.91 
#>  7 1     2       3        1.51 
#>  8 1     2       4       -1.68 
#>  9 1     3       1        0.978
#> 10 1     3       2        1.18 
#> # ℹ 790 more rows
#> # ℹ Use `print(n = ...)` to see more rows

formula <- brmsformula(
  response ~ time + group + group:time + unstr(time = time, gr = patient),
  sigma ~ 0 + time
)

formula
#> response ~ time + group + group:time + unstr(time = time, gr = patient) 
#> sigma ~ 0 + time

We always set sigma ~ 0 + time. Looking at the Stan code from code <- make_stancode(formula, data), I see a model block:

model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = rep_vector(0.0, N);
    // initialize linear predictor term
    vector[N] sigma = rep_vector(0.0, N);
    mu += Intercept + Xc * b;
    sigma += X_sigma * b_sigma;
    sigma = exp(sigma);
    target += normal_time_het_flex_lpdf(Y | mu, sigma, Lcortime, nobs_tg, begin_tg, end_tg, Jtime_tg);
  }
  // priors including constants
  target += lprior;
}

brms returns b_sigma, which is not what we want. Instead, we want sigma, which we cannot directly access because it is local to the model block. But I dug deeper:

debugonce(brm)
brm(
  data = data,
  formula = formula,
  prior = brms::prior("lkj_corr_cholesky(1)", class = "Lcortime")
)

and I looked at sdata$X_sigma in the debugger (see [where .make_standata() is called):

sdata$X_sigma
#>     time1 time2 time3 time4
#> 1       1     0     0     0
#> 2       0     1     0     0
#> 3       0     0     1     0
#> 4       0     0     0     1
#> 5       1     0     0     0
#> 6       0     1     0     0
#> 7       0     0     1     0
#> 8       0     0     0     1
#> 9       1     0     0     0
#> 10      0     1     0     0
#> 11      0     0     1     0
#> 12      0     0     0     1

And if I change the relevant part of the formula to just sigma ~ time, I get:

sdata$X_sigma
#>     Intercept time2 time3 time4
#> 1           1     0     0     0
#> 2           1     1     0     0
#> 3           1     0     1     0
#> 4           1     0     0     1
#> 5           1     0     0     0
#> 6           1     1     0     0
#> 7           1     0     1     0
#> 8           1     0     0     1
#> 9           1     0     0     0
#> 10          1     1     0     0
#> 11          1     0     1     0
#> 12          1     0     0     1

So it looks like sigma ~ ... indeed controls how b_sigma maps to sigma, and specifying sigma ~ 0 + time creates a 1:1 mapping. All we need to do to get sigma is exp(b_sigma). Indeed, Paul Buerkner confirmed this at for a different model at https://discourse.mc-stan.org/t/retrieve-sigma-on-original-scale-when-sigma-is-a-function-of-a-predictor-sigma-x-1/5058/8

[wlandau]

And then the X matrix straightforwardly corresponds to the first line of the formula.

I think I now know enough to write a methods vignette.

[wlandau]

Actually... I can start, but to be perfectly pedantic I will need to know exactly how emmeans magically transforms posterior samples on main effects into the marginal means, especially how it handles nuisance factors.