Hello! I’m designing a clinical trial of an intervention to improve emotion regulation. We plan to remotely measure participants’ mood >1x per day over the course of the study.

We expect that participants’ average mood will improve in the treatment condition, but also that the *variance* in mood will decrease over time. I’m trying to simulate the experiment ahead of time so that I can make an informed decision about how many mood measurements I should be taking.

Here is code that approximates the dataset I expect to end up with. (I haven’t yet tried to identify realistic values for the parameters, just want to make sure I can recover them at this point).

```
J <- 40 # number of people in experiment
weeks <- 12 # number of weeks in treatment (participants will attend one group treatment session per week)
n_per_person_per_week <- 14 # number of EMA measurements taken between each group treatment session
data <- tibble::as_tibble(expand.grid(id = 1:J, week = 1:weeks, time = 1:n_per_person_per_week))
data <- data[order(data$id, data$week),]
data$time <- rep(1:(weeks*n_per_person_per_week), J)
data$treatment <- with(data, ifelse(id %in% 1:(J/2),
1, 0))
a <- rnorm(J, 0, 1) # varying intercept
b <- rnorm(J, .1, .2) # general effect for being in either treatment or active control is centered around d=.1
theta <- rnorm(J, .4, .25) # varying effect per person in tx group
# centered around .4
# (this effect will get multiplied by zero for people in control)
sigma_y <- 1 # noise
y_pred <- a[data$id] + # varying intercept
b[data$id]*data$week + # treatment effect (common to both tx and active control)
theta[data$id]*data$treatment*data$week # varying slope for each person in treatment group
# we expect the variance in mood to decrease over time in the treatment group
emotion_dysregulation <- data$treatment * data$week * -.05
data$y <- rnorm(nrow(data), y_pred,
sigma_y + emotion_dysregulation)
data %>%
mutate(treatment = factor(treatment, levels=c(1,0), labels=c("treatment", "active control"))) %>%
ggplot() +
aes(x = time, y = y, group = id, color=treatment) +
geom_line(alpha=.5) +
ylab("mood") +
guides(color=FALSE) +
facet_wrap(~treatment)
```

The above code will produce data that looks something like this:

Now I use brms to fit a multilevel model of overall changes in mood, while also checking if the variance in mood decreases in the treatment group

```
fit <- brm(bf(
y ~ week + week:treatment +
(1 + week + week:treatment | id), # varying intercepts and slopes
sigma ~ week + week:treatment),
data=data,
cores=4)
```

The model does a good job recovering the parameters for the linear model. But look at the estimates related to sigma:

```
Estimate Est.Error Q2.5 Q97.5
sigma_Intercept 0.04 0.02 0.00 0.07
sigma_week 0.00 0.00 -0.01 0.00
sigma_week:treatment -0.07 0.00 -0.07 -0.06
```

The estimate for the effect of week (true value 0) and the week:treatment interaction (true value -.05) are okay, but the estimate for sigma_Intercept seems to imply essentially zero noise in the model, which is wrong and isn’t reflective of the plots derived from this model using marginal_effects().

So, I’m wondering if there is a problem with how the value of sigma_Intercept is reported for functions like print() or fixef(). Or, maybe I’m just totally misinterpreting something and there is no problem! I am prepared to be embarrassed and hopefully educated!

Thank you very much for your time.