Posterior predictive check for Bayesian mixed effects model for repeated measures

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

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)

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),
pp_check(stanreg1, plotfun = “violin_grouped”, group = visit , y_draw=“points”)

Oh, and here is my brms implementation.

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”)))

This sounds like a question for @jonah.

Currently pp_check is only going to let you group by one variable. But allowing two grouping variables for plots is on the to-do list. For now I think you’ll have to manually create the plot you want using ggplot2.

1 Like

A quick hack: Change the grouping variable to include information from both grouping variables.

pp_check(stanreg1, plotfun = "violin_grouped", group = interaction(visit, treatment), y_draw="points")

1 Like

Nice answer. One can also change around the order:
pp_check(stanreg1, plotfun = “violin_grouped”, group = factor( interaction(visit, treatment), levels = c(“1.0”,“1.1”,“2.0”, “2.1”, “3.0”, “3.1”, “4.0”, “4.1”, “5.0”, “5.1”)), y_draw=“points”)
But I did not manage to change to colors, but this is pretty nice already.