---
title: "Count-Outcome Equivalence: Parallel Design with 3 Arms and 3 Endpoints"
output: rmarkdown::html_vignette
bibliography: references.bib
link-citations: yes
vignette: >
  %\VignetteIndexEntry{Count-Outcome Equivalence: Parallel Design with 3 Arms and 3 Endpoints}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Overview

This example illustrates how to plan a parallel count-outcome equivalence
study with three treatment arms and three count endpoints. The arms are a test
product (`TEST`) and two reference products (`EU_REF` and `US_REF`). We use
three clinically interpretable, illustrative endpoint names: `Exacerbations`,
`Hospitalizations`, and `RescueEvents` (rescue-medication events). These labels
represent event counts collected over the same follow-up period. The general
distributional assumptions used by all
SimTOST outcomes are described in the companion vignette
[`methodological_assumptions.Rmd`](methodological_assumptions.Rmd).

The count module supports scalar or vector-valued endpoint inputs. A joint
simulation can require equivalence for `k` of the `m` endpoints and can apply
Bonferroni, Mielke's weak k-out-of-m (`adjust = "k"`),
Mielke's strong k-out-of-m (`adjust = "t"`), or Šidák adjustment
to the endpoint-wise one-sided significance level.

The joint count kernel can also simulate correlated endpoints through
the `cor_mat` argument. This dependence is generated with a Gaussian-copula
construction [@nelsen2006].

## Study assumptions

The assumed event rates are:

```{r 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
```

The rate-ratio equivalence interval is 0.80 to 1.25. We use a Poisson model
for the primary example and a one-sided significance level of 0.05 for each
TOST component.
The Poisson and negative-binomial planning framework follows published
count-outcome equivalence methods [@chang_sample_2017; @zhu_sample_2017].

## Power for a fixed sample size

Before estimating sample size, power can be examined for each endpoint and
comparison. Here we use 100 participants per treatment arm and a small number
of simulations to keep the vignette fast.

```{r 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 estimation for all comparisons and endpoints

We now estimate the smallest sample size per arm that reaches 80% power for
each individual comparison and endpoint.

```{r 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
```

For a simple conservative planning rule, select the maximum required sample
size per arm across all six comparison-endpoint combinations:

```{r 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)
```

The total is multiplied by three because the study has three treatment arms.
This rule ensures that each individual comparison and endpoint has at least
the target simulated power under its own assumptions. It is not a substitute
for a multiplicity-adjusted joint power calculation.

## Joint sample-size estimation with `k = 3`

The joint sample-size function receives all three arms, both comparison
families, and all three endpoints at once. Here, `k = 3` requires all three
endpoints to demonstrate equivalence for every comparison. Because `k = m = 3`,
all three endpoints must pass; no endpoint can be selected or omitted. An
endpoint-wise adjustment is needed when the rule allows the study to pass by
choosing only some endpoints, such as `k = 2` of `m = 3`, because there are
then several possible successful endpoint subsets. With `k = m`, there is only
one acceptable outcome: all three endpoints pass. The joint success rule
therefore already defines the required criterion, and no endpoint-wise
Bonferroni adjustment is required. The returned sample size is based on one
joint simulated success criterion rather than the maximum of separate
searches.

```{r 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
```

For a three-arm allocation, the reported total is three times the selected
number per arm. Both comparison families share the simulated test-arm counts,
and the endpoint correlation is generated through a Gaussian-copula count
model [@nelsen2006]. This is therefore a joint count simulation rather than a
maximum of separate comparison-specific searches.

## Why joint planning is important

The separate calculation above targets 80% power for each comparison-endpoint
combination. That does not mean that the complete trial has 80% probability of
passing all six requirements. For example, if two requirements each have 80%
power and are independent, their joint success probability is only
$0.80 \times 0.80 = 0.64$.

The following comparison evaluates the separate result under the actual joint
criterion. The joint calculations require all three endpoints for both
comparison families (`k = 3`). Since `k = m`, no endpoint-wise adjustment is
needed here; Bonferroni could be added only as a conservative sensitivity
analysis. The independent-endpoint scenario uses an identity correlation
matrix; the correlated scenario uses the matrix specified above.

```{r 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
```

The separate-search sample size is smaller because it guarantees only the
individual 80% targets. Its joint power can therefore be substantially below
80%, even without an endpoint-wise adjustment, whereas the joint search targets
the requested 80% trial-level power directly. This is the main advantage of joint
planning: it answers the question that matters for the confirmatory trial,
namely the probability that all required comparisons and endpoints succeed in
the same simulated study. Positive endpoint correlation can reduce the joint
sample size because endpoint successes tend to occur together, but it does not
remove the need for a joint calculation.

## Negative-binomial sensitivity analysis

If prior evidence suggests overdispersion, repeat the final joint Poisson
sample-size calculation with the negative-binomial model. All study settings
are kept the same as in `joint_result`; `update()` changes only the
distribution and dispersion. The `dispersion` parameter controls the amount
of overdispersion; larger values imply greater variability.

```{r 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
```