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.