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:
- 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 - Are prior predictve cheks a sensible approach in this context?
- Should I also model population level (fixed) effects or only the standard deviations of group-level (‘random’) effects?
- brms uses improper priors as default priors for the population level effects. How can I simulate these effects?
- 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?
- Is it correct that I take the exponent of the model for \sigma?
library(extraDistr)
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)
scale_default<-ifelse(mad_df>2.5,mad_df,2.5)
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(
rlst(Nsim,df=3,mu=0,sigma=log(mad_df))+
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.