Understanding the parameters of a hierarchical Dirichlet regression

Dear all,

Sorry if my question seems basic and not well stated - I am fairly new to working with Brms and only have basic statistical knowledge. Mainly I am struggling to get my model parameters (group and population level) from Dirichlet regression in an interpretable form on the probability scale.
I am running a Dirichlet regression on Brms for a response consisting of a composite over three particle size classes (y1, y2, y3). My predictor is a treatment variable with two levels (treatmentA & treatmentB) depending on group effects (27 groups).

brm(bind(y1,y2,y3) ~Treatment+(1+Treatment|Group) , data=data, dirichlet(refcat = “y3”),
chains = 4, iter = 4000)

The Dirichlet regression gives me the output shown below. I understand that all population parameters in the output are given as deviations from my reference category y3. However, the main purpose of my model is to quantify the impact of treatment on each dependent variable and to determine how much each group differs from the overall population. Therefore, I would like to extract my population and group parameters, if possible as probabilities and independently from the values of the reference category. Is there by any chance a way to do this?
Is setting refcat=NA the only option to obtain estimates for each of the three responses?

I would be very grateful for any recommendations.
Many thanks!

Family: dirichlet
Links: muy1 = logit; muy2 = logit; phi = identity
Formula: bind(y1, y2, y3) ~ Treatment + (1 + Treatment | Group)
Data: data (Number of observations: 80)
Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup draws = 8000

Group-Level Effects:
~Group (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(muy1_Intercept) 0.06 0.04 0.00 0.17 1.00 2095 3906
sd(muy1_TreatmentTreatmentB 0.15 0.10 0.01 0.36 1.00 1405 2306
sd(muy2_Intercept) 0.06 0.04 0.00 0.16 1.00 2571 3929
sd(muy2_TreatmentTreatmentB) 0.19 0.11 0.02 0.43 1.00 1436 2209
cor(muy1_Intercept,muy1_TreatmentTreatmentB) -0.39 0.52 -0.98 0.84 1.00 2538 3051
cor(muy2_Intercept,muy2_TreatmentTreatmentB) -0.32 0.54 -0.98 0.87 1.00 1757 3289

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
muy1_Intercept 2.94 0.05 2.84 3.05 1.00 5660 4115
muy2_Intercept 1.33 0.06 1.22 1.44 1.00 6802 6038
muy1_TreatmentTreatmentB 0.64 0.09 0.46 0.84 1.00 4121 2175
muy2_TreatmentTreatmentB 0.04 0.11 -0.18 0.26 1.00 4115 4141

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 283.37 33.63 221.02 353.89 1.00 7092 6126

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

1 Like

Hi, and sorry for the late response. posterior_epred will be yur friend in a model like this. For example:

df <- data.frame(y1 = runif(1000))
df$y2 <- runif(1000, 0, 1 - df$y1)
df$y3 <- 1 - df$y1 - df$y2
df$g <- c("a", "b", "c", "d")

mod <- brm(
  cbind(y1,y2,y3) ~ g, 
  data = df, 
  family = dirichlet(refcat = "y3"),
  backend = "cmdstanr",
  chains = 1
)

a <- posterior_epred(mod, newdata = data.frame(g = c("a", "b", "c", "d")))
dim(a)

#output for the first posterior iteration
a[1, , ]

# Note that probabilities always sum to 1 for every group
rowSums(a[1, ,])

Dear Jacob,

Thank you very much for your reply and suggestion. I have started to look into the posterior_epred function and it seems that it could be the solution to all my problems. Since I’ve been looking into it, I’ve come across a new question that is probably pretty basic again - still, I’m not entirely sure about it.

Is it a valid approach to apply non-linear hypothesis testing to the output of posterior_epred? When I apply the brms “hypothesis” function to the output of posterior_epred, I get results that look pretty reasonable in terms of numbers/values and are exactly what I was looking for, but I’m not sure if I’m not completely contradicting the intended purpose of the “hypothesis” function.

If possible, I would like to investigate the overall effect of the treatment by comparing the two treatments at the population level and further investigate how much/how often each group deviates from the overall population mean. Below is the code I tried for the model described above.

My main motivation for this approach is to try to make the results of the data analyses for these responses comparable to the other data I have from the same experiment. Most of them are simple continuous variables and I analyzed them using a mixed model with the same predictors and Gaussian distribution, e.g. brm(bind(cont1,cont2,cont3)~Treatment+(1+Treatment|Group) , data=data2, chains = 4, iter = 4000). Since this is a Gaussian model, I can easily work directly with the estimated model parameters. Thus, I can apply linear hypothesis tests to the treatment parameter (at the population level) and get an idea of the group differences from the overall population by looking at the posterior distributions of the random effects.

Is the intended use of posterior_epred and hypothesis acceptable or have I accidentally created a statistical monstrosity? If you can give me any feedback on this issue, I would be very grateful.

Code:
m<- brm(bind(y1,y2,y3) ~Treatment+(1+Treatment|Group) , data=data, dirichlet(refcat = “y3”), chains = 4, iter = 4000)

######Predictions on population level
pop_pred<- posterior_epred(m, newdata = data.frame(Treatment = c("treatmentA” , “treatmentB ")),re_formula = NA)

dim(pop_pred)
[1] 8000 2 3 #8000 posterior draws of the two treatments for the three responses

#hypothesis test treatment effect response y1
a<-pop_pred[,1]
a<-as.data.frame(a)
names(a)<- c("treatmentA” , “treatmentB ")

test<- hypothesis(a,hypothesis = c(“treatmentA - treatmentB = 0”))

test
Hypothesis Tests for class :
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (A-B) = 0 -0.08 0.01 -0.09 -0.07 NA NA *

‘CI’: 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
‘*’: For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

#######Predictions on group level
#Create dataframe with all combinations of factor level

new_dat<- data.frame(Treatment=c(“treatmentA” , “treatmentB”), Group=rep(Group,each=2) )
group_pred<- posterior_epred(m, newdata =new_dat)

dim(group_pred
[1] 8000 54 3 )
#8000 posterior draws of the 54 factor combinations (two treatments*27 groups) for the three responses

#hypothesis test response Y1
b<-group_pred[,1]
b<-as.data.frame(b)
names(b)<- c("Group1treatmentA” , “Group1treatmentB "……)
c<- cbind (a,b) #combine population predictions with group predictions in one data.frame

#Test if individual group differs from overall population prediction (similar to intercept in gaussian model: Treatment A)
test2<- hypothesis(c, hypothesis = c(“Group1treatmentA - treatmentA = 0”))

#test2
Hypothesis Tests for class :
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (Group1A - A) = 0 -0.04 0.02 -0.07 -0.01 NA NA *

‘CI’: 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
‘*’: For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

#Continue with other groups

#Test if individual group differs from overall population prediction (here parameter similar to slope in gaussian models: difference between Treatment A-TreatmentB)

#create “posterior” values of population slope
test<- hypothesis(a,hypothesis = c(“treatmentA - treatmentB = 0”))
pop_slope<-test$samples #dataframe with 8000 observations

#create “posterior” values of slope Group1
test3<- hypothesis(c, hypothesis = c(“Group1treatmentA - Group1treatmentB = 0”))
Group1_slope<-test3$samples #dataframe with 8000 observations

#Test if slope Group1 differs from overall population slope
slope<-cbind(pop_slope,Group1_slope)
test4<- hypothesis(slope,hypothesis = c(“pop_slope - Group1_slope = 0”))

#Continue with other groups