---
title: "Validation of bbssr Package Functions"
author: "Gosuke Homma"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Validation of bbssr Package Functions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  fig.align = "center",
  warning = FALSE,
  message = FALSE
)

# Load required packages
# Note: These packages must be installed separately before building this vignette
library(bbssr)

# Check if optional packages are available
required_packages <- c("dplyr", "ggplot2", "tidyr", "microbenchmark")
for (pkg in required_packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    message(paste("Package", pkg, "not available. Install with: install.packages('", pkg, "')", sep = ""))
  } else {
    library(pkg, character.only = TRUE)
  }
}

# Load base packages explicitly
library(stats)
library(utils)
library(base)

# Check if validation packages are available
has_exact <- requireNamespace("Exact", quietly = TRUE)
has_exact2x2 <- requireNamespace("exact2x2", quietly = TRUE)

if (!has_exact) {
  message("Note: Exact package not available. Install with: install.packages('Exact')")
} else {
  library(Exact)
}

if (!has_exact2x2) {
  message("Note: exact2x2 package not available. Install with: install.packages('exact2x2')")
} else {
  library(exact2x2)
}

if (!has_exact || !has_exact2x2) {
  knitr::opts_chunk$set(eval = FALSE)
  warning("Validation packages not available. Code will not be evaluated.")
}
```

## Introduction

This vignette validates the `BinaryPower()` and `BinarySampleSize()` functions in the `bbssr` package by comparing results with established R packages: `Exact` and `exact2x2`. We also compare computational performance to demonstrate the efficiency improvements in `bbssr`.

The validation covers all five statistical tests implemented in `bbssr`:
1. Pearson chi-squared test (`'Chisq'`)
2. Fisher exact test (`'Fisher'`)
3. Fisher mid-p test (`'Fisher-midP'`)
4. Z-pooled exact unconditional test (`'Z-pool'`)
5. Boschloo exact unconditional test (`'Boschloo'`)

```{r load_packages, eval=has_exact && has_exact2x2}
# Load external validation packages
library(Exact)
library(exact2x2)
```

```{r load_bbssr_only, eval=!has_exact || !has_exact2x2}
# If external packages are not available, just load bbssr
library(bbssr)
```

## Validation of BinaryPower Function

### Test Parameters

We use the following parameters for validation:

```{r validation_parameters}
# Test parameters for validation
p1 <- c(0.5, 0.6, 0.7, 0.8)  # Response rates in treatment group
p2 <- c(0.2, 0.2, 0.2, 0.2)  # Response rates in control group
N1 <- 10                     # Sample size in treatment group
N2 <- 40                     # Sample size in control group
alpha <- 0.025              # One-sided significance level
```

## Quick Verification Example

For users who want to see immediate results without running external packages:

```{r quick_verification}
# Quick verification using bbssr only
cat("=== Quick Verification Results ===\n")
cat("Test Parameters:\n")
cat("p1 =", paste(p1, collapse = ", "), "\n")
cat("p2 =", paste(p2, collapse = ", "), "\n") 
cat("N1 =", N1, ", N2 =", N2, ", alpha =", alpha, "\n\n")

# Show all test results
tests <- c('Chisq', 'Fisher', 'Fisher-midP', 'Z-pool', 'Boschloo')
for(test in tests) {
  powers <- BinaryPower(p1, p2, N1, N2, alpha, Test = test)
  cat(sprintf("%-12s: %s\n", test, paste(round(powers, 6), collapse = "  ")))
}
```

### 1. Pearson Chi-squared Test

```{r chisq_validation, eval=has_exact && has_exact2x2}
# bbssr results
bbssr_chisq <- BinaryPower(p1, p2, N1, N2, alpha, Test = 'Chisq')
print("bbssr results:")
print(round(bbssr_chisq, 7))

# Exact package results for comparison
exact_chisq <- sapply(seq(p1), function(i) {
  Exact::power.exact.test(p1[i], p2[i], N1, N2, method = 'pearson chisq', 'greater', alpha)$power
})
print("Exact package results:")
print(round(exact_chisq, 7))

# Create comparison table
chisq_comparison <- data.frame(
  p1 = p1,
  p2 = p2,
  bbssr = round(bbssr_chisq, 7),
  Exact = round(exact_chisq, 7),
  Difference = round(abs(bbssr_chisq - exact_chisq), 10)
)

knitr::kable(chisq_comparison, caption = "Pearson Chi-squared Test: bbssr vs Exact Package")
```

### 2. Fisher Exact Test

```{r fisher_validation, eval=has_exact && has_exact2x2}
# bbssr results
bbssr_fisher <- BinaryPower(p1, p2, N1, N2, alpha, Test = 'Fisher')

# Exact package results for comparison
exact_fisher <- sapply(seq(p1), function(i) {
  Exact::power.exact.test(p1[i], p2[i], N1, N2, method = 'fisher', 'greater', alpha)$power
})

# Create comparison table
fisher_comparison <- data.frame(
  p1 = p1,
  p2 = p2,
  bbssr = round(bbssr_fisher, 7),
  Exact = round(exact_fisher, 7),
  Difference = round(abs(bbssr_fisher - exact_fisher), 10)
)

knitr::kable(fisher_comparison, caption = "Fisher Exact Test: bbssr vs Exact Package")
```

### 3. Fisher Mid-p Test

```{r fisher_midp_validation, eval=has_exact && has_exact2x2}
# bbssr results
bbssr_midp <- BinaryPower(p1, p2, N1, N2, alpha, Test = 'Fisher-midP')

# exact2x2 package results for comparison
exact2x2_midp <- sapply(seq(p1), function(i) {
  exact2x2::Power2x2(N1, N2, p1[i], p2[i], alpha, pvalFunc = 
             function(x1, n1, x2, n2) {
               fisher.exact(matrix(c(x1, n1 - x1, x2, n2 - x2), 2), 
                          alt = 'greater', midp = TRUE)$p.value
             }
  )
})

# Create comparison table
midp_comparison <- data.frame(
  p1 = p1,
  p2 = p2,
  bbssr = round(bbssr_midp, 7),
  exact2x2 = round(exact2x2_midp, 7),
  Difference = round(abs(bbssr_midp - exact2x2_midp), 10)
)

knitr::kable(midp_comparison, caption = "Fisher Mid-p Test: bbssr vs exact2x2 Package")
```

### 4. Z-pooled Exact Unconditional Test

```{r zpool_validation, eval=has_exact && has_exact2x2}
# bbssr results
bbssr_zpool <- BinaryPower(p1, p2, N1, N2, alpha, Test = 'Z-pool')

# Exact package results for comparison
exact_zpool <- sapply(seq(p1), function(i) {
  Exact::power.exact.test(p1[i], p2[i], N1, N2, method = 'z-pooled', 'greater', alpha)$power
})

# Create comparison table
zpool_comparison <- data.frame(
  p1 = p1,
  p2 = p2,
  bbssr = round(bbssr_zpool, 7),
  Exact = round(exact_zpool, 7),
  Difference = round(abs(bbssr_zpool - exact_zpool), 10)
)

knitr::kable(zpool_comparison, caption = "Z-pooled Test: bbssr vs Exact Package")
```

### 5. Boschloo Exact Unconditional Test

```{r boschloo_validation, eval=has_exact && has_exact2x2}
# bbssr results
bbssr_boschloo <- BinaryPower(p1, p2, N1, N2, alpha, Test = 'Boschloo')

# Exact package results for comparison
exact_boschloo <- sapply(seq(p1), function(i) {
  Exact::power.exact.test(p1[i], p2[i], N1, N2, method = 'boschloo', 'greater', alpha)$power
})

# Create comparison table
boschloo_comparison <- data.frame(
  p1 = p1,
  p2 = p2,
  bbssr = round(bbssr_boschloo, 7),
  Exact = round(exact_boschloo, 7),
  Difference = round(abs(bbssr_boschloo - exact_boschloo), 10)
)

knitr::kable(boschloo_comparison, caption = "Boschloo Test: bbssr vs Exact Package")
```

## Performance Comparison

### Enhanced Computational Speed Benchmarks

We compare the computational speed using more demanding scenarios to highlight bbssr's performance advantages.

```{r performance_setup}
# Enhanced parameters for performance comparison
perf_scenarios <- expand.grid(
  p1 = c(0.5, 0.6, 0.7, 0.8),
  p2 = c(0.2, 0.3, 0.4),
  N1 = c(15, 20, 25),
  N2 = c(15, 20, 25)
) %>%
  filter(p1 > p2) %>%  # Only meaningful comparisons
  slice_head(n = 8)    # Select representative scenarios

perf_alpha <- 0.025

# Display scenarios for transparency
knitr::kable(perf_scenarios, caption = "Performance Benchmark Scenarios")
```

### BinaryPower Performance Comparison

```{r power_performance, eval=has_exact && has_exact2x2}
# Function to benchmark a single scenario
benchmark_scenario <- function(p1, p2, N1, N2) {
  microbenchmark::microbenchmark(
    bbssr_chisq = BinaryPower(p1, p2, N1, N2, perf_alpha, Test = 'Chisq'),
    exact_chisq = Exact::power.exact.test(p1, p2, N1, N2, 
                                         method = 'pearson chisq', 'greater', perf_alpha),
    
    bbssr_fisher = BinaryPower(p1, p2, N1, N2, perf_alpha, Test = 'Fisher'),
    exact_fisher = Exact::power.exact.test(p1, p2, N1, N2, 
                                          method = 'fisher', 'greater', perf_alpha),
    
    bbssr_zpool = BinaryPower(p1, p2, N1, N2, perf_alpha, Test = 'Z-pool'),
    exact_zpool = Exact::power.exact.test(p1, p2, N1, N2, 
                                         method = 'z-pooled', 'greater', perf_alpha),
    
    bbssr_boschloo = BinaryPower(p1, p2, N1, N2, perf_alpha, Test = 'Boschloo'),
    exact_boschloo = Exact::power.exact.test(p1, p2, N1, N2, 
                                            method = 'boschloo', 'greater', perf_alpha),
    
    times = 20
  )
}

# Run benchmarks for multiple scenarios
benchmark_results <- list()
for(i in 1:nrow(perf_scenarios)) {
  scenario <- perf_scenarios[i, ]
  cat(sprintf("Benchmarking scenario %d: p1=%.1f, p2=%.1f, N1=%d, N2=%d\n", 
              i, scenario$p1, scenario$p2, scenario$N1, scenario$N2))
  
  benchmark_results[[i]] <- benchmark_scenario(
    scenario$p1, scenario$p2, scenario$N1, scenario$N2
  ) %>%
    summary() %>%
    mutate(
      scenario_id = i,
      p1 = scenario$p1,
      p2 = scenario$p2,
      N1 = scenario$N1,
      N2 = scenario$N2,
      total_n = scenario$N1 + scenario$N2
    )
}

# Combine all benchmark results
all_benchmarks <- do.call(rbind, benchmark_results)

# Calculate speed improvements
speed_comparison <- all_benchmarks %>%
  mutate(
    test_name = case_when(
      grepl("chisq", expr) ~ "Chi-squared",
      grepl("fisher", expr) & !grepl("zpool|boschloo", expr) ~ "Fisher",
      grepl("zpool", expr) ~ "Z-pooled",
      grepl("boschloo", expr) ~ "Boschloo"
    ),
    package = if_else(grepl("bbssr", expr), "bbssr", "Exact")
  ) %>%
  group_by(scenario_id, test_name, package, p1, p2, N1, N2, total_n) %>%
  summarise(median_time = median(median), .groups = 'drop') %>%
  tidyr::pivot_wider(names_from = package, values_from = median_time) %>%
  mutate(
    speed_improvement = round(Exact / bbssr, 1),
    bbssr_ms = round(bbssr / 1000, 2),
    exact_ms = round(Exact / 1000, 2)
  )

# Display summary
speed_summary <- speed_comparison %>%
  group_by(test_name) %>%
  summarise(
    avg_improvement = round(mean(speed_improvement), 1),
    max_improvement = round(max(speed_improvement), 1),
    min_improvement = round(min(speed_improvement), 1),
    .groups = 'drop'
  )

knitr::kable(speed_summary, 
             col.names = c("Test", "Avg Speed-up", "Max Speed-up", "Min Speed-up"),
             caption = "Speed Improvement Summary: bbssr vs Exact Package (times faster)")
```

### Performance Visualization

```{r performance_plot, eval=has_exact && has_exact2x2, fig.width=8, fig.height=10, out.width="100%"}
# Create speed improvement comparison (how many times faster bbssr is)
speed_plot_data <- speed_comparison %>%
  mutate(
    test_name = factor(test_name, levels = c("Chi-squared", "Fisher", "Z-pooled", "Boschloo")),
    scenario_label = paste0("Scenario ", scenario_id, "\n(N=", total_n, ")")
  )

# Speed improvement plot
speed_plot <- ggplot(speed_plot_data, aes(x = scenario_label, y = speed_improvement)) +
  geom_col(fill = "#4CAF50", alpha = 0.8, width = 0.6) +
  geom_text(aes(label = paste0(speed_improvement, "x")), 
            vjust = -0.3, size = 3, fontface = "bold") +
  facet_wrap(~test_name, ncol = 1, scales = "free_y") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 1) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +  # Add 15% padding at top
  labs(
    title = "Speed Improvement: bbssr vs Exact Package",
    subtitle = "How many times faster bbssr is compared to Exact package",
    x = "Benchmark Scenario",
    y = "Speed Improvement Factor (times faster)",
    caption = "Red dashed line shows no improvement (factor = 1). Higher bars = better performance."
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.5, margin = margin(b = 5)),
    plot.subtitle = element_text(size = 11, hjust = 0.5, margin = margin(b = 15)),
    strip.text = element_text(size = 12, face = "bold", margin = margin(t = 8, b = 8)),
    strip.background = element_rect(fill = "gray95", color = "gray80"),
    axis.title.x = element_text(size = 11, margin = margin(t = 10)),
    axis.title.y = element_text(size = 11, margin = margin(r = 10)),
    axis.text.x = element_text(size = 8, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 9),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.margin = margin(t = 10, r = 10, b = 15, l = 10)
  )

print(speed_plot)

# Detailed execution time comparison
time_comparison_data <- speed_comparison %>%
  select(scenario_id, test_name, bbssr_ms, exact_ms, total_n) %>%
  tidyr::pivot_longer(
    cols = c(bbssr_ms, exact_ms),
    names_to = "package",
    values_to = "time_ms"
  ) %>%
  mutate(
    package = case_when(
      package == "bbssr_ms" ~ "bbssr",
      package == "exact_ms" ~ "Exact"
    ),
    test_name = factor(test_name, levels = c("Chi-squared", "Fisher", "Z-pooled", "Boschloo")),
    scenario_label = paste0("Scenario ", scenario_id, "\n(N=", total_n, ")")
  )

# Execution time plot
time_plot <- ggplot(time_comparison_data, aes(x = scenario_label, y = time_ms, fill = package)) +
  geom_col(position = "dodge", width = 0.7) +
  facet_wrap(~test_name, ncol = 1, scales = "free_y") +
  scale_fill_manual(
    values = c("bbssr" = "#2E8B57", "Exact" = "#CD5C5C"),
    name = "Package"
  ) +
  labs(
    title = "Execution Time Comparison: bbssr vs Exact Package",
    subtitle = "Actual execution times in milliseconds",
    x = "Benchmark Scenario",
    y = "Execution Time (milliseconds)",
    caption = "Lower bars indicate better performance"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.5, margin = margin(b = 5)),
    plot.subtitle = element_text(size = 11, hjust = 0.5, margin = margin(b = 15)),
    strip.text = element_text(size = 12, face = "bold", margin = margin(t = 8, b = 8)),
    strip.background = element_rect(fill = "gray95", color = "gray80"),
    legend.position = "bottom",
    legend.title = element_text(size = 11),
    legend.text = element_text(size = 10),
    axis.title.x = element_text(size = 11, margin = margin(t = 10)),
    axis.title.y = element_text(size = 11, margin = margin(r = 10)),
    axis.text.x = element_text(size = 8, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 9),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.margin = margin(t = 10, r = 10, b = 15, l = 10)
  )

print(time_plot)
```

### Speed Improvement Summary

```{r speed_improvement_summary, eval=has_exact && has_exact2x2}
# Overall performance summary
overall_summary <- speed_comparison %>%
  summarise(
    total_scenarios = n(),
    overall_avg_improvement = round(mean(speed_improvement), 1),
    overall_max_improvement = round(max(speed_improvement), 1),
    overall_min_improvement = round(min(speed_improvement), 1),
    scenarios_over_2x = sum(speed_improvement >= 2),
    scenarios_over_5x = sum(speed_improvement >= 5),
    scenarios_over_10x = sum(speed_improvement >= 10)
  )

cat("Performance Summary:\n")
cat(sprintf("- Total benchmark scenarios: %d\n", overall_summary$total_scenarios))
cat(sprintf("- Average speed improvement: %.1fx faster\n", overall_summary$overall_avg_improvement))
cat(sprintf("- Maximum speed improvement: %.1fx faster\n", overall_summary$overall_max_improvement))
cat(sprintf("- Minimum speed improvement: %.1fx faster\n", overall_summary$overall_min_improvement))
cat(sprintf("- Scenarios where bbssr is >2x faster: %d (%.1f%%)\n", 
            overall_summary$scenarios_over_2x, 
            100 * overall_summary$scenarios_over_2x / overall_summary$total_scenarios))
cat(sprintf("- Scenarios where bbssr is >5x faster: %d (%.1f%%)\n", 
            overall_summary$scenarios_over_5x, 
            100 * overall_summary$scenarios_over_5x / overall_summary$total_scenarios))
cat(sprintf("- Scenarios where bbssr is >10x faster: %d (%.1f%%)\n", 
            overall_summary$scenarios_over_10x, 
            100 * overall_summary$scenarios_over_10x / overall_summary$total_scenarios))
```

## Summary of Validation Results

### Accuracy Validation

The validation results demonstrate that `bbssr` functions produce **identical results** to established packages:

- **BinaryPower()**: All five statistical tests show perfect agreement with `Exact` and `exact2x2` packages
- **Numerical Precision**: Differences are within machine precision (< 1e-10)

### Performance Improvements

The `bbssr` package demonstrates significant computational efficiency across all statistical tests.

### Key Advantages of bbssr

1. **Computational Efficiency**: Significantly faster execution times across all statistical tests
2. **Numerical Accuracy**: Identical results to established packages with robust implementation
3. **Comprehensive Functionality**: Unified interface for multiple exact tests and BSSR methodology
4. **Optimized Implementation**: Efficient algorithms designed for clinical trial applications

## Conclusion

The validation results confirm that `bbssr` provides:

- **Accurate calculations** identical to established packages
- **Superior performance** with substantial speed improvements
- **Reliable implementation** suitable for production use in clinical trials

The package successfully combines statistical rigor with computational efficiency, making it an excellent choice for implementing blinded sample size re-estimation in clinical trials with binary endpoints.

## Session Information

```{r session_info}
sessionInfo()
```