---
title: "Verification of sample size calculation via simulation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Verification of sample size calculation via simulation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r, message=FALSE, warning=FALSE}
library(gsDesignNB)
library(data.table)
library(ggplot2)
library(gt)
```

## Introduction

This vignette verifies the accuracy of the `sample_size_nbinom()` function by comparing its theoretical predictions for average exposure,
 statistical information, and power against results from a large-scale simulation study.

We specifically test a scenario with:

*   Piecewise constant accrual rates.
*   Piecewise exponential dropout (constant in this example).
*   Negative binomial outcomes.
*   Event gaps (dead time after each event).
*   Fixed follow-up design.

The simulation uses the corrected sample size ($n = 436$) computed by `sample_size_nbinom()`, which includes the
second-order Taylor correction for Jensen's inequality in the event gap effective rate
(see `vignette("sample-size-nbinom")`).
Results for the first $n = 422$ subjects are also extracted from each simulation to compare
with the uncorrected (naive) sample size.

## Simulation design

### Parameters

*   **Rates**: $\lambda_1 = 0.4$ (Control), $\lambda_2 = 0.3$ (Experimental).
*   **Dispersion**: $k = 0.5$.
*   **Power**: 90%.
*   **Alpha**: 0.025 (one-sided).
*   **Dropout**: 10% per year adjusted to monthly rate ($\delta = 0.1 / 12$).
*   **Trial duration**: 24 months.
*   **Max follow-up**: 12 months.
*   **Event gap**: 20 days.
*   **Accrual**: Piecewise ramp-up over 12 months (Rate $R$ for 0-6mo, $2R$ for 6-12mo).

### Theoretical calculation

First, we calculate the required sample size and expected properties using `sample_size_nbinom()`.

```{r design}
# Parameters
lambda1 <- 0.4
lambda2 <- 0.3
dispersion <- 0.5
power <- 0.9
alpha <- 0.025
dropout_rate <- 0.1 / 12
max_followup <- 12
trial_duration <- 24
event_gap <- 20 / 30.42 # 20 days

# Accrual targeting 90% power
accrual_rate_rel <- c(1, 2)
accrual_duration <- c(6, 6)

design <- sample_size_nbinom(
  lambda1 = lambda1, lambda2 = lambda2, dispersion = dispersion,
  power = power,
  alpha = alpha, sided = 1,
  accrual_rate = accrual_rate_rel,
  accrual_duration = accrual_duration,
  trial_duration = trial_duration,
  dropout_rate = dropout_rate,
  max_followup = max_followup,
  event_gap = event_gap
)

print(design)
```

## Simulation results

We simulated 3,600 trials at the corrected sample size ($n = `r design$n_total`$).
For each trial, we also analyzed the first $n = 422$ subjects to provide a comparison with the 
naive (uncorrected) sample size.
This number of simulations achieves a standard error for the power estimate 
of approximately 0.005 when the true power is 90% ($\sqrt{0.9 \times 0.1 / 3600} = 0.005$). 
The simulation script is located in `data-raw/generate_simulation_data.R`.

```{r load_results}
# Load pre-computed simulation results
results_file <- system.file("extdata", "simulation_results.rds", package = "gsDesignNB")

if (results_file == "" && file.exists("../inst/extdata/simulation_results.rds")) {
  results_file <- "../inst/extdata/simulation_results.rds"
}

if (results_file != "") {
  sim_data <- readRDS(results_file)
  results <- sim_data$results
  results_naive <- sim_data$results_naive
  design_ref <- sim_data$design
  n_full <- sim_data$n_full
  n_naive <- sim_data$n_naive
} else {
  warning("Simulation results not found. Skipping verification plots.")
  results <- NULL
  results_naive <- NULL
  design_ref <- design
  n_full <- design$n_total
  n_naive <- 422L
}
```

### Power comparison

The key question is whether the corrected sample size achieves the nominal 90% power.

```{r power_comparison, eval=!is.null(results)}
power_full <- mean(results$p_value < alpha, na.rm = TRUE)
power_naive <- mean(results_naive$p_value < alpha, na.rm = TRUE)
power_full_score <- mean(results$p_value_score < alpha, na.rm = TRUE)
power_naive_score <- mean(results_naive$p_value_score < alpha, na.rm = TRUE)

ci_full <- binom.test(sum(results$p_value < alpha, na.rm = TRUE), nrow(results))$conf.int
ci_naive <- binom.test(sum(results_naive$p_value < alpha, na.rm = TRUE), nrow(results_naive))$conf.int
ci_full_score <- binom.test(sum(results$p_value_score < alpha, na.rm = TRUE), nrow(results))$conf.int
ci_naive_score <- binom.test(sum(results_naive$p_value_score < alpha, na.rm = TRUE), nrow(results_naive))$conf.int

power_df <- data.frame(
  Design = c(
    paste0("Corrected (n = ", n_full, ")"),
    paste0("Naive (n = ", n_naive, ")"),
    paste0("Corrected (n = ", n_full, ")"),
    paste0("Naive (n = ", n_naive, ")")
  ),
  Test = c("Wald", "Wald", "Score", "Score"),
  Theoretical = rep(design_ref$power, 4),
  Empirical = c(power_full, power_naive, power_full_score, power_naive_score),
  CI_Lower = c(ci_full[1], ci_naive[1], ci_full_score[1], ci_naive_score[1]),
  CI_Upper = c(ci_full[2], ci_naive[2], ci_full_score[2], ci_naive_score[2])
)

power_df |>
  gt() |>
  tab_header(
    title = md("**Power Comparison: Wald vs Score Test**"),
    subtitle = paste0("Based on ", nrow(results), " simulated trials")
  ) |>
  fmt_number(columns = c(Theoretical, Empirical, CI_Lower, CI_Upper), decimals = 4) |>
  cols_label(
    Design = "Design",
    Test = "Test",
    Theoretical = "Target",
    Empirical = "Empirical",
    CI_Lower = "95% CI Lower",
    CI_Upper = "95% CI Upper"
  )
```

### Summary of verification results (corrected design)

We compare the theoretical predictions from `sample_size_nbinom()` with the observed simulation results at the corrected sample size.

**Key distinction: Total Exposure vs Exposure at Risk**

*   **Total Exposure** (`tte_total`): The calendar time a subject is on study, from randomization 
    to the analysis cut date (or censoring). This is the same for both treatment arms by design.
*   **Exposure at Risk** (`tte`): The time during which a subject can experience a *new* event. 
    After each event, there is an "event gap" period during which new events are not counted 
    (e.g., representing recovery time or treatment effect). This differs by treatment group 
    because the group with more events loses more time to gaps.

```{r summary_table, eval=!is.null(results)}
# ---- Compute all metrics ----

# Number of simulations
n_sims <- sum(!is.na(results$estimate))

# Total Exposure (calendar follow-up time)
theo_exposure <- design_ref$exposure[1]

obs_exposure_ctrl <- mean(results$exposure_total_control)
obs_exposure_exp <- mean(results$exposure_total_experimental)
obs_exposure_at_risk_ctrl <- mean(results$exposure_at_risk_control)
obs_exposure_at_risk_exp <- mean(results$exposure_at_risk_experimental)

# Exposure at risk (time at risk excluding event gaps)
theo_exposure_at_risk_ctrl <- design_ref$exposure_at_risk_n1
theo_exposure_at_risk_exp <- design_ref$exposure_at_risk_n2

# Events by treatment group
theo_events_ctrl <- design_ref$events_n1
theo_events_exp <- design_ref$events_n2
obs_events_ctrl <- mean(results$events_control)
obs_events_exp <- mean(results$events_experimental)

# Treatment effect
true_rr <- lambda2 / lambda1
true_log_rr <- log(true_rr)
mean_log_rr <- mean(results$estimate, na.rm = TRUE)

# Variance
theo_var <- design_ref$variance
median_se_sq <- median(results$se^2, na.rm = TRUE)
emp_var <- var(results$estimate, na.rm = TRUE)

# Power
theo_power <- design_ref$power
emp_power <- mean(results$p_value < design_ref$inputs$alpha, na.rm = TRUE)

n_sim_total <- design_ref$n_total

# ---- Build summary table ----
summary_df <- data.frame(
  Metric = c(
    "Total Exposure (months) - Control",
    "Total Exposure (months) - Experimental",
    "Exposure at Risk (months) - Control",
    "Exposure at Risk (months) - Experimental",
    "Events per Subject - Control",
    "Events per Subject - Experimental",
    "Treatment Effect: log(RR)",
    "Variance of log(RR)",
    "Power"
  ),
  Theoretical = c(
    theo_exposure,
    theo_exposure,
    theo_exposure_at_risk_ctrl,
    theo_exposure_at_risk_exp,
    theo_events_ctrl / (n_sim_total / 2),
    theo_events_exp / (n_sim_total / 2),
    true_log_rr,
    theo_var,
    theo_power
  ),
  Simulated = c(
    obs_exposure_ctrl,
    obs_exposure_exp,
    obs_exposure_at_risk_ctrl,
    obs_exposure_at_risk_exp,
    obs_events_ctrl / (n_sim_total / 2),
    obs_events_exp / (n_sim_total / 2),
    mean_log_rr,
    median_se_sq,
    emp_power
  ),
  stringsAsFactors = FALSE
)

summary_df$Difference <- summary_df$Simulated - summary_df$Theoretical
summary_df$Rel_Diff_Pct <- 100 * summary_df$Difference / abs(summary_df$Theoretical)

summary_df |>
  gt() |>
  tab_header(
    title = md("**Verification of sample_size_nbinom() Predictions**"),
    subtitle = paste0("Based on ", n_sims, " simulated trials (n = ", n_sim_total, ")")
  ) |>
  tab_row_group(
    label = md("**Power**"),
    rows = Metric == "Power"
  ) |>
  tab_row_group(
    label = md("**Variance**"),
    rows = Metric == "Variance of log(RR)"
  ) |>
  tab_row_group(
    label = md("**Treatment Effect**"),
    rows = Metric == "Treatment Effect: log(RR)"
  ) |>
  tab_row_group(
    label = md("**Events**"),
    rows = grepl("Events", Metric)
  ) |>
  tab_row_group(
    label = md("**Exposure**"),
    rows = grepl("Exposure", Metric)
  ) |>
  row_group_order(groups = c("**Exposure**", "**Events**", "**Treatment Effect**",
                              "**Variance**", "**Power**")) |>
  fmt_number(columns = c(Theoretical, Simulated, Difference), decimals = 4) |>
  fmt_number(columns = Rel_Diff_Pct, decimals = 2) |>
  cols_label(
    Metric = "",
    Theoretical = "Theoretical",
    Simulated = "Simulated",
    Difference = "Difference",
    Rel_Diff_Pct = "Rel. Diff (%)"
  ) |>
  sub_missing(missing_text = "—")
```

**Notes:**

*   **Total Exposure** is the average calendar follow-up time per subject. The theoretical value is 
    the same for both arms (symmetric design), while simulated values may differ slightly due to 
    sampling variability in enrollment and dropout.
*   **Exposure at Risk** is the time during which new events can occur, after subtracting the 
    "event gap" periods following each event. The control group has lower exposure at risk because 
    it experiences more events (higher rate), each followed by a gap period.
*   **Events per Subject** is the mean number of events per subject in each treatment group.
*   **Variance** uses the median of the estimated variances (SE²) from each simulation, which is robust 
    to outliers from unstable model fits.

### Distribution of the test statistic

Here we plot the density of the Wald Z-scores from the simulation and compare it to the expected normal distribution centered at the theoretical mean Z-score.

```{r z_score_dist, eval=!is.null(results)}
z_scores <- results$estimate / results$se

theo_se <- sqrt(theo_var)
theo_mean_z <- log(lambda2 / lambda1) / theo_se
crit_val <- qnorm(design_ref$inputs$alpha)
prop_reject <- mean(z_scores < crit_val, na.rm = TRUE)

ggplot(data.frame(z = z_scores), aes(x = z)) +
  geom_density(aes(color = "Simulated Density"), linewidth = 1) +
  stat_function(
    fun = dnorm,
    args = list(mean = theo_mean_z, sd = 1),
    aes(color = "Theoretical Normal"),
    linewidth = 1,
    linetype = "dashed"
  ) +
  geom_vline(xintercept = crit_val, linetype = "dotted", color = "black") +
  annotate(
    "text",
    x = crit_val, y = 0.05,
    label = paste0("  Reject: ", round(prop_reject * 100, 1), "%"),
    hjust = 0, vjust = 0
  ) +
  labs(
    title = "Distribution of Wald Z-scores (corrected design)",
    subtitle = paste("Theoretical Mean Z =", round(theo_mean_z, 3)),
    x = "Z-score (Estimate / Estimated SE)",
    y = "Density",
    color = "Distribution"
  ) +
  theme_minimal() +
  theme(legend.position = "top")
```

### Wald vs score Z-score comparison

```{r z_score_comparison, eval=!is.null(results) && "z_score" %in% names(results)}
z_wald <- results$estimate / results$se
z_score <- results$z_score

z_df <- data.frame(
  z = c(z_wald, z_score),
  Test = rep(c("Wald", "Score"), each = length(z_wald))
)

ggplot(z_df, aes(x = z, color = Test)) +
  geom_density(linewidth = 1) +
  stat_function(
    fun = dnorm,
    args = list(mean = theo_mean_z, sd = 1),
    linewidth = 0.8,
    linetype = "dashed",
    color = "grey40"
  ) +
  geom_vline(xintercept = crit_val, linetype = "dotted", color = "black") +
  labs(
    title = "Z-score Distribution: Wald vs Score (corrected design)",
    subtitle = paste0("Dashed = theoretical N(", round(theo_mean_z, 2), ", 1)"),
    x = "Z-score",
    y = "Density",
    color = "Test"
  ) +
  theme_minimal() +
  theme(legend.position = "top")
```

### Distribution statistics for log(RR)

```{r log_rr_stats, eval=!is.null(results)}
theo_mean_log_rr <- log(lambda2 / lambda1)
theo_sd_log_rr <- sqrt(theo_var)
emp_mean_log_rr <- mean(results$estimate, na.rm = TRUE)
emp_sd_log_rr <- sd(results$estimate, na.rm = TRUE)
emp_median_log_rr <- median(results$estimate, na.rm = TRUE)

ests <- results$estimate[!is.na(results$estimate)]
q_low <- quantile(ests, 0.01)
q_high <- quantile(ests, 0.99)
ests_trimmed <- ests[ests >= q_low & ests <= q_high]
n_trimmed <- length(ests_trimmed)
m3 <- sum((ests_trimmed - mean(ests_trimmed))^3) / n_trimmed
m4 <- sum((ests_trimmed - mean(ests_trimmed))^4) / n_trimmed
s2 <- sum((ests_trimmed - mean(ests_trimmed))^2) / n_trimmed
emp_skew_log_rr <- m3 / s2^(3/2)
emp_kurt_log_rr <- m4 / s2^2

comparison_log_rr <- data.frame(
  Metric = c("Mean", "SD", "Median", "Skewness (trimmed)", "Kurtosis (trimmed)"),
  Theoretical = c(theo_mean_log_rr, theo_sd_log_rr, theo_mean_log_rr, 0, 3),
  Simulated = c(emp_mean_log_rr, emp_sd_log_rr, emp_median_log_rr, emp_skew_log_rr, emp_kurt_log_rr),
  Difference = c(
    emp_mean_log_rr - theo_mean_log_rr,
    emp_sd_log_rr - theo_sd_log_rr,
    emp_median_log_rr - theo_mean_log_rr,
    emp_skew_log_rr - 0,
    emp_kurt_log_rr - 3
  )
)
comparison_log_rr |>
  gt() |>
  tab_header(title = md("**Comparison of log(RR) Statistics**")) |>
  fmt_number(columns = where(is.numeric), decimals = 4)
```

## Type I error: Wald vs score test

A well-calibrated test should reject at most $\alpha = 0.025$ under the null ($\text{RR} = 1$).
We simulated trials under the null hypothesis using the same design parameters
(accrual, dropout, dispersion, sample size) but with $\lambda_2 = \lambda_1$ in both arms.
The Wald test (based on the ML estimate under the alternative) tends to be slightly
anti-conservative in finite samples, while the score test (evaluated entirely under $H_0$)
is expected to be better calibrated.

```{r load_null}
null_file <- system.file("extdata", "null_simulation_results.rds", package = "gsDesignNB")
if (null_file == "" && file.exists("../inst/extdata/null_simulation_results.rds")) {
  null_file <- "../inst/extdata/null_simulation_results.rds"
}
has_null <- null_file != ""
if (has_null) {
  null_data <- readRDS(null_file)
  results_null <- null_data$results_null
}
```

```{r type1_table, eval=has_null}
type1_wald <- mean(results_null$p_value_wald < alpha, na.rm = TRUE)
type1_score <- mean(results_null$p_value_score < alpha, na.rm = TRUE)

n_null <- sum(!is.na(results_null$p_value_wald))
ci_wald <- binom.test(sum(results_null$p_value_wald < alpha, na.rm = TRUE), n_null)$conf.int
ci_score <- binom.test(sum(results_null$p_value_score < alpha, na.rm = TRUE), n_null)$conf.int

type1_df <- data.frame(
  Test = c("Wald (ML)", "Score (null model)"),
  Nominal = c(alpha, alpha),
  Empirical = c(type1_wald, type1_score),
  CI_Lower = c(ci_wald[1], ci_score[1]),
  CI_Upper = c(ci_wald[2], ci_score[2])
)

type1_df |>
  gt() |>
  tab_header(
    title = md("**Type I Error: Wald vs Score Test Under RR = 1**"),
    subtitle = paste0("Based on ", n_null, " simulated null trials (n = ", null_data$n_total, ")")
  ) |>
  fmt_number(columns = c(Nominal, Empirical, CI_Lower, CI_Upper), decimals = 4) |>
  cols_label(
    Test = "Test",
    Nominal = "Nominal α",
    Empirical = "Empirical",
    CI_Lower = "95% CI Lower",
    CI_Upper = "95% CI Upper"
  )
```

```{r type1_z_density, eval=has_null}
z_null_df <- data.frame(
  z = c(results_null$z_wald, results_null$z_score),
  Test = rep(c("Wald", "Score"), each = nrow(results_null))
)

ggplot(z_null_df, aes(x = z, color = Test)) +
  geom_density(linewidth = 1) +
  stat_function(
    fun = dnorm, args = list(mean = 0, sd = 1),
    linewidth = 0.8, linetype = "dashed", color = "grey40"
  ) +
  geom_vline(xintercept = qnorm(alpha), linetype = "dotted", color = "black") +
  annotate("text", x = qnorm(alpha), y = 0.02,
           label = paste0("  α = ", alpha), hjust = 0) +
  labs(
    title = "Null Z-score Distribution: Wald vs Score",
    subtitle = "Dashed = N(0, 1) reference",
    x = "Z-score",
    y = "Density",
    color = "Test"
  ) +
  theme_minimal() +
  theme(legend.position = "top")
```

## Scenario sweep: impact of the Jensen correction

The Jensen correction for event gaps only applies when both the dispersion ($k > 0$) and the event gap ($g > 0$) are nonzero.
When either is zero, the corrected and naive formulas produce identical sample sizes.

To assess the correction's impact across a range of realistic parameter combinations, we ran a simulation sweep with 
10,000 replicates per scenario. For each scenario, trials were simulated at the corrected sample size ($n_\text{corr}$), 
and results were also extracted for the first $n_\text{naive}$ subjects (the uncorrected sample size). 
Because both analyses use the same simulated data, the **paired** comparison is very precise.

```{r load_sweep}
sweep_file <- system.file("extdata", "scenario_sweep_results.rds", package = "gsDesignNB")
if (sweep_file == "" && file.exists("../inst/extdata/scenario_sweep_results.rds")) {
  sweep_file <- "../inst/extdata/scenario_sweep_results.rds"
}
has_sweep <- sweep_file != ""
if (has_sweep) {
  sweep <- readRDS(sweep_file)
}
```

```{r sweep_table, eval=has_sweep}
sweep_display <- data.frame(
  Scenario = sprintf("k = %.1f, gap = %d days", sweep$k, sweep$gap_days),
  `n (corrected)` = sweep$n_corrected,
  `n (naive)` = sweep$n_naive,
  `Extra subjects` = sweep$n_diff,
  `Power (corrected)` = sweep$power_corrected,
  `Power (naive)` = sweep$power_naive,
  `Difference` = sprintf("%.2f ± %.2f", 100 * sweep$power_diff, 100 * 1.96 * sweep$power_diff_se),
  check.names = FALSE
)

sweep_display |>
  gt() |>
  tab_header(
    title = md("**Impact of Jensen Correction Across Scenarios**"),
    subtitle = "10,000 replicates per scenario; power difference is paired (95% CI)"
  ) |>
  fmt_number(columns = c(`Power (corrected)`, `Power (naive)`), decimals = 4) |>
  cols_label(
    Scenario = "Scenario",
    `n (corrected)` = "n (corr.)",
    `n (naive)` = "n (naive)",
    `Extra subjects` = "Δn",
    `Power (corrected)` = "Power (corr.)",
    `Power (naive)` = "Power (naive)",
    Difference = "Diff (pp) ± 95% CI"
  )
```

```{r sweep_interpretation, eval=has_sweep, results='asis'}
cat("**Key findings:**\n\n")
cat("*   The corrected design achieves power at or above the 90% target in **all** scenarios.\n")
cat("*   The naive design falls below 90% in **three of four** scenarios, with the shortfall increasing as $k$ and $g$ grow.\n")
cat(sprintf("*   The largest paired power difference is %.1f percentage points ($k = %.1f$, gap = %d days), ",
    100 * sweep$power_diff[which.max(sweep$power_diff)],
    sweep$k[which.max(sweep$power_diff)],
    sweep$gap_days[which.max(sweep$power_diff)]))
cat(sprintf("where the corrected design adds %d subjects.\n",
    sweep$n_diff[which.max(sweep$power_diff)]))
cat(sprintf("*   All paired differences are statistically significant (discordant pair ratios range from %.1f:1 to %.1f:1 in favor of the corrected design).\n",
    min(sweep$corr_only / sweep$naive_only),
    max(sweep$corr_only / sweep$naive_only)))
```

## Discussion: why we apply the correction

In the primary verification scenario above ($k = 0.5$, gap = 20 days), the corrected design slightly overshoots
the 90% target while the naive design is very close to 90%.
This pattern --- where the naive formula appears well-calibrated --- reflects a **partial cancellation of two biases**:

1.  **Design-stage bias** (hurts power): The naive effective rate $\lambda/(1+\lambda g)$ overestimates the true
    population-level effective rate, leading to an underestimated variance and thus a sample size that is too small.
    
2.  **Analysis-stage bias** (helps power): The model-based standard error from `MASS::glm.nb()` tends to
    slightly underestimate the true standard deviation of $\hat\theta$, making the Wald test slightly
    anti-conservative (z-statistic too extreme, more rejections).

These biases act in opposite directions.
In the $k = 0.5$, gap = 20 day scenario, they happen to nearly cancel, giving the naive design the appearance of correct calibration.

We apply the Jensen correction despite this for several reasons:

*   **The cancellation is fragile.** It depends on the specific test statistic (NB Wald via `glm.nb`). A different analysis method --- sandwich standard errors, bootstrap, permutation test, or a different link function --- would change the analysis-stage bias while the design-stage bias remains, breaking the cancellation.

*   **The correction is conservative.** It produces power at or slightly above 90% rather than at or slightly below.  For a pivotal trial, a slight overshoot of the power target is preferable to a slight undershoot.

*   **The correction is principled.** It fixes the effective rate approximation by accounting for the known distributional assumption (Gamma frailty). The remaining small overshoot is a separate finite-sample phenomenon.

*   **The scenario sweep demonstrates the benefit.** As the table above shows, the naive formula increasingly underpowers as $k$ and $g$ grow, while the corrected formula maintains nominal or conservative power across all scenarios tested.

*   **The cost is modest.** The extra subjects required by the correction range from 14 to 40 (3--6% of total), a small price for reliable power.
