I was looking into fitting a Bayesian mixed effects model for repeated measures (MMRM). The standard situation for this model would be that we have a clinical trial with multiple independent patients that are assigned into two treatment groups and the same variable is measured for each patient at multiple visits that occur at pre-specified times for the same patient (each time under the same treatment). As a result the observations within a patient are correlated, which in a standard frequentist analayis one would capture with an unstructured covariance matrix (either the same one for both treatment groups, or separate ones for each treatment group). I am attaching an plot of simulated data (see code below) illustrating what the data may look like.

Any comments on the model implementation in my rstanarm code below are of course welcome.

What I am trying to figure out is how I can get the posterior predictive plot to be grouped by both visit and treatment group.

patient <- rep(1:20, each=5)

visit <- rep(1:5, 20)## baseline (pre-treatment) value for each patient

base <- rep( rnorm(20, 0, 1), each=5)

epsilon <- rnorm(100, 0, 1)

treatment <- rep(0:1, each=50)## value at each visit

aval <- base + visit*treatment + treatment + epsilon

## change grom baseline at each visit

chg <- aval - base

testdata <- data.frame(patient=as.factor(patient), treatment=as.factor(treatment), visit=as.factor(visit), chg=chg, aval=aval, base=base)## Illustrate what the data look like

library(ggplot2)

ggplot(testdata, aes(visit, chg, col=treatment)) + theme_grey(base_size = 24) + scale_colour_manual(values=setNames(c(rgb(0,0,1,0.5),rgb(1,0,0,0.5)), c(0,1))) + geom_jitter(width=0.1, size=3)

library(rstanarm)

stanreg1 <- stan_glmer(chg ~ 0 + visit + (visit|patient) + base + treatment + visitbase + visittreatment,

data=testdata, family=â€śgaussianâ€ť,

prior_intercept = student_t(df=4, 0, 1000),

prior=student_t(df=4, 0, 10),

prior_aux = student_t(4, 0, 10),

prior_covariance = decov(regularization=2, concentration=1, shape=1, scale=1),

adapt_delta=0.99)

#summary(stanreg1)

#pp_check(stanreg1,lwd=3)

pp_check(stanreg1, plotfun = â€śviolin_groupedâ€ť, group = visit , y_draw=â€śpointsâ€ť)

Oh, and here is my brms implementation.

library(brms)

bsfit2 <- brm(chg ~ 0 + visit + (visit|patient) + base + treatment + visitbase + visittreatment,

data=testdata, family=â€śgaussianâ€ť,

prior=c( prior(â€ślkj(2)â€ť, class = â€ścorâ€ť),

prior(â€śstudent_t(4,0,10)â€ť, class = â€śbâ€ť),

prior(â€śstudent_t(4,0,10)â€ť, class = â€śsdâ€ť)))