Challenging dataset (for the frequentist mmrm package)

[wlandau]

It would be interesting to analyze the dataset from openpharma/mmrm#380 using brms.mmrm. Depending on the results, it might make a useful vignette. Thanks @yonicd for the idea.

[wlandau]

Predictably, this is a computationally demanding dataset because of the number of data points, number of time points, and heteroskedasticity mentioned in the original issue thread. I ran brms.mmrm (branch 102) on the data for 2000 iterations per chain, and it did not even come close to converging. I won't have time to come back to this for a while, but here is the code I used (mostly borrowed from the thread). Next step would be to try with more iterations and warmup.

library(brms.mmrm)
library(MASS)
library(mvtnorm)
library(mmrm)
library(dplyr)
library(tidyr)
library(ggplot2)
library(emmeans)
library(broom)
library(tictoc)

sim_data_fn = function(
    N=10, 
    n=200,
    mu=-100/52,
    delta = 50/52 ,
    mua = 2000,
    sigmaa = 300, 
    sigmab = 60, 
    corab = 0.2, 
    sigma = 10,
    times = c(0,2,6,12,24,36,52,70,88,104)
){
  nt = length(times)                                          # no of repeated measures
  out = as.data.frame(list(1,1,"Tmp",1,1))
  names(out) = c("Sim","Pts","Trt","Time","FEV1")                 # out is the output
  
  outi = as.data.frame(list(rep(0,n*nt),rep(0,n*nt),rep("Trt",n*nt),rep(0,n*nt),rep(0,n*nt)))
  names(outi) = names(out)                                        # outi is the output for each replication
  
  outi[,"Pts"] = rep(1:n,each=nt)                                 
  outi[,"Trt"] = c(rep("Placebo",n*nt/2),rep("Active",n*nt/2))    
  outi[,"Time"] = rep(times,n)                                   # the time of each FEV1 measurement
  covab = corab * sigmaa * sigmab                                # cov between a and b
  COV = matrix(c(sigmaa^2,covab,covab,sigmab^2),ncol=2)          # Cov matrix for the slope and intercept 
  
  Pts.ID = 1:n
  Trt = c(rep("Placebo",n/2),rep("Active",n/2))
  
  for (i in 1:N)                                                 # Loop over the N simulated studies
  {
    outi[,"Sim"] = rep(i,n*nt)
    si_list <- vector( length = n, mode = 'list') 
    for (j in 1:n)                                            # Loop over the n patients
    {
      si_list[[j]] = rmvnorm(1, mean = c(mua, mu + delta*(Trt[j]=="Active")),sigma=COV)   
      outi[((j-1)*nt+1):(j*nt),"FEV1"] = si_list[[j]][1] + si_list[[j]][2] *times + rnorm(nt,0,sd=sigma)
    }
    out = rbind(out,outi)                                      # Append study to the full data
  }
  
  out = out[-1,]                                                # remove temporary row
  return(out)
}     

set.seed(123)
sim_data <- sim_data_fn() |> tibble::as_tibble()

cov_params <- list(
  sigmaa = 300, 
  sigmab = 60, 
  corab = 0.2
)

si_params <- list(
  mu=-100/52,
  mua = 2000,
  delta = 50/52
)

set.seed(123)
sim_data2 <- tidyr::expand_grid(
  N = 10,
  n = 200,
  Sim = seq(N),
  Pts = seq(n),
  Time = c(0,2,6,12,24,36,52,70,88,104)
) |>
  dplyr::mutate(
    Trt = ifelse(Pts<=n/2,'Placebo','Active'),
    Trt = factor(Trt, levels = c('Placebo','Active'))
  ) |>
  dplyr::mutate(
    cov_params = list(cov_params),
    si_params = list(si_params),
    sigma = 10,
    # Define covariance matrix for the sim iteration
    cov_mat = purrr::map(cov_params,function(cov_params){
      list2env(cov_params, envir = environment())
      # cov between a and b
      covab = corab * sigmaa * sigmab
      # Cov matrix for the slope and intercept 
      matrix(c(sigmaa^2,covab,covab,sigmab^2),ncol=2)
    })
  ) |>
  tidyr::nest(data = -c(Sim, Trt, cov_mat, sigma, Pts, si_params)) |>
  dplyr::mutate(
    # draw random intercept/slope for each treatment arm
    si = purrr::pmap(list(trtj = Trt, cov_mat = cov_mat, si_params = si_params), function(trtj, cov_mat, si_params){
      list2env(si_params, envir = environment())
      rmvnorm(1, mean = c(mua, mu + delta*(trtj=="Active")), sigma = cov_mat)
    })
  ) |>
  tidyr::unnest(c(data)) |>
  dplyr::mutate(
    # generate fev outcomes per subject
    FEV1 = purrr::pmap_dbl(list(time_j = Time, si = si, sigma = sigma), function(time_j, si, sigma){
      si[1] + si[2] * time_j + rnorm(1, 0, sd = sigma)
    })
  ) |>
  mutate(
    Time = as.factor(Time),
    patient = rep(replicate(2000, basename(tempfile("x"))), each = 10)
  )

data <- brm_data(
  data = sim_data2,
  outcome = "FEV1",
  role = "response",
  group = "Trt",
  time = "Time",
  patient = "patient",
  reference_group = "Placebo",
  reference_time = "0"
)

formula <- brm_formula(
  data = data,
  intercept = TRUE,
  group = TRUE,
  time = TRUE,
  group_time = TRUE,
  correlation = "autoregressive",
  autoregressive_order = 1L,
  sigma = brm_formula_sigma(data = data, intercept = FALSE, time = TRUE)
)

fit <- brm_model(
  data = data,
  formula = formula,
  cores = 4,
  refresh = 10
)