## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, comment = "#>",
                      warning = FALSE, message = FALSE)
library(SimTOST)

## ----rates--------------------------------------------------------------------
rate_list <- list(
  TEST = c(Exacerbations = 0.19,
           Hospitalizations = 0.13,
           RescueEvents = 0.08),
  EU_REF = c(Exacerbations = 0.20,
             Hospitalizations = 0.13,
             RescueEvents = 0.08),
  US_REF = c(Exacerbations = 0.21,
             Hospitalizations = 0.12,
             RescueEvents = 0.08)
)

comparisons <- list(
  EU_comparison = c("TEST", "EU_REF"),
  US_comparison = c("TEST", "US_REF")
)

lower_margin <- 0.80
upper_margin <- 1.25
exposure <- 5

## ----fixed-power--------------------------------------------------------------
run_fixed_power <- function(comparison_name, endpoint_name) {
  arms <- comparisons[[comparison_name]]
  result <- simPower(
    n = 100,
    distribution = "pois",
    rate_list = setNames(list(
      setNames(rate_list[[arms[1]]][[endpoint_name]], endpoint_name),
      setNames(rate_list[[arms[2]]][[endpoint_name]], endpoint_name)
    ), arms),
    list_comparator = setNames(list(arms), comparison_name),
    list_lequi.tol = setNames(list(lower_margin), comparison_name),
    list_uequi.tol = setNames(list(upper_margin), comparison_name),
    exposure = exposure,
    dtype = "parallel",
    nsim = 1000,
    seed = 1234
  )
  data.frame(
    comparison = comparison_name, endpoint = endpoint_name,
    power = result$power, power_LCI = result$power_LCI,
    power_UCI = result$power_UCI
  )
}

# Each object is one transparent comparison-endpoint calculation.
fixed_power_EU_Exacerbations <- run_fixed_power("EU_comparison", "Exacerbations")
fixed_power_EU_Hospitalizations <- run_fixed_power("EU_comparison", "Hospitalizations")
fixed_power_EU_RescueEvents <- run_fixed_power("EU_comparison", "RescueEvents")
fixed_power_US_Exacerbations <- run_fixed_power("US_comparison", "Exacerbations")
fixed_power_US_Hospitalizations <- run_fixed_power("US_comparison", "Hospitalizations")
fixed_power_US_RescueEvents <- run_fixed_power("US_comparison", "RescueEvents")

fixed_power <- rbind(
  fixed_power_EU_Exacerbations,
  fixed_power_EU_Hospitalizations,
  fixed_power_EU_RescueEvents,
  fixed_power_US_Exacerbations,
  fixed_power_US_Hospitalizations,
  fixed_power_US_RescueEvents
)

fixed_power

## ----sample-size--------------------------------------------------------------
run_sample_size <- function(comparison_name, endpoint_name) {
  arms <- comparisons[[comparison_name]]
  result <- sampleSize(
    power = 0.80,
    distribution = "pois",
    rate_list = setNames(list(
      setNames(rate_list[[arms[1]]][[endpoint_name]], endpoint_name),
      setNames(rate_list[[arms[2]]][[endpoint_name]], endpoint_name)
    ), arms),
    list_comparator = setNames(list(arms), comparison_name),
    list_lequi.tol = setNames(list(lower_margin), comparison_name),
    list_uequi.tol = setNames(list(upper_margin), comparison_name),
    exposure = exposure,
    dtype = "parallel",
    nsim = 1000,
    seed = 1234,
    lower = 10,
    upper = 2000
  )
  data.frame(
    comparison = comparison_name, endpoint = endpoint_name,
    n_per_arm = result$n_per_arm, n_total_for_pair = result$n_total,
    achieved_power = result$power
  )
}

# Again, keep the six calculations as named objects so each result can be
# inspected independently before combining them into one table.
sample_size_EU_Exacerbations <- run_sample_size("EU_comparison", "Exacerbations")
sample_size_EU_Hospitalizations <- run_sample_size("EU_comparison", "Hospitalizations")
sample_size_EU_RescueEvents <- run_sample_size("EU_comparison", "RescueEvents")
sample_size_US_Exacerbations <- run_sample_size("US_comparison", "Exacerbations")
sample_size_US_Hospitalizations <- run_sample_size("US_comparison", "Hospitalizations")
sample_size_US_RescueEvents <- run_sample_size("US_comparison", "RescueEvents")

sample_size_results <- rbind(
  sample_size_EU_Exacerbations,
  sample_size_EU_Hospitalizations,
  sample_size_EU_RescueEvents,
  sample_size_US_Exacerbations,
  sample_size_US_Hospitalizations,
  sample_size_US_RescueEvents
)

sample_size_results

## ----conservative-size--------------------------------------------------------
required_per_arm <- max(sample_size_results$n_per_arm)
required_total <- 3 * required_per_arm

c(required_per_arm = required_per_arm, required_total = required_total)

## ----joint-sample-size--------------------------------------------------------
endpoint_corr <- matrix(c(
  1.0, 0.40, 0.25,
  0.40, 1.0, 0.35,
  0.25, 0.35, 1.0
), nrow = 3, byrow = TRUE)

joint_result <- sampleSize(
  power = 0.80,
  distribution = "pois",
  rate_list = rate_list,
  list_comparator = comparisons,
  list_lequi.tol = list(
    EU_comparison = rep(lower_margin, 3),
    US_comparison = rep(lower_margin, 3)
  ),
  list_uequi.tol = list(
    EU_comparison = rep(upper_margin, 3),
    US_comparison = rep(upper_margin, 3)
  ),
  exposure = rep(exposure, 3),
  cor_mat = endpoint_corr,
  dtype = "parallel",
  nsim = 500,
  seed = 1234,
  lower = 10,
  upper = 3000,
  k = 3,
  adjust = "none"
)

joint_sample_size <- data.frame(
  n_per_arm = joint_result$n_per_arm,
  n_total = joint_result$n_total,
  achieved_power = joint_result$power,
  k = joint_result$k,
  adjustment = joint_result$adjust
)

joint_sample_size

## ----joint-advantage----------------------------------------------------------
joint_power_at_separate <- simPower(
  n = required_per_arm,
  distribution = "pois",
  rate_list = rate_list,
  list_comparator = comparisons,
  list_lequi.tol = list(
    EU_comparison = rep(lower_margin, 3),
    US_comparison = rep(lower_margin, 3)
  ),
  list_uequi.tol = list(
    EU_comparison = rep(upper_margin, 3),
    US_comparison = rep(upper_margin, 3)
  ),
  exposure = rep(exposure, 3),
  cor_mat = endpoint_corr,
  dtype = "parallel",
  nsim = 1000,
  seed = 1234,
  k = 3,
  adjust = "none"
)

joint_independent_result <- update(
  joint_result,
  cor_mat = diag(3),
  nsim = 500,
  seed = 1234
)

joint_comparison <- data.frame(
  approach = c(
    "Separate searches; joint power evaluated afterward",
    "Joint search; independent endpoints",
    "Joint search; correlated endpoints"
  ),
  n_per_arm = c(
    required_per_arm,
    joint_independent_result$n_per_arm,
    joint_result$n_per_arm
  ),
  joint_power = c(
    joint_power_at_separate$power,
    joint_independent_result$power,
    joint_result$power
  )
)

joint_comparison

## ----negative-binomial-sensitivity--------------------------------------------
nb_dispersion <- 0.10

joint_negative_binomial_result <- update(
  joint_result,
  distribution = "nbinom",
  dispersion = nb_dispersion,
  nsim = 500,
  seed = 1234
)

nb_sample_size <- data.frame(
  n_per_arm = joint_negative_binomial_result$n_per_arm,
  n_total = joint_negative_binomial_result$n_total,
  achieved_power = joint_negative_binomial_result$power,
  k = joint_negative_binomial_result$k,
  adjustment = joint_negative_binomial_result$adjust,
  dispersion = nb_dispersion
)

nb_sample_size