inst/doc/count.R

In SimTOST: Sample Size Estimation for Bio-Equivalence Trials Through Simulation

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

## ----poisson-literature-benchmark---------------------------------------------
zhu_benchmark <- simPower(
  n = 2705,
  distribution = "pois",
  rate_list = list(TEST = 1, REF = 1),
  list_comparator = list(TEST_vs_REF = c("TEST", "REF")),
  list_lequi.tol = list(TEST_vs_REF = 0.9),
  list_uequi.tol = list(TEST_vs_REF = 1 / 0.9),
  exposure = 0.7,
  dtype = "parallel",
  alpha = 0.025,
  nsim = 1000,
  seed = 2024
)
zhu_benchmark
data.frame(reference_power = 0.80012,
           simulated_power = zhu_benchmark$power,
           simulated_lower = zhu_benchmark$power_LCI,
           simulated_upper = zhu_benchmark$power_UCI)

## ----poisson-inputs-----------------------------------------------------------
count_corr <- matrix(c(1, 0.5, 0.5, 1), nrow = 2,
                     dimnames = list(c("y1", "y2"), c("y1", "y2")))
count_rates <- list(TEST = c(y1 = 0.21, y2 = 0.24),
                    REF = c(y1 = 0.20, y2 = 0.22))
count_comparators <- list(TEST_vs_REF = c("TEST", "REF"))
count_lower <- list(TEST_vs_REF = c(y1 = 0.80, y2 = 0.80))
count_upper <- list(TEST_vs_REF = c(y1 = 1.25, y2 = 1.25))
count_exposure <- 10

## ----poisson-diagnostics------------------------------------------------------
poisson_diagnostic <- sampleSize(
  distribution = "pois", rate_list = count_rates,
  list_comparator = count_comparators,
  list_lequi.tol = count_lower, list_uequi.tol = count_upper,
  exposure = count_exposure, cor_mat = count_corr,
  dtype = "parallel", nsim = 500, seed = 1234,
  keep_sim_data = TRUE
)
poisson_diagnostic
confint(poisson_diagnostic)

## ----poisson-rate-distribution, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
plot_distribution(poisson_diagnostic, estimand = "rate")

## ----poisson-correlation-distribution, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
plot_distribution(poisson_diagnostic, estimand = "correlation",
                  arms = c("TEST", "REF"))

## ----poisson-outcome-distribution, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
plot_distribution(poisson_diagnostic, estimand = "outcome", type = "histogram")

## ----poisson-estimand-distribution, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
plot_distribution(poisson_diagnostic, estimand = "RR")

## ----poisson-stability, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
plot_stability(poisson_diagnostic)
plot_mc_error(poisson_diagnostic)

## ----poisson-type1, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
type1_one <- type1Error(x = poisson_diagnostic, null = "both", joint = TRUE)
plot(type1_one)

## ----poisson-sample-size------------------------------------------------------
poisson_sample_size <- sampleSize(
  power = 0.80, distribution = "pois", rate_list = count_rates,
  list_comparator = count_comparators,
  list_lequi.tol = count_lower, list_uequi.tol = count_upper,
  exposure = count_exposure, cor_mat = count_corr, dtype = "parallel",
  lower = 100, upper = 1000, nsim = 1000, seed = 1234, ncores = 1,
  keep_sim_data = TRUE
)
summary(poisson_sample_size)
confint(poisson_sample_size)

## ----poisson-result, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
plot(poisson_sample_size)

## ----negative-binomial-sample-size--------------------------------------------
nb_dispersion <- 0.50
negative_binomial_sample_size <- update(poisson_sample_size,
                                        distribution = "nbinom",
                                        dispersion = nb_dispersion)

summary(negative_binomial_sample_size)
confint(negative_binomial_sample_size)

## ----nb-outcome-distribution, fig.width=8, fig.height=5, out.width="100%", fig.align="center"----
plot_distribution(negative_binomial_sample_size, estimand = "outcome", type = "histogram")

## ----crossover-negative-binomial-sample-size----------------------------------

crossover_sample_size <- update(negative_binomial_sample_size,
                                dtype = "2x2",
                                sigmaB = 0.30,
                                Eper = c(0, 0.10),
                                Eco = c(0, 0),
                                dropout = c(0.10, 0.10))

summary(crossover_sample_size)
confint(crossover_sample_size)