I’m using two different remote systems to fit categorical multilevel models, each running a different version of brms (2.17.0 vs 2.18.7) and R (4.2.1 vs 4.2.2).
To my significant consternation, I have now discovered that identical models fit on these two different systems differ significantly in size, with those fit with the older versions twice the size of the newer (~200MB vs ~100MB).
Even more alarmingly, their posterior matrices have different dimensions. The posterior of the older version of brms has about a thousand more parameters. When I look into it, I find that most extra parameters in the “older” model have names like "z_1[1,1]", z_6[1,5]. In addition, the “old” model has three extra parameters which, judging from the names, look like the three population-level intercepts (one for each non-baseline response category), even though the names lack an initial b_ to identify them as population-level.
Despite the differences in size and posterior dimensions, the models have almost idential point estimates and almost identical elpd_loos. Also, the two model formulae are identical, even though identical(formula(old), formula(new)) returns FALSE.
I am therefore assuming that the “difference” is only an internal difference in the computation methods used by the two versions of BRMS, and that I can consider the models identical despite the fact that one is twice the size of the other, with almost twice the number of “parameters”. But it would be reassuring to get a confirmation from a developer.
(I’ll upload the models somewhere if I need to, but the files are huge, plus the data is somewhat sensitive)
If I remember correctly, the z_ parameters may be the standardized group-level effects for the non-centered parameterization for hierarchical models in brms. You could use the stancode() function to check that the two models are the same, but perhaps the latest version of brms doesn’t save the standardized effects but only the actual ones, and thus the file sizes are different? The actual ones are I think named something like r_ but all of this is just going off of memory and just a guess on my part. @paul.buerkner would know for sure
EDIT - I was able to find a brms 2.17.0 and 2.18.5 and try brm(count ~ 1 + Trt + (1|patient), data=epilepsy, cores=4) on both. What I describe above is correct in terms of naming of z_ and r_, but it doesn’t appear that my guess is the solution. The parameters you are describing in your older version model do sound like standardized group level effects. But, I actually can’t see much difference when I run the example model between the two versions. Maybe you should provide a generic brms formula for your model, so someone could try to replicate your scenario.
Thanks for looking into it. I think your diagnosis sounds plausible – basically that the newer version does things more efficiently, resulting in a less cumbersome model object with less output.
Calling stancode() on each model produces almost identical output, but lines 42 to 53 or so, under the heading initialize linear predictor term, are different:
old:
// initialize linear predictor term
// initialize linear predictor term
vector[N] muIndicative = Intercept_muIndicative
vector[N] muIndicative = rep_vector(0.0, N);
+ Xc_muIndicative[start : end] * b_muIndicative;
// initialize linear predictor term
// initialize linear predictor term
vector[N] muShould = Intercept_muShould
vector[N] muShould = rep_vector(0.0, N);
+ Xc_muShould[start : end] * b_muShould;
// initialize linear predictor term
// initialize linear predictor term
vector[N] muHarmony = Intercept_muHarmony
vector[N] muHarmony = rep_vector(0.0, N);
+ Xc_muHarmony[start : end] * b_muHarmony;
// linear predictor matrix
// linear predictor matrix
array[N] vector[ncat] mu;
Given that the results are practically the same, I’m pretty well satisfied that there’s no problem. It was just initially a shock to see such different sizes and posterior matrices.
In any case, I don’t think there is anything to worry about. I still wonder if it has to do with brms 2.17.0 saving the standardized group level effects, like rstan does. For example, if you run:
library(brms) #version 2.18.5
m1<-brm(count ~ 1 + Trt + (1|patient), data=epilepsy)
#create the same stan code and data as the above model to run in rstan
sc <- stancode(m1)
sda <- standata(m1)
library(rstan) #rstan version 2.26.13 (Stan version 2.26.1)
fit1 <- stan(model_code = sc, data = sda, cores = 4)
#compare parameters saved by brms and rstan
s1 <- as_draws_df(m1)
str(s1)
print(fit1)
#brms only saved the actual group-level effects (denoted by 'r_'), but rstan saved both the standardized and the actual (denoted by 'r_' and 's_')
#compare file size (note there will be some difference due to brms saving additional things like the dataset
object.size(m1)
object.size(fit1)
#around twice the size for rstan
However, I can’t see that brms does this for the 2.17 version that I ran, so not sure if that really is the solution.
Yep I think you’ve identified the difference. When I look at all those “extra” parameters named z[1,1] and so on, their number corresponds to the number of groups in each random-effect parameter. And the model fit by version 2.18.7 dispenses with the z_ terms but retains the r_ params, while the older version has both. That must be it. Thanks for your help!