# Simulating Data Using

**URL:** <https://discourse.mc-stan.org/t/simulating-data-using/5881>\
**Category:** brms\
**Created:** [October 4, 2018, 6:34pm UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881 "2018-10-04T18:34:54Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![jackbailey](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jackbailey/32/14978_2.png) [@jackbailey](https://discourse.mc-stan.org/u/jackbailey)\
**Post date:** [October 4, 2018, 6:34pm UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881/1 "2018-10-04T18:34:54Z")

</div>

In order to overcome some of my bad habits, I’ve begun to simulate data before fitting a model to my actual data.

Where the models are simple, this is fine. But it can become arduous quickly. In the case of some of brms’ functions, I am aware of what they are for, but have no clue how I could simulate data for them. This problem compounds itself when more than one of these functions are included in the model.

The obvious solution is: learn how to simulate the data. Yes, that’s true and worth doing. But there must still be a case for a convenience function that does the same job.

My question, then, is whether it’s possible to specify arbitrary parameters for any brms model, for example by specifying priors for each parameter in the model, from which I can simulate a data set for later use?

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [October 4, 2018, 6:42pm UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881/2 "2018-10-04T18:42:52Z")

</div>

The workflow is to specify priors and set `sample_prior = "only"` in `brm`. Then, use `posteror_predict` to obtain simulated responses.

---

<div class="post-metadata">

**Author:** ![jackbailey](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jackbailey/32/14978_2.png) [@jackbailey](https://discourse.mc-stan.org/u/jackbailey)\
**Post date:** [October 4, 2018, 7:03pm UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881/3 "2018-10-04T19:03:51Z")

</div>

Thanks, Paul. That’s excellent.

---

<div class="post-metadata">

**Author:** ![Jacquie](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jacquie/32/19488_2.png) [@Jacquie](https://discourse.mc-stan.org/u/Jacquie)\
**Post date:** [June 29, 2020, 9:22am UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881/4 "2020-06-29T09:22:45Z")

</div>

Hi there, I’m trying this out at the moment as well. For a **hierarchical** model. Something is not right… Is it meant to work for hierarchical models? I’ve put something like this:  
To make the priors:

```
bprior <- c(
prior(normal(-0.12,0.04),class = "Intercept"),
prior(normal(-0.14,0.01),class="b",coef="Coef1"),
prior(normal(0.66,0.03),class="sd",coef="Intercept",group="ParticipantID"),
prior(normal(0.07,0.01),class="sd",coef="Coef1",group="ParticipantID"),
prior(normal(0.58,0.01),class="sigma"))

```

Then I’ve made simulations like this:

```
simulModel= brm(y~1+Coef1+(1+Coef1|ParticipantID),data=df,sample_prior="only")

newY=predict(simulModel)

```

And finally I tried to fit simulation like this (first using adding newY to data frame df-\>‘newdf’)

```
fit= brm(newY ~ 1+Coef1+(1+Coef1|ParticipantID),data=newdf)

```

From the results it looks like the parameters that were used to make the simulation were without any noise:  
sd(Intercept) and sd(Coef1) is 0.01.  
Est.Error on Intercept and Coef1 are 0 and sigma is 0.01

Did I do something wrong trying to make the simulation or would this anyway not work for a hierarchical simulation? What I would try to do coding manually would be to first draw regression weights for each person from the normal distribution across people and then create predictions for each measurement per person.

---

<div class="post-metadata">

**Author:** ![julianquandt](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/julianquandt/32/5886_2.png) [@julianquandt](https://discourse.mc-stan.org/u/julianquandt)\
**Post date:** [July 8, 2020, 2:23pm UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881/5 "2020-07-08T14:23:07Z")

</div>

Hi Jacquie,

I think your approach should work for hierarchical models as well. One thing i noticed is that you dont use your priors in `simulModel` which is probably just a typo.

I think what might be the problem here might be the `newY = predict(simulModel)` uses the default `summary = TRUE` which gives the means of all samples and therefore does already average out the sampling-noise, which is probably why your parameters do not show the correct `SE`. Changing this to `newY=predict(simulModel, newdata = df, summary=FALSE, nsamples=1)[1,]` should do the trick.

> **Long answer + example**
>
> **Data generation**
> 
> ```no-highlight
> library(sjPlot)
> library(ggplot2)
> library(MASS)
> library(brms)
> library(plyr)
> library(gridExtra)
> 
> # data generation function
> gen_design <- function(n_participant, n_stimulus, n_levels){
> participants <- rep(1:n_participant, each = n_stimulus)
> #stims <- rep(1:n_stimulus, times = n_participant)
> x <- c()
> stims <- as.vector(sapply(seq(1, n_participant), function(x) sample(seq(1,n_stimulus), size = n_stimulus)))
> for(i in 1:n_participant){
> x <- append(x, seq(0,(n_levels-1)))
> #x <- append(x, sapply(seq(1,n_stimulus), function(x) sample(seq(1,n_levels), size = n_levels)))
> }
> expdesign <- data.frame(participant = participants, stimulus = stims, x = x)
> return(expdesign)
> }
> 
> sp <- list(n_participant = 30, n_stimulus = 20, n_levels =2) # data with 30 participants each seeing 20 stimuli (which will not be modelled in this example) with 2 levels (10 obs per participant per cell)
> exp1 <- gen_design(n_participant = sp$n_participant, n_stimulus = sp$n_stimulus, n_levels = sp$n_levels)
> 
> ```
> 
> **priors**
> 
> ```no-highlight
> all_priors <- c(set_prior("normal(10,1)", class = "Intercept"), 
> set_prior("normal(1,1)", class = "b"),
> set_prior("normal(5,1)", class = "sd", coef = "Intercept", group = "participant"),
> set_prior("normal(2,1)", class = "sd", coef = "x", group = "participant"),
> set_prior("normal(10,1)", class = "sigma"))
> 
> ```
> 
> **fit intial model**
> 
> ```no-highlight
> y <- rnorm(1:nrow(exp1), 1, 0) # just some arbitrary values for this column
> m1 <- brm(y ~ x + (1 + x | participant), data = exp1, sample_prior = "only", prior = all_priors)
> 
> ```
> 
> **make and fit predictions**
> 
> ```no-highlight
> new_y_sum <- predict(m1)[,1]
> new_y<- predict(m1, newdata = exp1, summary = F, nsamples = 1)[1,]
> exp1$new_y_sum <- new_y_sum
> exp1$new_y <- new_y
> m1_new_sum <- brm(new_y_sum ~ x + (1 + x | participant), data = exp1, chains = 4, cores = 4)
> m1_new <- brm(new_y ~ x + (1 + x | participant), data = exp1, chains = 4, cores = 4)
> 
> ```
> 
> **model output for predict(m1) - too small SEs**
> 
> ```no-highlight
> fixef(m1_new_sum)
> 
> ```
> 
> ```no-highlight
> Estimate Est.Error Q2.5 Q97.5
> Intercept 9.488095 0.01631371 9.4562210 9.520331
> x 1.009386 0.01550281 0.9784959 1.039939
> 
> ```
> 
> **model output for predict(m1, summary = F) intended SEs**
> 
> ```no-highlight
> fixef(m1_new)
> 
> ```
> 
> ```no-highlight
> Estimate Est.Error Q2.5 Q97.5
> Intercept 10.0320100 1.1200567 7.892897 12.145149
> x 0.3468883 0.8816665 -1.388313 2.085174
> 
> ```
> 
> > **More formal simualtion of SEs for predict(m1, summary = F)**
> >
> > ```no-highlight
> > set.seed(1)
> > fixef_vec_int <- list()
> > fixef_vec_x <- list()
> > fixef_vec_x_lb <- c()
> > 
> > for(i in 1:1000){
> > new_y <- predict(m1, newdata = exp1, summary = F, nsamples = 1)[1,]
> > exp1$new_y <- new_y
> > m1_tmp <- update(m1_new, newdata = exp1, chains = 1)
> > fixef_vec_int[[i]] <- c(fixef(m1_tmp)[1], fixef(m1_tmp)[3])
> > fixef_vec_x[[i]] <- c(fixef(m1_tmp)[2], fixef(m1_tmp)[4])
> > fixef_vec_x_lb[i] <- fixef(m1_tmp)[6]
> > }
> > 
> > fixef_x_means <- unlist(lapply(1:100, function(x) fixef_vec_x[[x]][1]))
> > fixef_x_sds <- unlist(lapply(1:100, function(x) fixef_vec_x[[x]][2]))
> > 
> > mean(fixef_x_means) # the mean of the coef-prior
> > mean(fixef_x_sds) # the mean of the sampled SE
> > 
> > ```

> [@Jacquie](#):
>
> What I would try to do coding manually would be to first draw regression weights for each person from the normal distribution across people and then create predictions for each measurement per person.

That sounds correct to me as well. [Schad, Betancourt, & Vasishth](https://osf.io/xs4zg/) have a good paper + script for these kind of simulations for mixed models.

Hope this helps,  
Julian

---

<div class="post-metadata">

**Author:** ![franzsf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/franzsf/32/8545_2.png) [@franzsf](https://discourse.mc-stan.org/u/franzsf)\
**Post date:** [July 8, 2020, 2:42pm UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881/6 "2020-07-08T14:42:45Z")

</div>

I would just add that you should thoroughly investigate the `re_formula` and `sample_new_levels` arguments to `predict` to ensure you’re getting what you expect.

---

<div class="post-metadata">

**Author:** ![Jacquie](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jacquie/32/19488_2.png) [@Jacquie](https://discourse.mc-stan.org/u/Jacquie)\
**Post date:** [July 13, 2020, 7:11am UTC](https://discourse.mc-stan.org/t/simulating-data-using/5881/7 "2020-07-13T07:11:17Z")

</div>

> [@julianquandt](#):
>
> newY=predict(simulModel, newdata = df, summary=FALSE, nsamples=1)[1,]

thank you, that’s really useful!