BRMS - use to summarize individual subjects, controlling for some covariates, but not others (ordered factors)

Sorry, this is probably a very trivial question, but I just can’t figure it out.

I want to use BRMS to give me a summary of a person’s data, controlling for some factors, but not others.

I have two models, one for each participant, and one hierarchical one. Each participant contributes ~50 measurements that vary over time and people differ according to a group assignment. I want to use BRMS to extract each person’s mean controlling for time, but not for group.

Here is the model I run separately on each participant:

bf(y ~ 1+time, sigma ~1)

And from this I can get the values of each person using fixef(). This works well in that I get a value for each person that has still their individual variability due to group assignment.

However, now I have simulated data (the ys are used as parameters in another model) and it turned out that parameter recovery was better with a hierarchical model, like this:

bf(y ~1+time+ mo(group) + (1+time|ID), sigma ~ mo(group) + (1|ID))

If I extract the intercepts (and sigma) values of each person with coef(), it gets rid of effects of group. Is there a way to add this back in? When I look at fixef(), there is a value for group, but it’s listed as the same for each person and I’m not sure then how to then multiply it with the categorical (ordered) group assignment.

Many thanks

I am not entirely sure that I understand what you are trying to do, but if you want to get the predicted mean response values and summary from the model, for specific individuals and values of time and group, then you can create a dataframe with those specific characteristics and then use fitted() and the newdata= argument in brms. So, you could make a new dataframe with the value of time, group, and ID that you want, and then feed that into the fitted() function like, fitted(model, newdata=newdata, summary=TRUE), to obtain the predicted mean response with credible intervals.

I thought first I should use fitted(), but then I was unsure how I would get the average for one participant. But now that you’re saying it, maybe this is possible after all, thanks!

What I’m not quite sure about is how to deal with ‘time’ as I want basically the average value for the person (independent of time). So maybe I put time as the average value and the same for everyone and that would get rid of it?

So what I want to get is given that I generated data using

y ~ normal(b_intercept+b_time*list_of_times,b_sigma)

for each person and I fitted this now with the hierachical model (which also has mo(group) as additional factor), what are the outputs of the model that relate back again to b_intercept, b_time and b_sigma?

When I used ‘coef’, there were good correlations among these values with the ones I used to generate the simulations, but it did not include differences between people due to group assignment. (which I’ve not done in the simulation, but if I were to simulate it would be like this:

y ~ normal(b_intercept+b_time*list_of_times+group,b_sigma)

It would help if you posted code for a reproducible example.
But if I am understanding your model correctly, there won’t be “differences between people due to group assignment” in the coef’s, because you modeled group using a monotonic effect in the fixed effects part of the model. Group is not included in the varying effects. The coefficient for the effect of mo(group) in the model should be the same for each person. As it says in the vignette for monotonic effects, “This is realized by parameterizing as follows: One parameter, b, takes care of the direction and size of the effect similar to an ordinary regression parameter. If the monotonic effect is used in a linear model, b can be interpreted as the expected average difference between two adjacent categories of the ordinal predictor.”

Doing a reproducible example was a great suggestion. I’m sorry the code is long, I’ve put the full code in the end in case this description is not clear.

Basically, I understand now how brms’ fitted() works and this is exactly what I wanted. However, I’m still not clear why in my real data my two approaches give different results - as they do not differ in simulations. I think the careful approach here then is to find the result untrustworthy and not report them.

So, what I do is: I simulate 3 groups of participants. Each group differs in their mean and standard deviation (b2 and b5 below). In each group I have 30 participants. Each participant differs in their mean and standard deviation (beyond effects of group) (b0 and b4 below). There is also another covariate of no interest - Age (b3). Each person contributes 50 measurements (‘days’), the measurements can vary linearly over days (b1 below):

p ~ normal(b0 + b1Xday + b2Xgroup1 + b3Xage, sd ~ b4+b5XGroup2)

I analyse the data with two types of models.
[1] Hierarchical model: = brm(bf( y    ~ 1+day+mo(group) + Age + (1+day|ID),
                      sigma~ mo(group) + (1|ID)), data =

I then extract results with fitted(). For this I first generate ‘’ which has the effect of ‘day’ averaged out: = %>%

And then feed this into fitted() to get mean and standard deviation values for each person (not corrected for group or age):

output$  =,,dpar='sigma'))$Estimate

[2] Nonhierarchical model:

df = %>% 
  dplyr::filter(ID==1) # Example for participant with ID=1
fit.singleSubject= brm(y ~ 1+ day,data=df)

And extract parameters like this:

t2 =lapply(,mean)
pars.temp$nonhier.mean = t2$b_Intercept
pars.temp$   = t2$sigma

And finally, I make scatter plots comparing the true value (intercept.mean and with the fitted values (par.mean & for hierarchical model; and nonhier.mean & for nonhierarchical model). I define the true values like this:

temp.sub$intercept.mean = temp.sub$b0 + temp.sub$b1*mean(my.days) + temp.sub$b2 + temp.sub$b3*this.sub.age
temp.sub$ temp.sub$b4 + temp.sub$b5

The result of this is that recovery is perfect and that it’s the same for ther hierachical or non-hierarchical model.


Here is the full code that I’ve used (uploaded is the result of running the code); sorry this is long, I’m still quite far away from full proficiency with R (as well as Stan/BRMS obviously…):

# Simulate some data  -----
# We simulate 3 groups with 30 participants per group and each participant contributes 50 data points
nsub = 30 # subjects per group
my.ranges$b0 =c(-1,1) # intercept
my.ranges$b1 = c(-1,1) # day
my.ranges$b2 = c(-5,0,5) # group 1
my.ranges$b3 = 0.05 # age - fixed effect for each subject
my.ranges$b4 = c(0.6,1) # intercept SD
my.ranges$b5 = c(-0.5,0,0.5) # group SD

age.range = c(20,50)


all.people=data.frame() =data.frame()
for (igroup in 1:3){
  for (isub in 1:nsub){
    temp.sub = list()
    temp.sub$b0 = runif(min=my.ranges$b0[1], max=my.ranges$b0[2],n=1) # intercept
    temp.sub$b1 = runif(min=my.ranges$b1[1], max=my.ranges$b1[2],n=1) # day
    temp.sub$b2 = my.ranges$b2[[igroup]] # group
    temp.sub$b3 = my.ranges$b3 #AGE
    temp.sub$b4 = runif(min=my.ranges$b4[1], max=my.ranges$b4[2],n=1) # SD
    temp.sub$b5 = my.ranges$b5[[igroup]] # group SD
    this.sub.age= runif(min=age.range[1],max=age.range[2],n=1)
    temp.sub$ID = (igroup-1)*nsub + isub
    temp.sub$intercept.mean = temp.sub$b0 + temp.sub$b1*mean(my.days) + temp.sub$b2 + temp.sub$b3*this.sub.age
    temp.sub$ temp.sub$b4 + temp.sub$b5
    # append to subject list
    all.people = rbind(all.people,
    for (iday in my.days){$Age =  this.sub.age$ID = (igroup-1)*nsub + isub$day= iday$group=igroup$y = rnorm(mean = temp.sub$b0 + temp.sub$b1*iday + temp.sub$b2 + temp.sub$b3*$Age,
                        sd   = temp.sub$b4 + temp.sub$b5,
                        n=1) = rbind(,
}$group = factor($group, ordered=TRUE, levels=c(1,2,3))$ID    = factor($ID,ordered=TRUE, levels=unique($ID)) = %>%

# Run BRMS and extract results----
# Hierarchical BRMS model, allowing individual differences in mean and standard deviation = brm(bf( y    ~ 1+day+mo(group) + Age + (1+day|ID),
                 sigma~ mo(group) + (1|ID)),
              data =, silent=TRUE,refresh=0)

# We extract results using brms' fitted(). For this, we generate new data in which for each person, we set the regressor 'day' to the average, to get their average parameter independent of the linear effect of day
# Generate 'new data' to extract results = %>%
names(brmsterms($formula)$dpar) # checking what we want to extract ('mu' is the same as not specifying anything)

# We extract each person's mean and standard deviation
output$ID      = unique($data$ID)
output$group   =$data %>% 
  group_by(ID) %>% dplyr::summarize(group=unique(group)) %>% dplyr::pull(group)
output$  =,,dpar='sigma'))$Estimate

# Single-subject models -----
# For comparison, we also have models that are fitted separately for each person. From these we also extract participants' mean and standard deviation
df = %>% 
fit.singleSubject= brm(y ~ 1+ day,data=df,silent=TRUE,refresh=0)
output.nonHier = data.frame()
for (isub in unique($ID)){ = %>% 
  out = update(fit.singleSubject,,silent=TRUE,refresh=0)
  t2 =lapply(,mean)
  pars.temp= list()
  pars.temp$ID = isub
  pars.temp$nonhier.mean = t2$b_Intercept
  pars.temp$   = t2$sigma
  pars.temp$  = t2$b_day
  output.nonHier = rbind(output.nonHier,

# Store list of data frames if we want to rerun these checks -----
my.checks = list()
all.people$ID = factor(all.people$ID,ordered=TRUE, levels=unique(all.people$ID))
my.checks$output = output

# Compare the results to the simulation ----
to.plot = left_join(my.checks$output,my.checks$output.nonHier,by='ID') %>%
  # add columns for y 
ps$p1=ggscatter(to.plot,x='nonhier.mean',y='par.mean',color='group') # nonhierarchically vs. hierarchically fitted mean
ps$p2=ggscatter(to.plot,x='',y='',color='group') # nonhierarchically vs. hierarchically fitted sd
ps$p3=ggscatter(to.plot,x='intercept.mean',y='par.mean',color='group') # true vs. hierarchically fitted mean
ps$p4=ggscatter(to.plot,x='',y='',color='group') # true vs. hierarchically fitted sd
ps$p5=ggscatter(to.plot,x='intercept.mean',y='nonhier.mean',color='group') # true vs. nonhierarchically fitted mean
ps$p6=ggscatter(to.plot,x='',y='',color='group') # true vs. nonhierarchically fitted sd

brms_checks.RData (43.7 KB)

Ok. I am having a little trouble following your simulation code, but I am not convinced that you are simulating a hierarchical model like you think you are (or perhaps something else is wrong). It looks like you may be creating your varying slopes and intercepts by sampling from a random uniform in a certain range?

However, I’m still not clear why in my real data my two approaches give different results - as they do not differ in simulations.

Unless there is no variation at all based on the level of ID, then the results for each individual taken from the hierarchical model should not be the same as the results of running 90 separate models. When you include varying intercepts and/or varying slopes, your model will be partially pooling information, so the estimates you get from running 90 different individual models on each ID, should not give you the same results as running a hierarchical model and then getting results for each ID. The former assumes no pooling and the latter partial pooling. This vignette from the rstanarm package explains the concept pretty well.
I’m not sure I have time to try to go through all of this sim code, but if you are getting the same results, something doesn’t seem right with the sim…

Thanks! About pooling: I am wondering whether the results are the same for the non-hierarchical and the hierarchical model because as you say the hierarchy in the simulation is a bit weird: the parameters for each person are drawn from a uniform distribution (and the values for each of the measurements for a single person are drawn from a normal distribution).

I did this because usually people (in my field) show plots like what I have that scatter a uniform range of plausible ‘true parameters’ against ‘fitted parameters’. Rather than the true parameters following a normal distribution.

Do you think this is a bad thing and I need to simulate parameters following a normal (rather than uniform distribution) for it to make sense to use a hierarchical model to fit?

EDIT: I have worked out why there were divergent results in my data with different methods:
In the simulation, every participant had the same number of days, so it didn’t really matter that to generate the for fitted(), I had set the day to the average day for each person (rather than setting it to zero to remove it). But it did make a difference in my real data!