Simple model not compiling after 15 minutes, because code is not vectorized?

[potash]

I am trying to fit the multilevel model using brms 2.22.0 and backend cmdstanr 2.36.0:

y ~ x_1 +... x_K + (x_1 + ... + x_K || g1 * g2) 

However, when K = 32 the model does not compile after 15 minutes. Note the following cases do compile quickly:

• replace || with | -- however with K=32 I do not want to fit a 32 x 32 covariance matrix
• K = 2
• replacing g1 * g2 with just g1

Looking at the generated stan code-- when switching from | to ||, things that were vectorized get unvectorized. For example the priors on the random effects involve 32+1 calls to std_normal_lpdf

Is there a way to have uncorrelated effects but still vectorized code?

Here's a script to generate fake data, print the model, and try to compile (you can change || to | to get it to work):

N = 100 # number of observations
K = 32 # number of predictors
J1 = 2 # number of groups g1
J2 = 2 # number of groups g2

df = as.data.frame(matrix(rnorm(N*K), nrow=N))
x_names = paste0("x_", 1:K)
colnames(df) = x_names

df$g1 = factor(sample(1:J1, N, replace=TRUE))
df$g2 = factor(sample(1:J2, N, replace=TRUE))

# y has no relationship to x and g, just fake data for compiling model
df$y = rnorm(N) 

sum_x_str = paste(x_names, collapse="+")
formula = as.formula(sprintf("y ~ %s + (%s||g1*g2)", sum_x_str, sum_x_str))

brms::make_stancode(formula, df, backend="cmdstanr", chains=0)
brms::brm(formula, df, backend="cmdstanr", chains=0)

[paul-buerkner]

Fair point. The code looks pretty intense in places. However, also with | it takes a long time to compile for me. Also there the code looks pretty intense. Yet, I am wondering why compilation time is so slow. Perhaps ask on the stan discourse forums?

In the meantime, if | works for you, you can fix the cholesky factors factor parameters to identity using the constant prior feature and stanvars. See ?constant for related examples

[potash]

Thanks @paul-buerkner for the quick reply. I posted on the discourse here. And thanks for pointing me to the constant feature, this code works for me:

brms::brm(formula, df, backend="cmdstanr", chains=0,
                prior=prior(constant(Id, broadcast=TRUE), class="cor"),
                stanvars= stanvar(diag(1, K+1) , name="Id"),
                )

[paul-buerkner]

I am glad we found a workaround for now :-)

[potash]

Suppose I actually am happy to have a single sd parameter for all the random effects, I can do that like so, but is there a way to put a prior() on the parameter sd_global created by stanvar() or do I need to manually insert it into the stan code?

brms::stancode(formula, df, backend="cmdstanr", chains=1, iter=100,
                prior=c(prior(constant(Id, broadcast=TRUE), class="cor"),
                        prior(constant(sd_global), class="sd")),
                stanvars= c(stanvar(diag(1, K+1) , name="Id"),
                            stanvar(scode="real<lower=0> sd_global;", name="sd_global", block="parameters")),
                )

[paul-buerkner]

You have to add the prior manually also via stanvar by entering it into the transformed parameters block.

[potash]

Whoops, I assumed stanvar was just for variable declarations not arbitrary code. Thanks!