sbc

vignettes/sbc.Rmd

@@ -46,6 +46,7 @@ library(tibble)
 library(tidyr)
 
 source("sbc/R/prior.R")
+source("sbc/R/response.R")
 source("sbc/R/scenarios.R")
 
 read_ranks <- function(path) {
@@ -152,7 +153,7 @@ setup_prior(unstructured) |>
   as.data.frame()
 ```
 
-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 uniformly distributed.
+SBC checking rank statistics:
 
 ```{r}
 ranks_unstructured <- read_ranks("sbc/results/unstructured.fst")
@@ -182,9 +183,10 @@ ranks_unstructured |>
   plot_ranks()
 ```
 
-# Autoregressive scenario
 
-This scenario uses an autoregressive model with order 1 for the correlation structure. Assumptions:
+# 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
@@ -195,6 +197,55 @@ This scenario uses an autoregressive model with order 1 for the correlation stru
 
 Model formula:
 
+```{r, echo = FALSE}
+autoregressive_moving_average()$formula
+```
+
+The prior was randomly generated and used for both simulation and analysis:
+
+```{r, paged.print = FALSE}
+setup_prior(autoregressive_moving_average) |>
+  select(prior, class, coef, dpar) |>
+  as.data.frame()
+```
+
+SBC checking rank statistics:
+
+```{r}
+ranks_autoregressive_moving_average <- read_ranks(
+  "sbc/results/autoregressive_moving_average.fst"
+)
+```
+
+Fixed effect parameter ranks:
+
+```{r}
+ranks_autoregressive_moving_average |>
+  filter(grepl("^b_", parameter)) |>
+  filter(!grepl("^b_sigma", parameter)) |>
+  plot_ranks()
+```
+
+Log-scale standard deviation parameter ranks:
+
+```{r}
+ranks_autoregressive_moving_average |>
+  filter(grepl("b_sigma", parameter)) |>
+  plot_ranks()
+```
+
+Correlation parameter ranks:
+
+```{r}
+ranks_autoregressive_moving_average |>
+  filter(parameter %in% c("ar[1]", "ar[2]", "ma[1]", "ma[2]")) |>
+  plot_ranks()
+```
+
+# Autoregressive scenario
+
+This scenario is the same as above, but the correlation structure is AR(1). Model formula:
+
 ```{r, echo = FALSE}
 autoregressive()$formula
 ```
@@ -207,7 +258,7 @@ setup_prior(autoregressive) |>
   as.data.frame()
 ```
 
-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 uniformly distributed.
+SBC checking rank statistics:
 
 ```{r}
 ranks_autoregressive <- read_ranks("sbc/results/autoregressive.fst")
@@ -237,18 +288,55 @@ ranks_autoregressive |>
   plot_ranks()
 ```
 
-# Compound symmetry scenario
+# Moving average scenario
 
-This scenario uses compound symmetry for the correlation structure. Assumptions:
+This scenario is the same as above, but it uses a moving average correlation structure with order 2. Model formula:
 
-* 2 treatment groups
-* No subgroup
-* 4 time points
-* 50 patients per treatment group
-* No covariate adjustment
-* No missing responses
+```{r, echo = FALSE}
+moving_average()$formula
+```
 
-Model formula:
+The prior was randomly generated and used for both simulation and analysis:
+
+```{r, paged.print = FALSE}
+setup_prior(moving_average) |>
+  select(prior, class, coef, dpar) |>
+  as.data.frame()
+```
+
+SBC checking rank statistics:
+
+```{r}
+ranks_moving_average <- read_ranks("sbc/results/moving_average.fst")
+```
+
+Fixed effect parameter ranks:
+
+```{r}
+ranks_moving_average |>
+  filter(grepl("^b_", parameter)) |>
+  filter(!grepl("^b_sigma", parameter)) |>
+  plot_ranks()
+```
+Log-scale standard deviation parameter ranks:
+
+```{r}
+ranks_moving_average |>
+  filter(grepl("b_sigma", parameter)) |>
+  plot_ranks()
+```
+
+Correlation parameter ranks:
+
+```{r}
+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:
 
 ```{r, echo = FALSE}
 compound_symmetry()$formula
@@ -262,7 +350,7 @@ setup_prior(compound_symmetry) |>
   as.data.frame()
 ```
 
-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 uniformly distributed.
+SBC checking rank statistics:
 
 ```{r}
 ranks_compound_symmetry <- read_ranks("sbc/results/compound_symmetry.fst")
@@ -294,16 +382,7 @@ ranks_compound_symmetry |>
 
 # Diagonal scenario
 
-This scenario assumes time points within patients are independent (diagonal correlation matrix). Other assumptions:
-
-* 2 treatment groups
-* No subgroup
-* 3 time points
-* 50 patients per treatment group
-* No covariate adjustment
-* No missing responses
-
-In addition, the standard deviations are parameterized differently. Model formula:
+This scenario is the same as above, but it uses a diagonal correlation structure (independent time points within patients). Model formula:
 
 ```{r, echo = FALSE}
 diagonal()$formula
@@ -317,7 +396,7 @@ setup_prior(diagonal) |>
   as.data.frame()
 ```
 
-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 uniformly distributed.
+SBC checking rank statistics:
 
 ```{r}
 ranks_diagonal <- read_ranks("sbc/results/diagonal.fst")

vignettes/sbc/R/response.R

@@ -33,6 +33,80 @@ simulate_unstructured <- function(data, formula, prior) {
   list(data = data, parameters = parameters)
 }
 
+simulate_arma22 <- function(data, formula, prior) {
+  beta <- simulate_beta(data = data, formula = formula, prior = prior)
+  x_beta <- derive_x_beta(
+    data = data,
+    formula = formula,
+    prior = prior,
+    beta = beta
+  )
+  b_sigma <- simulate_b_sigma(data = data, formula = formula, prior = prior)
+  sigma <- derive_sigma(
+    data = data,
+    formula = formula,
+    prior = prior,
+    b_sigma = b_sigma
+  )
+  ar <- replicate(2L, eval(parse(text = prior[prior$class == "ar", "r"])))
+  ma <- replicate(2L, eval(parse(text = prior[prior$class == "ma", "r"])))
+  residuals <- rnorm(n = length(sigma), mean = 0, sd = sigma)
+  n_time <- length(unique(data[[attr(data, "brm_time")]]))
+  n_patient <- nrow(data) / n_time
+  for (patient in seq_len(n_patient)) {
+    rows <- seq_len(n_time) + n_time * (patient - 1L)
+    e <- arima.sim(
+      n = n_time,
+      model = list(ar = ar, ma = ma),
+      sd = sigma[rows]
+    )
+    data[[attr(data, "brm_outcome")]][rows] <- x_beta[rows] + as.numeric(e)
+  }
+  data$response[data$missing] <- NA_real_
+  names(beta) <- paste0("b_", names(beta))
+  names(b_sigma) <- paste0("b_sigma_", names(b_sigma))
+  parameters <- c(
+    beta,
+    b_sigma,
+    `ar[1]` = ar[1L],
+    `ar[2]` = ar[2L],
+    `ma[1]` = ma[1L],
+    `ma[2]` = ma[2L]
+  )
+  list(data = data, parameters = parameters)
+}
+
+simulate_ma2 <- function(data, formula, prior) {
+  beta <- simulate_beta(data = data, formula = formula, prior = prior)
+  x_beta <- derive_x_beta(
+    data = data,
+    formula = formula,
+    prior = prior,
+    beta = beta
+  )
+  b_sigma <- simulate_b_sigma(data = data, formula = formula, prior = prior)
+  sigma <- derive_sigma(
+    data = data,
+    formula = formula,
+    prior = prior,
+    b_sigma = b_sigma
+  )
+  ma <- replicate(2L, eval(parse(text = prior[prior$class == "ma", "r"])))
+  residuals <- rnorm(n = length(sigma), mean = 0, sd = sigma)
+  n_time <- length(unique(data[[attr(data, "brm_time")]]))
+  n_patient <- nrow(data) / n_time
+  for (patient in seq_len(n_patient)) {
+    rows <- seq_len(n_time) + n_time * (patient - 1L)
+    e <- arima.sim(n = n_time, model = list(ma = ma), sd = sigma[rows])
+    data[[attr(data, "brm_outcome")]][rows] <- x_beta[rows] + as.numeric(e)
+  }
+  data$response[data$missing] <- NA_real_
+  names(beta) <- paste0("b_", names(beta))
+  names(b_sigma) <- paste0("b_sigma_", names(b_sigma))
+  parameters <- c(beta, b_sigma, `ma[1]` = ma[1L], `ma[2]` = ma[2L])
+  list(data = data, parameters = parameters)
+}
+
 simulate_ar1 <- function(data, formula, prior) {
   beta <- simulate_beta(data = data, formula = formula, prior = prior)
   x_beta <- derive_x_beta(

vignettes/sbc/R/scenarios.R

@@ -65,6 +65,61 @@ unstructured <- function() {
   list(data = data, formula = formula, simulate = simulate_unstructured)
 }
 
+autoregressive_moving_average <- function() {
+  n_group <- 2L
+  n_patient <- 50L
+  n_time <- 4L
+  data <- brms.mmrm::brm_simulate_outline(
+    n_group = n_group,
+    n_patient = n_patient,
+    n_time = n_time,
+    rate_dropout = 0,
+    rate_lapse = 0
+  )
+  data[[attr(data, "brm_outcome")]] <- seq_len(nrow(data))
+  formula <- brms.mmrm::brm_formula(
+    data = data,
+    intercept = FALSE,
+    baseline = FALSE,
+    group = TRUE,
+    time = TRUE,
+    baseline_time = FALSE,
+    group_time = FALSE,
+    covariates = FALSE,
+    correlation = "autoregressive_moving_average",
+    autoregressive_order = 2L,
+    moving_average_order = 2L
+  )
+  list(data = data, formula = formula, simulate = simulate_arma22)
+}
+
+moving_average <- function() {
+  n_group <- 2L
+  n_patient <- 50L
+  n_time <- 4L
+  data <- brms.mmrm::brm_simulate_outline(
+    n_group = n_group,
+    n_patient = n_patient,
+    n_time = n_time,
+    rate_dropout = 0,
+    rate_lapse = 0
+  )
+  data[[attr(data, "brm_outcome")]] <- seq_len(nrow(data))
+  formula <- brms.mmrm::brm_formula(
+    data = data,
+    intercept = FALSE,
+    baseline = FALSE,
+    group = TRUE,
+    time = TRUE,
+    baseline_time = FALSE,
+    group_time = FALSE,
+    covariates = FALSE,
+    correlation = "moving_average",
+    moving_average_order = 2L
+  )
+  list(data = data, formula = formula, simulate = simulate_ma2)
+}
+
 autoregressive <- function() {
   n_group <- 2L
   n_patient <- 50L