# Prior predictive check in distributional mixed model including uniform default priors

Hi,

I am trying to do prior predictive checks to compare my own priors to the default priors in brms. I have a distributional model with random effects. Below you find a minimal example adapted from the brms manual (Estimating Distributional Models with brms).

My model of the patients symptoms should look like this:

y_{jk}=\alpha+\beta \cdot 1(\text{group}_{jk}=\text{treatment}) +\tau_j+\epsilon_{jk}

With random effects of 5 hospitals \gamma_j \sim \mathcal{N}(0,\sigma_{\tau}) and residuals with treatment specific standard deviations \epsilon_{jk} \sim \mathcal{N}(0,\sigma_\text{group}) where

\sigma_\text{group}=\exp(\psi+\lambda \cdot 1(\text{group}_{jk}=\text{treatment}) ) are modelled in brms via the log-link function as an exponent of a linear predictor.

I have a couple of questions:

1. I have some problems in matching the output of the default prior via get_prior() with my formula.
If I understand it correctly, then
• The line starting with “student_t(3, 1, 2.5) Intercept” represents \alpha in my model.
• The line “(flat) b Grouptreat” with an empty “dpar” is \beta in my model and \lambda for the line with “dpar”=sigma.
• The line starting with “student_t(3, 0, 2.5) Intercept” and “dpar”=sigma represents \psi in my model, the logarithm of the mean standard deviation in the first group.
• The line starting with “student_t(3, 0, 2.5) sd Intercept Hospital” is \sigma_{\tau} in my model.
Are these interpretations correct and to what model parameters do the other lines with class =”sd” and class=”sigma” belong to
2. Are prior predictve cheks a sensible approach in this context?
3. Should I also model population level (fixed) effects or only the standard deviations of group-level (‘random’) effects?
4. brms uses improper priors as default priors for the population level effects. How can I simulate these effects?
5. If I understand it correctly, then Stan samples variance parameters from truncated distributions. Do I have to respect this in my prior predictive simulation somehow?
6. Is it correct that I take the exponent of the model for \sigma?

library(crch)
set.seed(seed = 9048203)

## Examplary Data ---------------------------------------

n_per_group<-30
n_hospital<-5
group_effect<-1
hospital_sd<-group_effect/4

groups <-c(“treat”, “placebo”)
N<-length(groups)*n_per_group

group<-rep(groups,each=n_per_group)
hospital<-factor(sample(1:n_hospital,N,replace=T))

hospital_effect<-rnorm(n_hospital,0,hospital_sd)
symptom_post ← c(rnorm(30, mean = group_effect, sd = 2), rnorm(30, mean = 0, sd = 1))+
hospital_effect[hospital]

df<-data.table(“ID”=1:N,
“Group”=group,
“Hospital”=hospital,
“Symptom_post”=symptom_post
)

reg.form<-bf(Symptom_post~Group+(1|Hospital),sigma~Group)

get_prior(reg.form,data=df) # Default brms priors

## Prior predictive checks ---------------------------

Nsim<-nrow(df)
mad_df<-round(mad(df$Symptom_post),2) mean_default<-round(median(df$Symptom_post),1)

## Default prior

alpha<-rnorm(1,mean = mean_default,sd=scale_default) # intercept
tau<-rlst(length(unique(df$Hospital)),df=3,mu=0,sigma=scale_default) # hospital effect sigma_group_treat<-exp(rlst(Nsim,df=3,mu=0,sigma=scale_default)) # model of standard deviation ## ??question 4: how to add flat priors for group effect in expected value and standard devation sim = alpha+ tau[df$Hospital]+
sigma_group_treat

data_default ← data.frame(
Symptom_post = df$Symptom_post, Sim=sim ) ggplot(data_default, aes(x = symptom_post, y = Sim)) + geom_point(alpha = 0.1, color = “red”) ## Weakly-informative prior alpha<-rnorm(1,mean = mean_default,sd=mad_df) # intercept beta<-rnorm(1,mean=0,sd=mad_df) # agnostic (zero-centered) prior for the effect tau<-rnorm(length(unique(df$Hospital)),mean=0,sd=mad_df) # hospital effect

beta_sd<-rnorm(1,mean=0,sd=mad_df/4) # effect of group on standrd deviation
sigma_group_treat<-exp(
ifelse(df$Group==“treat”,beta_sd,0) ) summary(sigma_group_treat) sim = alpha+ ifelse(df$Group==“treat”,beta,0)+
tau[df$Hospital]+ sigma_group_treat data_wip ← data.frame( Symptom_post = df$Symptom_post,
Sim=sim
)

ggplot(data_wip, aes(x = symptom_post, y = Sim)) +
geom_point(alpha = 0.1, color = “red”)
}

Thank you a lot for your help and for reading my question and rather long code.

I realized by now that I misinterpreted the output of get_prior() so I think I should change the simulations to the following:

## Prior predictive checks ---------------------------
Nsim<-nrow(df)
mad_df<-round(mad(df$Symptom_post),2) mean_default<-round(median(df$Symptom_post),1)

## Default prior
alpha<-rlst(1,df=3,mu = mean_default,sigma=scale_default) # intercept
psi<-rlst(1,df=3,mu=0,sigma=scale_default) # log standard deviation in control group

infinity<-10^100
beta<-runif(1,min=-infinity,max=infinity) # representing flat prior for the effect (although infinity is not really infinity)
lambda<-runif(1,min=-infinity,max=infinity)
## ??question 4: how to add flat priors for group effect in expected value and standard devation

sigma_tau<-abs(rlst(1,df=3,mu=0,sigma=scale_default)) # sd of hospital effect  (half t)

tau<-rnorm(length(unique(df$Hospital)),0,sigma_tau) epsilon<-rnorm(Nsim,0,exp(psi+ifelse(df$Group=="treat",lambda,0)))

sim = alpha+
ifelse(df$Group=="treat",beta,0)+ tau[df$Hospital]+
epsilon

data_default <- data.frame(
Symptom_post = df$Symptom_post, Sim=sim ) ggplot(data_default, aes(x = symptom_post, y = Sim)) + geom_point(alpha = 0.1, color = "red") ## Weakly-informative prior alpha<-rnorm(1,mean = mean_default,sd=mad_df) # intercept psi<-rnorm(1,mean = mean_default,sd=mad_df) # log standard deviation in control group beta<-rnorm(1,mean=0,sd=mad_df) # agnostic (zero-centered) prior for the effect lambda<-rnorm(1,mean=0,sd=mad_df) sigma_tau<-abs(rnorm(1,0,scale_default)) # sd of hospital effect (half normal) tau<-rnorm(length(unique(df$Hospital)),0,sigma_tau)
epsilon<-rnorm(Nsim,0,exp(psi+ifelse(df$Group=="treat",lambda,0))) sim = alpha+ ifelse(df$Group=="treat",beta,0)+
tau[df$Hospital]+ epsilon data_wip <- data.frame( Symptom_post = df$Symptom_post,
Sim=sim
)

ggplot(data_wip, aes(x = symptom_post, y = Sim)) +
geom_point(alpha = 0.1, color = "red")