Simulation-based calibration checking

About

This vignette shows the results of a simulation-based calibration (SBC) checking study to validate the implementation of the models in brms.mmrm. SBC checking tests the ability of a Bayesian model to recapture the parameters used to simulate prior predictive data. For details on SBC checking, please read Modrák et al. (2024) and the SBC R package (Kim et al. 2022). This particular SBC checking study uses the targets pipeline in the viggnettes/sbc/ subdirectory of the brms.mmrm package source code.

Conclusion

From the results below, the SBC rank statistics are approximately uniformly distributed. In other words, the posterior distribution from the brms/Stan MMRM modeling code matches the prior from which the datasets were simulated. This is evidence that both the subgroup and non-subgroup models in brms.mmrm are implemented correctly.

Setup

To show the SBC checking results in this vignette, we first load code from the SBC checking study, and we use custom functions to read and plot rank statistics:

library(dplyr)
library(ggplot2)
library(tibble)
library(tidyr)

source("sbc/R/prior.R")
source("sbc/R/response.R")
source("sbc/R/scenarios.R")

read_ranks <- function(path) {
  fst::read_fst(path) |>
    tibble::as_tibble() |>
    pivot_longer(
      cols = everything(),
      names_to = "parameter",
      values_to = "rank"
    )
}

plot_ranks <- function(ranks) {
  ggplot(ranks) +
    geom_histogram(
      aes(x = rank),
      breaks = seq(from = 0, to = max(ranks$rank), length.out = 10)
    ) +
    facet_wrap(~parameter)
}

Each section below is its own SBC checking study based on a given modeling scenario. Each scenario shows results from 1000 independent simulations from the prior.

Subgroup scenario

The subgroup scenario distinguishes itself from the others by the presence of a subgroup factor. Assumptions:

  • 2 treatment groups
  • 2 subgroup levels
  • 3 time points
  • 150 patients per treatment group
  • Baseline covariates (2 continuous and 2 categorical)
  • 30% dropout
  • 8% rate of independent/sporadic missing response.

Model formula:

#> response ~ group + group:subgroup + group:subgroup:time + group:time + subgroup + subgroup:time + time + continuous1 + continuous2 + balanced + unbalanced + unstr(time = time, gr = patient) 
#> sigma ~ 0 + time

The prior was randomly generated and used for both simulation and analysis:

setup_prior(subgroup) |>
  select(prior, class, coef, dpar) |>
  as.data.frame()
#>                      prior     class                                       coef
#> 1   normal(0.0915, 2.1422) Intercept                                           
#> 3   normal(0.1261, 2.2076)         b                             balancedlevel2
#> 4    normal(0.246, 2.5987)         b                             balancedlevel3
#> 5   normal(-0.0079, 0.543)         b                                continuous1
#> 6  normal(-0.0015, 2.5463)         b                                continuous2
#> 7   normal(0.1521, 2.3486)         b                               groupgroup_2
#> 8   normal(0.2304, 2.2328)         b            groupgroup_2:subgroupsubgroup_2
#> 9   normal(0.2223, 1.3432)         b groupgroup_2:subgroupsubgroup_2:timetime_2
#> 10 normal(-0.0966, 1.6992)         b groupgroup_2:subgroupsubgroup_2:timetime_3
#> 11  normal(-0.1506, 1.602)         b                    groupgroup_2:timetime_2
#> 12  normal(0.1017, 1.0989)         b                    groupgroup_2:timetime_3
#> 13  normal(0.0671, 2.0957)         b                         subgroupsubgroup_2
#> 14  normal(0.0501, 1.1963)         b              subgroupsubgroup_2:timetime_2
#> 15 normal(-0.1145, 1.5875)         b              subgroupsubgroup_2:timetime_3
#> 16  normal(0.0792, 0.8305)         b                                 timetime_2
#> 17 normal(-0.0788, 0.5149)         b                                 timetime_3
#> 18  normal(0.1882, 2.3792)         b                           unbalancedlevel2
#> 19  normal(0.0803, 2.7501)         b                           unbalancedlevel3
#> 20             lkj(1.3486)   cortime                                           
#> 22  normal(-0.0487, 2.384)         b                                 timetime_1
#> 23 normal(-0.0795, 0.6806)         b                                 timetime_2
#> 24 normal(-0.0439, 2.1926)         b                                 timetime_3
#>     dpar
#> 1       
#> 3       
#> 4       
#> 5       
#> 6       
#> 7       
#> 8       
#> 9       
#> 10      
#> 11      
#> 12      
#> 13      
#> 14      
#> 15      
#> 16      
#> 17      
#> 18      
#> 19      
#> 20      
#> 22 sigma
#> 23 sigma
#> 24 sigma

The following histograms show the SBC rank statistics which compare the prior parameter draws draws to the posterior draws. If the data simulation code and modeling code are both correct and consistent, then the rank statistics should be approximately uniformly distributed.

ranks_subgroup <- read_ranks("sbc/results/subgroup.fst")

Fixed effect parameter ranks:

ranks_subgroup |>
  filter(grepl("^b_", parameter)) |>
  filter(!grepl("^b_sigma", parameter)) |>
  plot_ranks()

Standard deviation parameter ranks:

ranks_subgroup |>
  filter(grepl("b_sigma", parameter)) |>
  plot_ranks()

Correlation parameter ranks:

ranks_subgroup |>
  filter(grepl("cortime_", parameter)) |>
  plot_ranks()

Unstructured scenario

This scenario uses unstructured correlation and does not use a subgroup variable. Assumptions:

  • 3 treatment groups
  • No subgroup
  • 4 time points
  • 100 patients per treatment group
  • No covariate adjustment
  • No missing responses

Model formula:

#> response ~ 0 + group + time + unstr(time = time, gr = patient) 
#> sigma ~ 0 + time

The prior was randomly generated and used for both simulation and analysis:

setup_prior(unstructured) |>
  select(prior, class, coef, dpar) |>
  as.data.frame()
#>                      prior   class         coef  dpar
#> 2    normal(0.105, 0.4728)       b groupgroup_1      
#> 3  normal(-0.1704, 2.5406)       b groupgroup_2      
#> 4   normal(0.2469, 1.2809)       b groupgroup_3      
#> 5  normal(-0.1022, 0.2827)       b   timetime_2      
#> 6   normal(0.1312, 2.0238)       b   timetime_3      
#> 7  normal(-0.1573, 1.1288)       b   timetime_4      
#> 8              lkj(1.1293) cortime                   
#> 10  normal(0.1593, 1.4516)       b   timetime_1 sigma
#> 11  normal(-0.1321, 1.431)       b   timetime_2 sigma
#> 12 normal(-0.1661, 1.9156)       b   timetime_3 sigma
#> 13   normal(0.032, 0.5265)       b   timetime_4 sigma

SBC checking rank statistics:

ranks_unstructured <- read_ranks("sbc/results/unstructured.fst")

Fixed effect parameter ranks:

ranks_unstructured |>
  filter(grepl("^b_", parameter)) |>
  filter(!grepl("^b_sigma", parameter)) |>
  plot_ranks()

Log-scale standard deviation parameter ranks:

ranks_unstructured |>
  filter(grepl("b_sigma", parameter)) |>
  plot_ranks()

Correlation parameter ranks:

ranks_unstructured |>
  filter(grepl("cortime_", parameter)) |>
  plot_ranks()

Autoregressive moving average scenario

This scenario uses an autoregressive moving average (ARMA) model with autoregressive order 1 and moving average order 1. Assumptions:

  • 2 treatment groups
  • No subgroup
  • 3 time points
  • 100 patients per treatment group
  • No covariate adjustment
  • No missing responses

Model formula:

#> response ~ 0 + group + time + arma(time = time, gr = patient, p = 1L, q = 1L, cov = FALSE) 
#> sigma ~ 0 + time

The prior was randomly generated and used for both simulation and analysis:

setup_prior(autoregressive_moving_average) |>
  select(prior, class, coef, dpar) |>
  as.data.frame()
#>                      prior class         coef  dpar
#> 1        uniform(0.1, 0.9)    ar                   
#> 3  normal(-0.1382, 2.0354)     b groupgroup_1      
#> 4  normal(-0.0841, 0.5827)     b groupgroup_2      
#> 5    normal(-0.208, 2.734)     b   timetime_2      
#> 6   normal(0.0233, 0.3739)     b   timetime_3      
#> 7        uniform(0.1, 0.9)    ma                   
#> 9  normal(-0.0289, 2.5524)     b   timetime_1 sigma
#> 10   normal(0.1457, 0.457)     b   timetime_2 sigma
#> 11 normal(-0.0895, 2.3152)     b   timetime_3 sigma

SBC checking rank statistics:

ranks_autoregressive_moving_average <- read_ranks(
  "sbc/results/autoregressive_moving_average.fst"
)

Fixed effect parameter ranks:

ranks_autoregressive_moving_average |>
  filter(grepl("^b_", parameter)) |>
  filter(!grepl("^b_sigma", parameter)) |>
  plot_ranks()

Log-scale standard deviation parameter ranks:

ranks_autoregressive_moving_average |>
  filter(grepl("b_sigma", parameter)) |>
  plot_ranks()

Correlation parameter ranks:

ranks_autoregressive_moving_average |>
  filter(parameter %in% c("ar[1]", "ma[1]")) |>
  plot_ranks()

Autoregressive scenario

This scenario is the same as above, but the correlation structure is autoregressive with order 2. Model formula:

#> response ~ 0 + group + time + ar(time = time, gr = patient, p = 2L, cov = FALSE) 
#> sigma ~ 0 + time

The prior was randomly generated and used for both simulation and analysis:

setup_prior(autoregressive) |>
  select(prior, class, coef, dpar) |>
  as.data.frame()
#>                      prior class         coef  dpar
#> 1        uniform(0.1, 0.9)    ar                   
#> 3  normal(-0.1453, 1.2003)     b groupgroup_1      
#> 4   normal(0.0259, 2.2736)     b groupgroup_2      
#> 5   normal(0.0017, 1.6365)     b   timetime_2      
#> 6   normal(0.0326, 2.3636)     b   timetime_3      
#> 8   normal(0.1486, 2.3464)     b   timetime_1 sigma
#> 9   normal(0.0837, 0.9759)     b   timetime_2 sigma
#> 10  normal(0.0207, 1.2563)     b   timetime_3 sigma

SBC checking rank statistics:

ranks_autoregressive <- read_ranks("sbc/results/autoregressive.fst")

Fixed effect parameter ranks:

ranks_autoregressive |>
  filter(grepl("^b_", parameter)) |>
  filter(!grepl("^b_sigma", parameter)) |>
  plot_ranks()

Log-scale standard deviation parameter ranks:

ranks_autoregressive |>
  filter(grepl("b_sigma", parameter)) |>
  plot_ranks()

Correlation parameter ranks:

ranks_autoregressive |>
  filter(parameter %in% c("ar[1]", "ar[2]")) |>
  plot_ranks()

Moving average scenario

This scenario is the same as above, but it uses a moving average correlation structure with order 2. Model formula:

#> response ~ 0 + group + time + ma(time = time, gr = patient, q = 2L, cov = FALSE) 
#> sigma ~ 0 + time

The prior was randomly generated and used for both simulation and analysis:

setup_prior(moving_average) |>
  select(prior, class, coef, dpar) |>
  as.data.frame()
#>                      prior class         coef  dpar
#> 2   normal(0.1644, 1.9148)     b groupgroup_1      
#> 3   normal(0.0564, 1.6111)     b groupgroup_2      
#> 4    normal(0.0013, 2.095)     b   timetime_2      
#> 5    normal(0.1472, 0.791)     b   timetime_3      
#> 6        uniform(0.1, 0.9)    ma                   
#> 8   normal(0.0203, 0.9309)     b   timetime_1 sigma
#> 9   normal(0.2004, 1.4374)     b   timetime_2 sigma
#> 10 normal(-0.1997, 1.8687)     b   timetime_3 sigma

SBC checking rank statistics:

ranks_moving_average <- read_ranks("sbc/results/moving_average.fst")

Fixed effect parameter ranks:

ranks_moving_average |>
  filter(grepl("^b_", parameter)) |>
  filter(!grepl("^b_sigma", parameter)) |>
  plot_ranks()

Log-scale standard deviation parameter ranks:

ranks_moving_average |>
  filter(grepl("b_sigma", parameter)) |>
  plot_ranks()

Correlation parameter ranks:

ranks_moving_average |>
  filter(parameter %in% c("ma[1]", "ma[2]")) |>
  plot_ranks()

Compound symmetry scenario

This scenario is the same as above, but it uses a compound symmetry correlation structure. Model formula:

#> response ~ 0 + group + time + cosy(time = time, gr = patient) 
#> sigma ~ 0 + time

The prior was randomly generated and used for both simulation and analysis:

setup_prior(compound_symmetry) |>
  select(prior, class, coef, dpar) |>
  as.data.frame()
#>                      prior class         coef  dpar
#> 2   normal(0.0941, 1.4367)     b groupgroup_1      
#> 3   normal(0.1114, 1.1299)     b groupgroup_2      
#> 4   normal(0.2303, 0.7995)     b   timetime_2      
#> 5   normal(0.1673, 2.5556)     b   timetime_3      
#> 6        uniform(0.1, 0.9)  cosy                   
#> 8  normal(-0.1927, 1.6853)     b   timetime_1 sigma
#> 9   normal(0.1181, 0.3831)     b   timetime_2 sigma
#> 10  normal(0.2009, 0.4912)     b   timetime_3 sigma

SBC checking rank statistics:

ranks_compound_symmetry <- read_ranks("sbc/results/compound_symmetry.fst")

Fixed effect parameter ranks:

ranks_compound_symmetry |>
  filter(grepl("^b_", parameter)) |>
  filter(!grepl("^b_sigma", parameter)) |>
  plot_ranks()

Log-scale standard deviation parameter ranks:

ranks_compound_symmetry |>
  filter(grepl("b_sigma", parameter)) |>
  plot_ranks()

Correlation parameter ranks:

ranks_compound_symmetry |>
  filter(parameter == "cosy") |>
  plot_ranks()

Diagonal scenario

This scenario is the same as above, but it uses a diagonal correlation structure (independent time points within patients). Model formula:

#> response ~ 0 + group + time 
#> sigma ~ group + group:time + time

The prior was randomly generated and used for both simulation and analysis:

setup_prior(diagonal) |>
  select(prior, class, coef, dpar) |>
  as.data.frame()
#>                      prior     class                    coef  dpar
#> 2  normal(-0.0901, 1.1291)         b            groupgroup_1      
#> 3    normal(0.0094, 1.342)         b            groupgroup_2      
#> 4   normal(0.0518, 2.7578)         b              timetime_2      
#> 5    normal(0.068, 1.0882)         b              timetime_3      
#> 6    normal(0.203, 2.2932) Intercept                         sigma
#> 8    normal(-0.0675, 0.43)         b            groupgroup_2 sigma
#> 9   normal(0.1825, 2.9551)         b groupgroup_2:timetime_2 sigma
#> 10 normal(-0.0941, 2.6258)         b groupgroup_2:timetime_3 sigma
#> 11  normal(0.0708, 1.9776)         b              timetime_2 sigma
#> 12 normal(-0.0574, 1.2542)         b              timetime_3 sigma

SBC checking rank statistics:

ranks_diagonal <- read_ranks("sbc/results/diagonal.fst")

Fixed effect parameter ranks:

ranks_diagonal |>
  filter(grepl("^b_", parameter)) |>
  filter(!grepl("^b_sigma", parameter)) |>
  plot_ranks()

Log-scale standard deviation parameter ranks:

ranks_diagonal |>
  filter(grepl("b_sigma", parameter)) |>
  plot_ranks()

References

Kim, S., Moon, H., Modrák, M., and Säilynoja, T. (2022), SBC: Simulation based calibration for rstan/cmdstanr models.
Modrák, M., Moon, A. H., Kim, S., Bürkner, P., Huurre, N., Faltejsková, K., Gelman, A., and Vehtari, A. (2024), “Simulation-Based Calibration Checking for Bayesian Computation: The Choice of Test Quantities Shapes Sensitivity,” Bayesian Analysis, forthcoming. https://doi.org/10.1214/23-BA1404.