R/Simulations-class.R

In crmPack: Object-Oriented Implementation of Dose Escalation Designs

#' @include helpers.R
#' @include Data-class.R
#' @include Simulations-validity.R
#' @include CrmPackClass-class.R
NULL

# GeneralSimulations ----

## class ----

#' `GeneralSimulations`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This class captures trial simulations.
#' Here also the random generator state before starting the simulation is
#' saved, in order to be able to reproduce the outcome. For this just use
#' [`set.seed`] with the `seed` as argument before running
#' [`simulate,Design-method`].
#'
#' @slot data (`list`)\cr produced [`Data`] objects.
#' @slot doses (`numeric`)\cr final dose recommendations.
#' @slot seed (`integer`)\cr random generator state before starting the simulation.
#'
#' @aliases GeneralSimulations
#' @export
.GeneralSimulations <-
  setClass(
    Class = "GeneralSimulations",
    slots = c(
      data = "list",
      doses = "numeric",
      seed = "integer"
    ),
    prototype = prototype(
      data = list(
        Data(
          x = 1:2,
          y = 0:1,
          doseGrid = 1:2,
          ID = 1L:2L,
          cohort = 1L:2L
        ),
        Data(
          x = 3:4,
          y = 0:1,
          doseGrid = 3:4,
          ID = 1L:2L,
          cohort = 1L:2L
        )
      ),
      doses = c(1, 2),
      seed = 1L
    ),
    contains = "CrmPackClass",
    validity = v_general_simulations
  )

## constructor ----

#' @rdname GeneralSimulations-class
#'
#' @param data (`list`)\cr see slot definition.
#' @param doses (`numeric`)\cr see slot definition.
#' @param seed (`integer`)\cr see slot definition.
#'
#' @example examples/Simulations-class-GeneralSimulations.R
#' @export
GeneralSimulations <- function(data, doses, seed) {
  assert_integerish(seed)
  .GeneralSimulations(
    data = data,
    doses = doses,
    seed = as.integer(seed)
  )
}


## default constructor

#' @rdname GeneralSimulations-class
#' @note Typically, end users will not use the `.DefaultGeneralSimulations()` function.
#' @export
.DefaultGeneralSimulations <- function() {
  GeneralSimulations(
    data = list(
      Data(x = 1:3, y = c(0, 1, 0), doseGrid = 1:3, ID = 1L:3L, cohort = 1L:3L),
      Data(x = 4:6, y = c(0, 1, 0), doseGrid = 4:6, ID = 1L:3L, cohort = 1L:3L)
    ),
    doses = c(1, 2),
    seed = 123
  )
}


# Simulations ----

## class ----

#' `Simulations`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This class captures the trial simulations from model based designs.
#' Additional slots `fit`, `stop_reasons`, `stop_report`,`additional_stats` compared to
#' the general class [`GeneralSimulations`].
#'
#' @slot fit (`list`)\cr final fits
#' @slot stop_reasons (`list`)\cr stopping reasons for each simulation run
#' @slot stop_report matrix of stopping rule outcomes
#' @slot additional_stats list of additional statistical summary
#' @aliases Simulations
#' @export
.Simulations <-
  setClass(
    Class = "Simulations",
    slots = c(
      fit = "list",
      stop_report = "matrix",
      stop_reasons = "list",
      additional_stats = "list"
    ),
    prototype = prototype(
      fit = list(
        c(0.1, 0.2),
        c(0.1, 0.2)
      ),
      stop_report = matrix(TRUE, nrow = 2),
      stop_reasons = list("A", "A"),
      additional_stats = list(a = 1, b = 1)
    ),
    contains = "GeneralSimulations",
    validity = v_simulations
  )

## constructor ----

#' @rdname Simulations-class
#'
#' @param fit (`list`)\cr see slot definition.
#' @param stop_reasons (`list`)\cr see slot definition.
#' @param stop_report see [`Simulations`]
#' @param additional_stats (`list`)\cr see slot definition.
#' @param \dots additional parameters from [`GeneralSimulations`]
#'
#' @example examples/Simulations-class-Simulations.R
#' @export
Simulations <- function(fit, stop_reasons, stop_report, additional_stats, ...) {
  start <- GeneralSimulations(...)
  .Simulations(
    start,
    fit = fit,
    stop_report = stop_report,
    stop_reasons = stop_reasons,
    additional_stats = additional_stats
  )
}

## default constructor ----

#' @rdname Simulations-class
#' @note Typically, end users will not use the `.DefaultSimulations()` function.
#' @export
.DefaultSimulations <- function() {
  design <- .DefaultDesign()
  myTruth <- probFunction(design@model, alpha0 = 7, alpha1 = 8)

  simulate(
    design,
    args = NULL,
    truth = myTruth,
    nsim = 1,
    seed = 819,
    mcmcOptions = .DefaultMcmcOptions(),
    parallel = FALSE
  )
}

# HierarchicalSimulations ----

## class ----

#' `HierarchicalSimulations`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' This class captures trial simulations from hierarchical designs
#' ([`HierarchicalDesign`]).
#'
#' @slot data (`list`)\cr produced [`HierarchicalData`] objects.
#' @slot doses (`list`)\cr final dose recommendations, one named list per
#'   simulation run with entries for each hierarchical arm.
#' @slot samples (`list`)\cr final [`HierarchicalSamples`] objects.
#' @slot fit (`list`)\cr final arm-specific fits for each simulation run.
#' @slot stop_reasons (`list`)\cr arm-specific stopping reasons for each
#'   simulation run.
#' @slot stop_report (`list`)\cr arm-specific stopping rule outcomes for each
#'   simulation run.
#' @slot additional_stats (`list`)\cr additional arm-specific statistics.
#' @slot seed (`integer`)\cr random generator state before starting the
#'   simulation.
#'
#' @aliases HierarchicalSimulations
#' @export
.HierarchicalSimulations <- setClass(
  Class = "HierarchicalSimulations",
  slots = c(
    data = "list",
    doses = "list",
    samples = "list",
    fit = "list",
    stop_report = "list",
    stop_reasons = "list",
    additional_stats = "list",
    seed = "integer"
  ),
  prototype = prototype(
    data = list(.DefaultHierarchicalData()),
    doses = list(list(arm1 = 1, combo = c(drug1 = 10, drug2 = 20))),
    samples = list(.HierarchicalSamples()),
    fit = list(list()),
    stop_report = list(list()),
    stop_reasons = list(list()),
    additional_stats = list(list()),
    seed = 1L
  ),
  contains = "CrmPackClass",
  validity = v_hierarchical_simulations
)

## constructor ----

#' @rdname HierarchicalSimulations-class
#'
#' @param data (`list`)\cr see slot definition.
#' @param doses (`list`)\cr see slot definition.
#' @param samples (`list`)\cr see slot definition.
#' @param fit (`list`)\cr see slot definition.
#' @param stop_reasons (`list`)\cr see slot definition.
#' @param stop_report (`list`)\cr see slot definition.
#' @param additional_stats (`list`)\cr see slot definition.
#' @param seed (`integer`)\cr see slot definition.
#'
#' @export
HierarchicalSimulations <- function(
  data,
  doses,
  samples,
  fit,
  stop_reasons,
  stop_report,
  additional_stats,
  seed
) {
  assert_integerish(seed)
  .HierarchicalSimulations(
    data = data,
    doses = doses,
    samples = samples,
    fit = fit,
    stop_report = stop_report,
    stop_reasons = stop_reasons,
    additional_stats = additional_stats,
    seed = as.integer(seed)
  )
}

## default constructor ----

#' @rdname HierarchicalSimulations-class
#' @note Typically, end users will not use the
#'   `.DefaultHierarchicalSimulations()` function directly.
#' @export
.DefaultHierarchicalSimulations <- function() {
  design <- .DefaultHierarchicalDesign()
  truth <- lapply(
    design@arms,
    function(arm) {
      probFunction(arm@design@model, alpha0 = 7, alpha1 = 8)
    }
  )

  simulate(
    object = design,
    args = NULL,
    truth = truth,
    nsim = 1,
    seed = 819,
    mcmcOptions = .DefaultMcmcOptions(),
    parallel = FALSE
  )
}

# ComboSimulations ----

## class ----

#' `ComboSimulations`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' This class captures trial simulations from two-drug combination designs
#' ([`DesignCombo`]).
#'
#' It is intentionally separate from [`Simulations`] because
#' [`DataCombo`] is not a subclass of [`Data`], and final recommendations are
#' dose combinations.
#'
#' @slot data (`list`)
#'   produced [`DataCombo`] objects.
#' @slot doses (`matrix`)
#'   final recommended dose combinations; one row per simulation run and
#'   two columns (one per drug).
#' @slot fit (`list`)
#'   final fitted toxicity surfaces for each simulation run.
#' @slot stop_reasons (`list`)
#'   stopping reasons for each simulation run.
#' @slot stop_report (`matrix`)
#'   matrix of stopping rule outcomes.
#' @slot additional_stats (`list`)
#'   list of additional statistical summary values.
#' @slot seed (`integer`)
#'   random generator state before starting the simulation.
#' @aliases ComboSimulations
#' @export
.ComboSimulations <-
  setClass(
    Class = "ComboSimulations",
    slots = c(
      data = "list",
      doses = "matrix",
      fit = "list",
      stop_report = "matrix",
      stop_reasons = "list",
      additional_stats = "list",
      seed = "integer"
    ),
    prototype = prototype(
      data = list(
        DataCombo(
          x = cbind(drug1 = c(10, 10, 20), drug2 = c(20, 20, 20)),
          y = c(0L, 0L, 1L),
          ID = 1L:3L,
          cohort = c(1L, 1L, 2L),
          doseGrid = list(drug1 = c(10, 20, 30), drug2 = c(20, 40, 60))
        ),
        DataCombo(
          x = cbind(drug1 = c(10, 20, 20), drug2 = c(20, 20, 40)),
          y = c(0L, 0L, 0L),
          ID = 1L:3L,
          cohort = c(1L, 2L, 3L),
          doseGrid = list(drug1 = c(10, 20, 30), drug2 = c(20, 40, 60))
        )
      ),
      doses = cbind(drug1 = c(20, 20), drug2 = c(20, 40)),
      fit = list(
        data.frame(
          drug1 = 10,
          drug2 = 20,
          middle = 0.1,
          lower = 0.05,
          upper = 0.2
        ),
        data.frame(
          drug1 = 20,
          drug2 = 40,
          middle = 0.2,
          lower = 0.1,
          upper = 0.3
        )
      ),
      stop_report = matrix(TRUE, nrow = 2),
      stop_reasons = list("A", "B"),
      additional_stats = list(list(), list()),
      seed = 1L
    ),
    contains = "CrmPackClass",
    validity = v_combo_simulations
  )

## constructor ----

#' @rdname ComboSimulations-class
#'
#' @param data (`list`)
#'   see slot definition.
#' @param doses (`matrix`)
#'   see slot definition.
#' @param fit (`list`)
#'   see slot definition.
#' @param stop_reasons (`list`)
#'   see slot definition.
#' @param stop_report (`matrix`)
#'   see slot definition.
#' @param additional_stats (`list`)
#'   see slot definition.
#' @param seed (`integer`)
#'   see slot definition.
#' @export
ComboSimulations <- function(
  data,
  doses,
  fit,
  stop_reasons,
  stop_report,
  additional_stats,
  seed
) {
  if (!is.matrix(doses)) {
    doses <- as.matrix(doses)
  }

  assert_integerish(seed)

  .ComboSimulations(
    data = data,
    doses = doses,
    fit = fit,
    stop_report = stop_report,
    stop_reasons = stop_reasons,
    additional_stats = additional_stats,
    seed = as.integer(seed)
  )
}

## default constructor ----

#' @rdname ComboSimulations-class
#' @note Typically, end users will not use the `.DefaultComboSimulations()` function.
#' @export
.DefaultComboSimulations <- function() {
  ComboSimulations(
    data = list(
      DataCombo(
        x = cbind(drug1 = c(10, 10, 20), drug2 = c(20, 20, 20)),
        y = c(0L, 0L, 1L),
        ID = 1L:3L,
        cohort = c(1L, 1L, 2L),
        doseGrid = list(drug1 = c(10, 20, 30), drug2 = c(20, 40, 60))
      ),
      DataCombo(
        x = cbind(drug1 = c(10, 20, 20), drug2 = c(20, 20, 40)),
        y = c(0L, 0L, 0L),
        ID = 1L:3L,
        cohort = c(1L, 2L, 3L),
        doseGrid = list(drug1 = c(10, 20, 30), drug2 = c(20, 40, 60))
      )
    ),
    doses = cbind(drug1 = c(20, 20), drug2 = c(20, 40)),
    fit = list(
      data.frame(
        drug1 = 10,
        drug2 = 20,
        middle = 0.1,
        lower = 0.05,
        upper = 0.2
      ),
      data.frame(drug1 = 20, drug2 = 40, middle = 0.2, lower = 0.1, upper = 0.3)
    ),
    stop_report = matrix(c(TRUE, FALSE), nrow = 2),
    stop_reasons = list("A", "B"),
    additional_stats = list(list(), list()),
    seed = 123L
  )
}

# DualSimulations ----

## class ----

#' `DualSimulations`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This class captures the trial simulations from dual-endpoint model based
#' designs. In comparison to the parent class [`Simulations`],
#' it contains additional slots to capture the dose-biomarker `fits`, and the
#' `sigma2W` and `rho` estimates.
#'
#' @slot rho_est (`numeric`)\cr vector of final posterior median rho estimates
#' @slot sigma2w_est (`numeric`)\cr vector of final posterior median sigma2W estimates
#' @slot fit_biomarker (`list`)\cr with the final dose-biomarker curve fits
#' @aliases DualSimulations
#' @export
.DualSimulations <-
  setClass(
    Class = "DualSimulations",
    slots = c(
      rho_est = "numeric",
      sigma2w_est = "numeric",
      fit_biomarker = "list"
    ),
    prototype = prototype(
      rho_est = c(0.2, 0.3),
      sigma2w_est = c(0.2, 0.3),
      fit_biomarker = list(
        c(0.1, 0.2),
        c(0.1, 0.2)
      )
    ),
    contains = "Simulations",
    validity = v_dual_simulations
  )


## constructor ----

#' @rdname DualSimulations-class
#'
#' @param rho_est (`numeric`)\cr see [`DualSimulations`]
#' @param sigma2w_est (`numeric`)\cr [`DualSimulations`]
#' @param fit_biomarker (`list`)\cr see [`DualSimulations`]
#' @param \dots additional parameters from [`Simulations`]
#'
#' @example examples/Simulations-class-DualSimulations.R
#' @export
DualSimulations <- function(rho_est, sigma2w_est, fit_biomarker, ...) {
  start <- Simulations(...)
  .DualSimulations(
    start,
    rho_est = rho_est,
    sigma2w_est = sigma2w_est,
    fit_biomarker = fit_biomarker
  )
}

## default constructor ----

#' @rdname DualSimulations-class
#' @note Typically, end users will not use the `.DefaultDualSimulations()` function.
#' @export
.DefaultDualSimulations <- function() {
  DualSimulations(
    rho_est = c(0.25, 0.35),
    sigma2w_est = c(0.15, 0.25),
    fit_biomarker = list(c(0.3, 0.4), c(0.4, 0.5)),
    fit = list(
      c(0.1, 0.2),
      c(0.3, 0.4)
    ),
    stop_report = matrix(c(TRUE, FALSE), nrow = 2),
    stop_reasons = list("A", "B"),
    additional_stats = list(a = 1, b = 1),
    data = list(
      Data(
        x = 1:2,
        y = 0:1,
        doseGrid = 1:2,
        ID = 1L:2L,
        cohort = 1L:2L
      ),
      Data(
        x = 3:4,
        y = 0:1,
        doseGrid = 3:4,
        ID = 1L:2L,
        cohort = 1L:2L
      )
    ),
    doses = c(1, 2),
    seed = 123L
  )
}

#' `GeneralSimulationsSummary`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This class captures the summary of general simulations output. Note that objects
#' should not be created by users, therefore no initialization
#' function is provided for this class.
#'
#' @slot target (`numeric`)\cr target toxicity interval
#' @slot target_dose_interval (`numeric`)\cr corresponding target dose interval
#' @slot nsim (`integer`)\cr number of simulations
#' @slot prop_dlts (`ANY`)\cr A numeric array (multi-dimensional) or list representing proportions of DLTs in the trials
#' @slot mean_tox_risk (`numeric`)\cr mean toxicity risks for the patients
#' @slot dose_selected (`numeric`)\cr doses selected as MTD
#' @slot tox_at_doses_selected (`numeric`)\cr true toxicity at doses selected
#' @slot prop_at_target (`numeric`)\cr Proportion of trials selecting target MTD
#' @slot dose_most_selected (`numeric`)\cr dose most often selected as MTD
#' @slot obs_tox_rate_at_dose_most_selected (`numeric`)\cr observed toxicity rate at dose most often selected
#' @slot n_obs (`ANY`)\cr A numeric array (multi-dimensional) or list representing number of patients overall.
#' @slot n_above_target (`integer`)\cr number of patients treated above target tox interval
#' @slot dose_grid (`numeric`)\cr the dose grid that has been used
#' @slot placebo (`logical`)\cr set to TRUE (default is FALSE) for a design with placebo
#' @slot any_backfilled (`flag`)\cr indicates if any backfill cohorts were used
#' @slot n_backfill (`ANY`)\cr number of patients in backfill cohorts (only if `any_backfilled=TRUE`)
#' @slot backfill_doses (`ANY`)\cr list with doses used in backfill cohorts (only if `any_backfilled=TRUE`)
#' @aliases GeneralSimulationsSummary
#' @export
.GeneralSimulationsSummary <-
  setClass(
    Class = "GeneralSimulationsSummary",
    slots = c(
      target = "numeric",
      target_dose_interval = "numeric",
      nsim = "integer",
      prop_dlts = "ANY",
      mean_tox_risk = "numeric",
      dose_selected = "numeric",
      tox_at_doses_selected = "numeric",
      prop_at_target = "numeric",
      dose_most_selected = "numeric",
      obs_tox_rate_at_dose_most_selected = "numeric",
      n_obs = "ANY",
      n_above_target = "integer",
      dose_grid = "numeric",
      placebo = "logical",
      any_backfilled = "logical",
      n_backfill = "ANY",
      backfill_doses = "ANY"
    )
  )

## default constructor ----

#' @rdname GeneralSimulationsSummary-class
#' @note Typically, end users will not use the `.DefaultGeneralSimulationsSummary()` function.
#' @export
.DefaultGeneralSimulationsSummary <- function() {
  stop(
    paste(
      "Class GeneralSimulationsSummary cannot be instantiated directly.",
      "Please use one of its subclasses instead."
    )
  )
}

## SimulationsSummary ----

## class ----

#' `SimulationsSummary`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' In addition to the slots in the parent class [`GeneralSimulationsSummary`],
#' it contains two slots with model fit information.
#'
#' @slot stop_report (`matrix`)\cr matrix of stopping rule outcomes
#' @slot fit_at_dose_most_selected (`numeric`)\cr fitted toxicity rate at dose most often selected
#' @slot additional_stats (`list`)\cr list of additional statistical summary
#' @slot mean_fit (`list`)\cr list with the average, lower (2.5%) and upper (97.5%)
#' quantiles of the mean fitted toxicity at each dose level
#'
#' @aliases SimulationsSummary
#' @export
.SimulationsSummary <-
  setClass(
    Class = "SimulationsSummary",
    slots = c(
      stop_report = "matrix",
      fit_at_dose_most_selected = "numeric",
      additional_stats = "list",
      mean_fit = "list"
    ),
    contains = "GeneralSimulationsSummary"
  )

## default constructor ----

#' @rdname SimulationsSummary-class
#' @note Typically, end users will not use the `.DefaultSimulationsSummary()` function.
#' @export
.DefaultSimulationsSummary <- function() {
  stop(paste(
    "Class SimulationsSummary cannot be instantiated directly.",
    "Please use one of its subclasses instead."
  ))
}

# HierarchicalSimulationsSummary ----

## class ----

#' `HierarchicalSimulationsSummary`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' This class captures arm-specific summaries from [`HierarchicalSimulations`]
#' objects.
#'
#' @slot arms (`list`)\cr named list of arm-level simulation summary objects.
#' @slot nsim (`integer`)\cr number of simulations.
#' @aliases HierarchicalSimulationsSummary
#' @export
.HierarchicalSimulationsSummary <-
  setClass(
    Class = "HierarchicalSimulationsSummary",
    slots = c(
      arms = "list",
      nsim = "integer"
    ),
    contains = "CrmPackClass"
  )

## default constructor ----

#' @rdname HierarchicalSimulationsSummary-class
#' @note Typically, end users will not use the `.DefaultHierarchicalSimulationsSummary()` function.
#' @export
.DefaultHierarchicalSimulationsSummary <- function() {
  stop(paste(
    "Class HierarchicalSimulationsSummary cannot be instantiated directly.",
    "Please use summary(HierarchicalSimulations, ...) to create objects of this class."
  ))
}

# ComboSimulationsSummary ----

## class ----

#' `ComboSimulationsSummary`
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' This class captures summary output from [`ComboSimulations`] objects.
#'
#' @slot target (`numeric`)
#'   target toxicity interval. Empty if no truth function was supplied.
#' @slot nsim (`integer`)
#'   number of simulations.
#' @slot n_obs (`integer`)
#'   number of patients in each simulation.
#' @slot prop_dlts (`numeric`)
#'   observed proportion of DLTs in each simulation.
#' @slot mean_tox_risk (`numeric`)
#'   average fitted toxicity risk in each simulation.
#' @slot dose_selected (`matrix`)
#'   selected dose combinations; one row per simulation.
#' @slot tox_at_doses_selected (`numeric`)
#'   true toxicity at selected dose combinations (if truth supplied).
#' @slot prop_at_target (`numeric`)
#'   proportion of selected combinations within target interval (if truth supplied).
#' @slot dose_most_selected (`numeric`)
#'   most frequently selected dose combination.
#' @slot obs_tox_rate_at_dose_most_selected (`numeric`)
#'   observed toxicity rate at the most frequently selected combination.
#' @slot dose_grid (`list`)
#'   dose grid for each drug.
#' @slot stop_report (`matrix`)
#'   matrix of stopping rule outcomes.
#' @slot stop_reasons (`list`)
#'   stopping reasons by simulation.
#' @slot additional_stats (`list`)
#'   additional statistics.
#' @aliases ComboSimulationsSummary
#' @export
.ComboSimulationsSummary <-
  setClass(
    Class = "ComboSimulationsSummary",
    slots = c(
      target = "numeric",
      nsim = "integer",
      n_obs = "integer",
      prop_dlts = "numeric",
      mean_tox_risk = "numeric",
      dose_selected = "matrix",
      tox_at_doses_selected = "numeric",
      prop_at_target = "numeric",
      dose_most_selected = "numeric",
      obs_tox_rate_at_dose_most_selected = "numeric",
      dose_grid = "list",
      stop_report = "matrix",
      stop_reasons = "list",
      additional_stats = "list"
    )
  )

## default constructor ----

#' @rdname ComboSimulationsSummary-class
#' @note Typically, end users will not use the `.DefaultComboSimulationsSummary()` function.
#' @export
.DefaultComboSimulationsSummary <- function() {
  stop(
    paste(
      "Class ComboSimulationsSummary cannot be instantiated directly.",
      "Please use summary(ComboSimulations, ...) to create objects of this class."
    )
  )
}

# DualSimulationsSummary ----

# class ----

#' `DualSimulationsSummary`
#'
#' @description `r lifecycle::badge("stable")`
#' This class captures the summary of dual-endpoint simulations output.
#' In comparison to its parent class [`SimulationsSummary`], it has additional slots.
#'
#' @slot biomarker_fit_at_dose_most_selected (`numeric`)\cr fitted biomarker level at most often selected dose.
#' @slot mean_biomarker_fit (`list`)\cr list with average, lower (2.5%) and upper (97.5%) quantiles of
#' mean fitted biomarker level at each dose
#' @aliases DualSimulationsSummary
#' @export
.DualSimulationsSummary <-
  setClass(
    Class = "DualSimulationsSummary",
    slots = c(
      biomarker_fit_at_dose_most_selected = "numeric",
      mean_biomarker_fit = "list"
    ),
    contains = "SimulationsSummary"
  )

# default constructor

#' @rdname DualSimulationsSummary-class
#' @note Typically, end users will not use the `.DefaultDualSimulationsSummary()` function.
#' @export
.DefaultDualSimulationsSummary <- function() {
  empty_data <- DataDual(doseGrid = c(1, 3, 5, 10, 15, 20, 25, 30))

  my_model <- DualEndpointRW(
    mean = c(0, 1),
    cov = matrix(c(1, 0, 0, 1), nrow = 2),
    sigma2betaW = 0.01,
    sigma2W = c(a = 0.1, b = 0.1),
    rho = c(a = 1, b = 1),
    rw1 = TRUE
  )

  my_next_best <- NextBestDualEndpoint(
    target = c(0.9, 1),
    overdose = c(0.35, 1),
    max_overdose_prob = 0.25
  )

  my_size1 <- CohortSizeRange(
    intervals = c(0, 30),
    cohort_size = c(1, 3)
  )
  my_size2 <- CohortSizeDLT(
    intervals = c(0, 1),
    cohort_size = c(1, 3)
  )
  my_size <- maxSize(my_size1, my_size2)

  my_stopping1 <- StoppingTargetBiomarker(
    target = c(0.9, 1),
    prob = 0.5
  )

  my_stopping <- my_stopping1 | StoppingMinPatients(10) | StoppingMissingDose()

  my_increments <- IncrementsRelative(
    intervals = c(0, 20),
    increments = c(1, 0.33)
  )

  my_design <- DualDesign(
    model = my_model,
    data = empty_data,
    nextBest = my_next_best,
    stopping = my_stopping,
    increments = my_increments,
    cohort_size = CohortSizeConst(3),
    startingDose = 3
  )

  beta_mod <- function(dose, e0, eMax, delta1, delta2, scal) {
    maxDens <- (delta1^delta1) *
      (delta2^delta2) /
      ((delta1 + delta2)^(delta1 + delta2))
    dose <- dose / scal
    e0 + eMax / maxDens * (dose^delta1) * (1 - dose)^delta2
  }

  true_biomarker <- function(dose) {
    beta_mod(
      dose,
      e0 = 0.2,
      eMax = 0.6,
      delta1 = 5,
      delta2 = 5 * 0.5 / 0.5,
      scal = 100
    )
  }

  true_tox <- function(dose) {
    pnorm((dose - 60) / 10)
  }

  x <- simulate(
    object = my_design,
    trueTox = true_tox,
    trueBiomarker = true_biomarker,
    sigma2W = 0.01,
    rho = 0,
    nsim = 1,
    parallel = FALSE,
    seed = 3,
    startingDose = 6,
    mcmcOptions = .DefaultMcmcOptions()
  )
}

# PseudoSimulations ----

## class ----

#' `PseudoSimulations`
#'
#' @description `r lifecycle::badge("stable")`
#' This class captures trial simulations from designs using pseudo model.
#' It has additional slots `fit` and `stop_reasons` compared to the
#' general class [`GeneralSimulations`].
#'
#' @slot fit (`list`)\cr final fit values.
#' @slot final_td_target_during_trial_estimates (`numeric`)\cr final estimates of the `td_target_during_trial`.
#' @slot final_td_target_end_of_trial_estimates (`numeric`)\cr final estimates of the `td_target_end_of_trial`.
#' @slot final_td_target_during_trial_at_dose_grid (`numeric`)
#'        \cr dose levels at dose grid closest below the final `td_target_during_trial` estimates.
#' @slot final_td_target_end_of_trial_at_dose_grid (`numeric`)
#'        \cr dose levels at dose grid closest below the final `td_target_end_of_trial` estimates.
#' @slot final_tdeot_cis (`list`)\cr 95% credibility intervals of the final estimates for `td_target_end_of_trial`.
#' @slot final_tdeot_ratios (`numeric`)\cr ratio of the upper to the lower 95%
#'        credibility intervals for `td_target_end_of_trial`.
#' @slot final_cis (`list`)\cr final 95% credibility intervals for `td_target_end_of_trial` estimates.
#' @slot final_ratios (`numeric`)\cr final ratios of the upper to the lower 95%
#'        credibility interval for `td_target_end_of_trial`.
#' @slot stop_report (`matrix`)\cr outcomes of stopping rules.
#' @slot stop_reasons (`list`)\cr reasons for stopping each simulation run.
#'
#' @aliases PseudoSimulations
#' @export
.PseudoSimulations <-
  setClass(
    Class = "PseudoSimulations",
    slots = c(
      fit = "list",
      final_td_target_during_trial_estimates = "numeric",
      final_td_target_end_of_trial_estimates = "numeric",
      final_td_target_during_trial_at_dose_grid = "numeric",
      final_td_target_end_of_trial_at_dose_grid = "numeric",
      final_tdeot_cis = "list",
      final_tdeot_ratios = "numeric",
      final_cis = "list",
      final_ratios = "numeric",
      stop_report = "matrix",
      stop_reasons = "list"
    ),
    prototype = prototype(
      final_td_target_during_trial_estimates = c(0.1, 0.1),
      final_td_target_end_of_trial_estimates = c(0.1, 0.1),
      final_td_target_during_trial_at_dose_grid = c(0.1, 0.1),
      final_td_target_end_of_trial_at_dose_grid = c(0.1, 0.1),
      final_tdeot_cis = list(c(0.1, 0.2), c(0.1, 0.2)),
      final_tdeot_ratios = c(0.1, 0.1),
      final_cis = list(c(0.1, 0.2), c(0.1, 0.2)),
      final_ratios = c(0.1, 0.1),
      stop_report = matrix(TRUE, nrow = 2),
      stop_reasons = list("A", "A")
    ),
    contains = "GeneralSimulations",
    validity = v_pseudo_simulations
  )

## constructor ----

#' @rdname PseudoSimulations-class
#'
#' @param fit (`list`)\cr see slot definition.
#' @param final_td_target_during_trial_estimates (`numeric`)\cr see slot definition.
#' @param final_td_target_end_of_trial_estimates (`numeric`)\cr see slot definition.
#' @param final_td_target_during_trial_at_dose_grid (`numeric`)\cr see slot definition.
#' @param final_td_target_end_of_trial_at_dose_grid (`numeric`)\cr see slot definition.
#' @param final_tdeot_cis (`list`)\cr see slot definition.
#' @param final_tdeot_ratios (`numeric`)\cr see slot definition.
#' @param final_cis (`list`)\cr see slot definition.
#' @param final_ratios (`numeric`)\cr see slot definition.
#' @param stop_report see [`PseudoSimulations`]
#' @param stop_reasons (`list`)\cr see slot definition.
#' @param \dots additional parameters from [`GeneralSimulations`]
#'
#' @export
PseudoSimulations <- function(
  fit,
  final_td_target_during_trial_estimates,
  final_td_target_end_of_trial_estimates,
  final_td_target_during_trial_at_dose_grid,
  final_td_target_end_of_trial_at_dose_grid,
  final_tdeot_cis,
  final_tdeot_ratios,
  final_cis,
  final_ratios,
  stop_report,
  stop_reasons,
  ...
) {
  start <- GeneralSimulations(...)
  .PseudoSimulations(
    start,
    fit = fit,
    final_td_target_during_trial_estimates = final_td_target_during_trial_estimates,
    final_td_target_end_of_trial_estimates = final_td_target_end_of_trial_estimates,
    final_td_target_during_trial_at_dose_grid = final_td_target_during_trial_at_dose_grid,
    final_td_target_end_of_trial_at_dose_grid = final_td_target_end_of_trial_at_dose_grid,
    final_tdeot_cis = final_tdeot_cis,
    final_tdeot_ratios = final_tdeot_ratios,
    final_cis = final_cis,
    final_ratios = final_ratios,
    stop_report = stop_report,
    stop_reasons = stop_reasons
  )
}

## default constructor ----

#' @rdname PseudoSimulations-class
#' @note Typically, end users will not use the `.DefaultPseudoSimulations()` function.
#' @export
.DefaultPseudoSimulations <- function() {
  stop(
    "Class PseudoSimulations cannot be instantiated directly. Please use one of its subclasses instead."
  )
}

# PseudoDualSimulations ----

## class ----

#' `PseudoDualSimulations`
#'
#' @description `r lifecycle::badge("stable")`
#' This class conducts trial simulations for designs using both the
#' DLE and efficacy responses. It defines final values for
#' efficacy fit and DLE, estimates of Gstar, optimal dose and sigma2.
#'
#' @slot fit_eff (`list`)\cr final values of efficacy fit.
#' @slot final_gstar_estimates (`numeric`)\cr final Gstar estimates.
#' @slot final_gstar_at_dose_grid (`numeric`)\cr final Gstar estimates at dose grid.
#' @slot final_gstar_cis (`list`)\cr list of 95% confidence interval for Gstar estimates.
#' @slot final_gstar_ratios (`numeric`)\cr ratios of confidence intervals for Gstar estimates.
#' @slot final_optimal_dose (`numeric`)\cr final optimal dose.
#' @slot final_optimal_dose_at_dose_grid (`numeric`)\cr final optimal dose at dose grid.
#' @slot sigma2_est (`numeric`)\cr final sigma2 estimates.
#'
#' @aliases PseudoDualSimulations
#' @export
.PseudoDualSimulations <-
  setClass(
    Class = "PseudoDualSimulations",
    slots = c(
      fit_eff = "list",
      final_gstar_estimates = "numeric",
      final_gstar_at_dose_grid = "numeric",
      final_gstar_cis = "list",
      final_gstar_ratios = "numeric",
      final_optimal_dose = "numeric",
      final_optimal_dose_at_dose_grid = "numeric",
      sigma2_est = "numeric"
    ),
    prototype = prototype(
      final_gstar_estimates = c(0.1, 0.1),
      final_gstar_at_dose_grid = c(0.1, 0.1),
      final_gstar_cis = list(c(0.1, 0.2), c(0.1, 0.2)),
      final_gstar_ratios = c(0.01, 0.01),
      final_optimal_dose = c(0.01, 0.01),
      final_optimal_dose_at_dose_grid = c(0.01, 0.01),
      sigma2_est = c(0.001, 0.002)
    ),
    contains = "PseudoSimulations",
    validity = v_pseudo_dual_simulations
  )

## constructor ----

#' @rdname PseudoDualSimulations-class
#'
#' @param fit_eff (`list`)\cr see slot definition.
#' @param final_gstar_estimates (`numeric`)\cr see slot definition.
#' @param final_gstar_at_dose_grid (`numeric`)\cr see slot definition.
#' @param final_gstar_cis (`list`)\cr see slot definition.
#' @param final_gstar_ratios (`numeric`)\cr see slot definition.
#' @param final_optimal_dose (`numeric`)\cr see slot definition.
#' @param final_optimal_dose_at_dose_grid (`numeric`)\cr see slot definition.
#' @param sigma2_est (`numeric`)\cr see slot definition.
#' @param \dots additional parameters from [`PseudoSimulations`]
#' @export
PseudoDualSimulations <- function(
  fit_eff,
  final_gstar_estimates,
  final_gstar_at_dose_grid,
  final_gstar_cis,
  final_gstar_ratios,
  final_optimal_dose,
  final_optimal_dose_at_dose_grid,
  sigma2_est,
  ...
) {
  start <- PseudoSimulations(...)
  .PseudoDualSimulations(
    start,
    fit_eff = fit_eff,
    final_gstar_estimates = final_gstar_estimates,
    final_gstar_at_dose_grid = final_gstar_at_dose_grid,
    final_gstar_cis = final_gstar_cis,
    final_gstar_ratios = final_gstar_ratios,
    final_optimal_dose = final_optimal_dose,
    final_optimal_dose_at_dose_grid = final_optimal_dose_at_dose_grid,
    sigma2_est = sigma2_est
  )
}

## default constructor ----

#' @rdname PseudoDualSimulations-class
#' @note Do not use the `.DefaultPseudoDualSimulations()` function.
#' @export
.DefaultPseudoDualSimulations <- function() {
  stop(
    "Class PseudoDualSimulations cannot be instantiated directly. Please use a subclass."
  )
}

# PseudoDualFlexiSimulations ----

## class ----

#' `PseudoDualFlexiSimulations`
#'
#' @description `r lifecycle::badge("stable")`
#' This class captures the trial simulations design using both the DLE and
#' efficacy responses using [`EffFlexi`] efficacy model.
#' It extends [`PseudoDualSimulations`] by adding the capability to capture the sigma2betaW estimates.
#'
#' @slot sigma2_beta_w_est (`numeric`)\cr the vector of the final posterior mean sigma2betaW estimates
#' @aliases PseudoDualFlexiSimulations
#' @export
.PseudoDualFlexiSimulations <-
  setClass(
    Class = "PseudoDualFlexiSimulations",
    slots = c(sigma2_beta_w_est = "numeric"),
    prototype = prototype(sigma2_beta_w_est = c(0.001, 0.002)),
    contains = "PseudoDualSimulations"
  )

## constructor ----

#' @rdname PseudoDualFlexiSimulations-class
#'
#' @param sigma2_beta_w_est (`numeric`)\cr the vector of the final posterior mean sigma2betaW estimates
#' @param \dots additional parameters from [`PseudoDualSimulations`]
#'
#' @export
PseudoDualFlexiSimulations <- function(sigma2_beta_w_est, ...) {
  start <- PseudoDualSimulations(...)
  .PseudoDualFlexiSimulations(start, sigma2_beta_w_est = sigma2_beta_w_est)
}

## default constructor ----

#' @rdname PseudoDualFlexiSimulations-class
#' @note Typically, end users will not use the `.DefaultPseudoFlexiSimulations()` function.
#' @export
.DefaultPseudoDualFlexiSimulations <- function() {
  stop(
    "Class PseudoFlexiSimulations cannot be instantiated directly. Please use one of its subclasses instead."
  )
}

# PseudoSimulationsSummary ----

## class ----

#' `PseudoSimulationsSummary`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This class captures the summary of pseudo-models simulations output.
#' Note that objects should not be created by users, therefore no
#' initialization function is provided for this class.
#'
#' @slot target_end_of_trial (`numeric`)\cr the target probability of DLE wanted at the end of a trial
#' @slot target_dose_end_of_trial (`numeric`)\cr the dose level corresponds to the target probability
#'   of DLE wanted at the end of a trial, TDEOT
#' @slot target_dose_end_of_trial_at_dose_grid (`numeric`)\cr the dose level at dose grid corresponds to the
#'   target probability of DLE wanted at the end of a trial
#' @slot target_during_trial (`numeric`)\cr the target probability of DLE wanted during a trial
#' @slot target_dose_during_trial (`numeric`)\cr the dose level corresponds to the target probability of DLE
#'   wanted during the trial. TDDT
#' @slot target_dose_during_trial_at_dose_grid (`numeric`)\cr the dose level at dose grid corresponds to the
#'   target probability of DLE wanted during a trial
#' @slot tdeot_summary (`table`)\cr the six-number table summary, include the lowest, the 25th percentile
#'   (lower quartile), the 50th percentile (median), the mean, the 75th percentile and the highest values of the
#'   final dose levels obtained corresponds to the target probability of DLE
#'   want at the end of a trial across all simulations
#' @slot tddt_summary (`table`)\cr the six-number table summary, include the lowest, the 25th percentile
#'   (lower quartile), the 50th percentile (median), the mean, the 75th percentile and the highest values of the
#'   final dose levels obtained corresponds to the target probability of DLE
#'   want during a trial across all simulations
#' @slot final_dose_rec_summary (`table`)\cr the six-number table summary, include the lowest, the 25th percentile
#'   (lower quartile), the 50th percentile (median), the mean, the 75th percentile and the highest values of the
#'   final optimal doses, which is either the TDEOT when only DLE response are incorporated into
#'   the escalation procedure or the minimum of the TDEOT and Gstar when DLE and efficacy responses are
#'   incorporated, across all simulations
#' @slot ratio_tdeot_summary (`table`)\cr the six-number summary table of the final ratios of the upper to the
#'   lower 95% credibility intervals of the final TDEOTs across all simulations
#' @slot final_ratio_summary (`table`)\cr the six-number summary table of the final ratios of the upper to the
#'   lower 95% credibility intervals of the final optimal doses across all simulations
#' @slot nsim (`integer`)\cr number of simulations
#' @slot prop_dle (`numeric`)\cr proportions of DLE in the trials
#' @slot mean_tox_risk (`numeric`)\cr mean toxicity risks for the patients
#' @slot dose_selected (`numeric`)\cr doses selected as MTD (target_dose_end_of_trial)
#' @slot tox_at_doses_selected (`numeric`)\cr true toxicity at doses selected
#' @slot prop_at_target_end_of_trial (`numeric`)\cr Proportion of trials selecting at the dose_grid closest below the
#'   MTD, the target_dose_end_of_trial
#' @slot prop_at_target_during_trial (`numeric`)\cr Proportion of trials selecting at the dose_grid closest below
#'   the target_dose_during_trial
#' @slot dose_most_selected (`numeric`)\cr dose most often selected as MTD
#' @slot obs_tox_rate_at_dose_most_selected (`numeric`)\cr observed toxicity rate at dose most often
#'   selected
#' @slot n_obs (`integer`)\cr number of patients overall
#' @slot n_above_target_end_of_trial (`integer`)\cr number of patients treated above target_dose_end_of_trial
#' @slot n_above_target_during_trial (`integer`)\cr number of patients treated above target_dose_during_trial
#' @slot dose_grid (`numeric`)\cr the dose grid that has been used
#' @slot fit_at_dose_most_selected (`numeric`)\cr fitted toxicity rate at dose most often selected
#' @slot mean_fit (`list`)\cr list with the average, lower (2.5%) and upper (97.5%)
#'   quantiles of the mean fitted toxicity at each dose level
#' @slot stop_report (`matrix`)\cr matrix of stopping rule outcomes
#'
#' @aliases PseudoSimulationsSummary
#' @export
.PseudoSimulationsSummary <-
  setClass(
    Class = "PseudoSimulationsSummary",
    slots = c(
      target_end_of_trial = "numeric",
      target_dose_end_of_trial = "numeric",
      target_dose_end_of_trial_at_dose_grid = "numeric",
      target_during_trial = "numeric",
      target_dose_during_trial = "numeric",
      target_dose_during_trial_at_dose_grid = "numeric",
      tdeot_summary = "table",
      tddt_summary = "table",
      final_dose_rec_summary = "table",
      ratio_tdeot_summary = "table",
      final_ratio_summary = "table",
      nsim = "integer",
      prop_dle = "numeric",
      mean_tox_risk = "numeric",
      dose_selected = "numeric",
      tox_at_doses_selected = "numeric",
      prop_at_target_end_of_trial = "numeric",
      prop_at_target_during_trial = "numeric",
      dose_most_selected = "numeric",
      obs_tox_rate_at_dose_most_selected = "numeric",
      n_obs = "integer",
      n_above_target_end_of_trial = "integer",
      n_above_target_during_trial = "integer",
      dose_grid = "numeric",
      fit_at_dose_most_selected = "numeric",
      mean_fit = "list",
      stop_report = "matrix"
    )
  )

## default constructor ----

#' @rdname PseudoSimulationsSummary-class
#' @note Typically, end users will not use the `.DefaultPseudoSimulationsSummary()` function.
#' @export
.DefaultPseudoSimulationsSummary <- function() {
  stop(
    "Class PseudoSimulationsSummary cannot be instantiated directly. Please use one of its subclasses instead."
  )
}

# PseudoDualSimulationsSummary ----

## class ----

#' `PseudoDualSimulationsSummary`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This class captures the summary of the dual responses simulations using pseudo models.
#' It contains all slots from [`PseudoSimulationsSummary`] object. In addition to
#' the slots in the parent class [`PseudoSimulationsSummary`], it contains additional
#' slots for the efficacy model fit information.
#'
#' Note that objects should not be created by users, therefore no initialization function
#' is provided for this class.
#'
#' @slot target_gstar (`numeric`)\cr the target dose level such that its gain value is at maximum
#' @slot target_gstar_at_dose_grid (`numeric`)\cr the dose level at dose Grid closest and below Gstar
#' @slot gstar_summary (`table`)\cr the six-number table summary (lowest, 25th, 50th (median), 75th percentile, mean
#'   and highest value) of the final Gstar values obtained across all simulations
#' @slot ratio_gstar_summary (`table`)\cr the six-number summary table of the ratios of the upper to the lower 95%
#'   credibility intervals of the final Gstar across all simulations
#' @slot eff_fit_at_dose_most_selected (`numeric`)\cr fitted expected mean efficacy value at dose most often
#'   selected
#' @slot mean_eff_fit (`list`)\cr list with mean, lower (2.5%) and upper (97.5%) quantiles of the fitted expected
#'   efficacy value at each dose level.
#'
#' @aliases PseudoDualSimulationsSummary
#' @export
.PseudoDualSimulationsSummary <-
  setClass(
    Class = "PseudoDualSimulationsSummary",
    contains = "PseudoSimulationsSummary",
    slots = c(
      target_gstar = "numeric",
      target_gstar_at_dose_grid = "numeric",
      gstar_summary = "table",
      ratio_gstar_summary = "table",
      eff_fit_at_dose_most_selected = "numeric",
      mean_eff_fit = "list"
    )
  )

## default constructor ----

#' @rdname PseudoDualSimulationsSummary-class
#' @note Typically, end users will not use the `.DefaultPseudoDualSimulationsSummary()` function.
#' @export
.DefaultPseudoDualSimulationsSummary <- function() {
  stop(
    "Class PseudoDualSimulationsSummary cannot be instantiated directly. Please use one of its subclasses instead."
  )
}

# DASimulations ----

## class ----

#' `DASimulations`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This class captures the trial simulations from DA based designs.
#' In comparison to the parent class [`Simulations`],
#' it contains additional slots to capture the time to DLT fits, additional
#' parameters and the trial duration.
#'
#' @slot trial_duration (`numeric`)\cr the vector of trial duration values for all simulations.
#'
#' @aliases DASimulations
#' @export
.DASimulations <-
  setClass(
    Class = "DASimulations",
    slots = c(trial_duration = "numeric"),
    prototype = prototype(trial_duration = rep(0, 2)),
    contains = "Simulations",
    validity = v_da_simulations
  )

## constructor ----

#' @rdname DASimulations-class
#'
#' @param trial_duration (`numeric`)\cr see [`DASimulations`]
#' @param \dots additional parameters from [`Simulations`]
#'
#' @export
DASimulations <- function(trial_duration, ...) {
  start <- Simulations(...)
  .DASimulations(start, trial_duration = trial_duration)
}


## default constructor ----

#' @rdname DASimulations-class
#' @note Typically, end users will not use the `.DefaultDASimulations()` function.
#' @export
.DefaultDASimulations <- function() {
  design <- .DefaultDADesign()
  myTruth <- probFunction(design@model, alpha0 = 2, alpha1 = 3)
  exp_cond_cdf <- function(x, onset = 15) {
    a <- stats::pexp(28, 1 / onset, lower.tail = FALSE)
    1 - (stats::pexp(x, 1 / onset, lower.tail = FALSE) - a) / (1 - a)
  }

  simulate(
    design,
    args = NULL,
    truthTox = myTruth,
    truthSurv = exp_cond_cdf,
    trueTmax = 80,
    nsim = 2,
    seed = 819,
    mcmcOptions = .DefaultMcmcOptions(),
    firstSeparate = TRUE,
    deescalate = FALSE,
    parallel = FALSE
  )
}

# tidy

## tidy-Simulations ----

#' @rdname tidy
#' @aliases tidy-Simulations
#' @example examples/Simulations-method-tidy.R
#' @export
setMethod(
  f = "tidy",
  signature = signature(x = "Simulations"),
  definition = function(x, ...) {
    slot_names <- slotNames(x)
    rv <- list()
    for (nm in slot_names) {
      if (!is.function(slot(x, nm))) {
        if (nm %in% c("stop_reasons", "additional_stats")) {} else {
          rv[[nm]] <- h_tidy_slot(x, nm)
        }
      }
    }
    # Column bind of all list elements have the same number of rows
    if (length(rv) > 1 & length(unique(sapply(rv, nrow))) == 1) {
      rv <- rv %>% dplyr::bind_cols()
    }
    rv %>% h_tidy_class(x)
  }
)

## tidy-ComboSimulations ----

#' @rdname tidy
#' @aliases tidy-ComboSimulations
#' @export
setMethod(
  f = "tidy",
  signature = signature(x = "ComboSimulations"),
  definition = function(x, ...) {
    nsim <- length(x@data)
    simulations <- seq_len(nsim)

    data <- lapply(
      simulations,
      function(i) {
        tidy(x@data[[i]]) %>%
          tibble::add_column(Simulation = i, .before = 1L)
      }
    ) %>%
      dplyr::bind_rows()

    doses <- tibble::as_tibble(x@doses) %>%
      tibble::add_column(Simulation = simulations, .before = 1L)

    fit <- lapply(
      simulations,
      function(i) {
        tibble::as_tibble(x@fit[[i]]) %>%
          tibble::add_column(Simulation = i, .before = 1L)
      }
    ) %>%
      dplyr::bind_rows()

    stop_report <- tibble::as_tibble(x@stop_report, .name_repair = "unique")
    if (is.null(names(stop_report))) {
      names(stop_report) <- paste0("Rule", seq_along(stop_report))
    }
    stop_report <- stop_report %>%
      tibble::add_column(Simulation = simulations, .before = 1L)

    stop_reasons <- tibble::tibble(
      Simulation = simulations,
      StopReason = x@stop_reasons
    )

    additional_stats <- tibble::tibble(
      Simulation = simulations,
      AdditionalStats = x@additional_stats
    )

    list(
      data = data,
      doses = doses,
      fit = fit,
      stop_report = stop_report,
      stop_reasons = stop_reasons,
      additional_stats = additional_stats,
      seed = tibble::tibble(Seed = list(x@seed))
    ) %>%
      h_tidy_class(x)
  }
)

## tidy-HierarchicalSimulations ----

#' @rdname tidy
#' @aliases tidy-HierarchicalSimulations
#' @export
setMethod(
  f = "tidy",
  signature = signature(x = "HierarchicalSimulations"),
  definition = function(x, ...) {
    nsim <- length(x@data)
    simulations <- seq_len(nsim)

    data <- lapply(
      simulations,
      function(i) {
        tidy(x@data[[i]]) %>%
          tibble::add_column(Simulation = i, .before = 1L)
      }
    ) %>%
      dplyr::bind_rows()

    doses <- lapply(
      simulations,
      function(i) {
        dose_list <- x@doses[[i]]
        tibble::tibble(
          Simulation = i,
          Arm = names(dose_list),
          Dose = unname(dose_list)
        )
      }
    ) %>%
      dplyr::bind_rows()

    samples <- tibble::tibble(
      Simulation = simulations,
      Samples = lapply(x@samples, tidy)
    )

    fit <- tibble::tibble(
      Simulation = simulations,
      Fit = x@fit
    )

    stop_report <- tibble::tibble(
      Simulation = simulations,
      StopReport = x@stop_report
    )

    stop_reasons <- tibble::tibble(
      Simulation = simulations,
      StopReason = x@stop_reasons
    )

    additional_stats <- tibble::tibble(
      Simulation = simulations,
      AdditionalStats = x@additional_stats
    )

    list(
      data = data,
      doses = doses,
      samples = samples,
      fit = fit,
      stop_report = stop_report,
      stop_reasons = stop_reasons,
      additional_stats = additional_stats,
      seed = tibble::tibble(Seed = list(x@seed))
    ) %>%
      h_tidy_class(x)
  }
)