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

**URL:** <https://discourse.mc-stan.org/t/prior-predictive-check-in-distributional-mixed-model-including-uniform-default-priors/27527>\
**Category:** Modeling\
**Tags:** prior-predictive, brms\
**Created:** [May 18, 2022, 11:08am UTC](https://discourse.mc-stan.org/t/prior-predictive-check-in-distributional-mixed-model-including-uniform-default-priors/27527 "2022-05-18T11:08:08Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Theresa\_U](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/theresa_u/32/15539_2.png) [@Theresa\_U](https://discourse.mc-stan.org/u/Theresa_U)\
**Post date:** [May 18, 2022, 11:08am UTC](https://discourse.mc-stan.org/t/prior-predictive-check-in-distributional-mixed-model-including-uniform-default-priors/27527/1 "2022-05-18T11:08:08Z")

</div>

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](https://cran.r-project.org/web/packages/brms/vignettes/brms_distreg.html)).

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(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.

---

<div class="post-metadata">

**Author:** ![Theresa\_U](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/theresa_u/32/15539_2.png) [@Theresa\_U](https://discourse.mc-stan.org/u/Theresa_U)\
**Post date:** [May 25, 2022, 1:44pm UTC](https://discourse.mc-stan.org/t/prior-predictive-check-in-distributional-mixed-model-including-uniform-default-priors/27527/2 "2022-05-25T13:44:24Z")

</div>

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

```r
## 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<-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") 
``
```