Group-level effect on the response scale as conditional effect of a new group level?

Interfaces brms
, ,
1

Dear all,
I recently came across this interesting blog by Andrew Heiss, which nicely explains marginal effects and conditional effects for existing and new groups (A guide to correctly calculating posterior predictions and average marginal effects with multilievel Bayesian models | Andrew Heiss). He gave me some interesting ideas for my own analyses. However, I’m not quite sure if my interpretation and application of conditional effects for new groups is correct. So I would be very grateful for any feedback.

I am currently working with several Brms multilevel models with non-Gaussian response distributions (e.g. family dirichlet, lognormal, hurdle_lognormal). I’m struggling with how to communicate the magnitude of the group-level effects in a meaningful way for “regular readers”/non-statisticians. Just stating the sd parameter on the link scale is not really easy to interpret or at all meaningful when it stands alone.

My models all have the general form of:
Y ~ treatment+ (treatment|species).
Where treatment is a two-level factor as a fixed effect and a random slope and intercept for different species (n=35). I included my species as random factors because I am not interested in differences between individual species, but in quantifying the overall variability of the treatment effect across the range of potential species. Thus, I am mainly interested in interpreting the variability of the random slopes.

I will report the population treatment effect as a marginal effect by posterior predictions as mean and 95% confidence interval.
Now I wonder. if it would be possible to specify the variability added by my group level parameters as the confidence interval of a new group level created by sample_new_levels. If I understand the above blog correctly, simulating a new group level results in the following:
a. sample_new_levels = “uncertainty”: a conditional treatment effect for a new species that summarizes the variability of all the species I studied.
b. sample_new_levels = “gaussian”: a conditional treatment effect for a new species that might exist as suggested by my data.
Would it therefore be valid to communicate the increase in confidence interval between the marginal effect and the conditional effect of a new group as variability contributed by my group-level parameter “species”?

So, is it possible to derive hypothetical statements like the following from these comparisons? "On average, Y was +5 [3;5] higher in treatment A than in treatment B (marginal effect). However, the treatment effect was clearly modified by the species studied, with predicted treatment effects ranging from +2 to +6 (95% CI of sample_new_levels = “uncertainty”). Alternatively: The treatment effect was substantially modified by species, since it is predicted to range from +0.8 to +8.2 within the individual species (95% CI of sample_new_levels = “gaussian”).

If this is a generally valid approach, is it sufficient to obtain posterior predictions for one new group level, or would it be better to obtain predictions for, say, 100 new groups and summarize their posterior predictions? If this is not the right way, are there other ways I can easily describe group level effects in terms of their magnitude on the response scale?
Sorry for the long post. Any feedback would be greatly appreciated.

Cheers.

2

Here’s a short example of what the code would look like:

###Marginal treatment effects

#calculate posterior predictions
t<-fitted(model,
            newdata = data.frame(Treatment = c("A", "B")),
            summary = FALSE,
            re_formula=NA)
dim(t)
# 8000    2
#Column 1: Treatment A, Column 2: Treatment B

#Calculate difference between treatments
marg_effect<- t[, 2] - t[,1]
marg_effect<-c(
  mean = mean(marg_effect),
  quantile(marg_effect, c(0.025, 0.975)))


#####Conditional treatment effect

#calculate posterior predictions for a new group-level 
new_phen<-fitted(fit_zip8,
                 newdata = data.frame(Treatment = c("A", "B"),
                                      Genotype=rep("Brandnew_phenotype",each=2)),
                 re_formula = NULL, 
                 allow_new_levels = TRUE,
                 sample_new_levels = "gaussian",
                 summary = FALSE)
dim(new_phen)
# 8000    2

#Calculate difference between treatments 
con_effect<-new_phen[, 2] - new_phen[,1]
con_effect<-c(
  mean = mean(con_effect),
  quantile(con_effect, c(0.025, 0.975)))

###Result
round(marg_effect,2)
# mean   2.5%   97.5%   
#-0.76  -0.98  -0.54  

round(con_effect,2)
#mean    2.5%   97.5%  
#-0.75   -1.75  0.31

Why are in the above described approach the means of the conditional and marginal treatment effects the same when the response distribution is Gaussian, but not for other response distributions?