## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----setup, message=FALSE, warning=FALSE--------------------------------------
library(twoCoprimary)
library(dplyr)
library(tidyr)
library(knitr)

## ----basic_example------------------------------------------------------------
# Design parameters
result <- ss2Continuous(
  delta1 = 0.5,      # Effect size for endpoint 1
  delta2 = 0.5,      # Effect size for endpoint 2
  sd1 = 1,           # Standard deviation for endpoint 1
  sd2 = 1,           # Standard deviation for endpoint 2
  rho = 0.5,         # Correlation between endpoints
  r = 1,             # Balanced allocation
  alpha = 0.025,     # One-sided significance level
  beta = 0.2,        # Type II error (80% power)
  known_var = TRUE
)

print(result)

## ----correlation_impact-------------------------------------------------------
# Calculate sample sizes for different correlations
correlations <- c(0, 0.3, 0.5, 0.8)
sample_sizes <- sapply(correlations, function(rho) {
  ss2Continuous(
    delta1 = 0.5, delta2 = 0.5,
    sd1 = 1, sd2 = 1,
    rho = rho, r = 1,
    alpha = 0.025, beta = 0.2,
    known_var = TRUE
  )$N
})

# Create summary table
correlation_table <- data.frame(
  Correlation = correlations,
  Total_N = sample_sizes,
  Reduction = c(0, round((1 - sample_sizes[-1]/sample_sizes[1]) * 100, 1))
)

kable(correlation_table,
      caption = "Sample Size vs Correlation (delta = 0.5, alpha = 0.025, power = 0.8)",
      col.names = c("Correlation (rho)", "Total N", "Reduction (%)"))

## ----plot_correlation, fig.height=5, fig.width=7------------------------------
# Use plot method to visualize sample size vs correlation
plot(result, type = "sample_size_rho")

## ----plot_contour, fig.height=5, fig.width=7----------------------------------
# Create contour plot for effect sizes
plot(result, type = "effect_contour")

## ----table1_sozu--------------------------------------------------------------
# Create parameter grid (delta1 <= delta2)
param_grid <- expand.grid(
  delta1 = c(0.2, 0.25, 0.3, 0.35, 0.4),
  delta2 = c(0.2, 0.25, 0.3, 0.35, 0.4),
  sd1 = 1,
  sd2 = 1
) %>% 
  arrange(delta1, delta2) %>% 
  filter(delta2 >= delta1)

# Calculate sample sizes for different correlations
result_table <- design_table(
  param_grid = param_grid,
  rho_values = c(0, 0.3, 0.5, 0.8),
  r = 1,
  alpha = 0.025,
  beta = 0.2,
  endpoint_type = "continuous"
) %>% 
  mutate_at(vars(starts_with("rho_")), ~ . / 2)  # Per-group sample size

# Display table
kable(result_table,
      caption = "Table 1: Sample Sizes Per Group (Sozu et al. 2011, alpha = 0.025, power = 0.8)",
      digits = 2)

## ----power_calculation--------------------------------------------------------
# Calculate power with n1 = n2 = 100
power_result <- power2Continuous(
  n1 = 100, n2 = 100,
  delta1 = 0.5, delta2 = 0.5,
  sd1 = 1, sd2 = 1,
  rho = 0.5,
  alpha = 0.025,
  known_var = TRUE
)

print(power_result)

## ----power_verification-------------------------------------------------------
# Calculate sample size
ss_result <- ss2Continuous(
  delta1 = 0.5, delta2 = 0.5,
  sd1 = 1, sd2 = 1,
  rho = 0.5, r = 1,
  alpha = 0.025, beta = 0.2,
  known_var = TRUE
)

# Verify power with calculated sample size
power_check <- power2Continuous(
  n1 = ss_result$n1, n2 = ss_result$n2,
  delta1 = 0.5, delta2 = 0.5,
  sd1 = 1, sd2 = 1,
  rho = 0.5,
  alpha = 0.025,
  known_var = TRUE
)

cat("Calculated sample size per group:", ss_result$n2, "\n")
cat("Target power: 0.80\n")
cat("Achieved power:", round(power_check$powerCoprimary, 4), "\n")

## ----unified_interface--------------------------------------------------------
# Sample size calculation mode
twoCoprimary2Continuous(
  delta1 = 0.5, delta2 = 0.5,
  sd1 = 1, sd2 = 1,
  rho = 0.5, power = 0.8, r = 1,
  alpha = 0.025, known_var = TRUE
)

# Power calculation mode
twoCoprimary2Continuous(
  n1 = 100, n2 = 100,
  delta1 = 0.5, delta2 = 0.5,
  sd1 = 1, sd2 = 1,
  rho = 0.5,
  alpha = 0.025, known_var = TRUE
)

## ----unknown_variance---------------------------------------------------------
# Sample size calculation with unknown variance
ss_unknown <- ss2Continuous(
  delta1 = 0.5, delta2 = 0.5,
  sd1 = 1, sd2 = 1,
  rho = 0.5, r = 1,
  alpha = 0.025, beta = 0.2,
  known_var = FALSE,
  nMC = 10000  # Number of Monte Carlo simulations
)

print(ss_unknown)

## ----sensitivity_analysis-----------------------------------------------------
# Test robustness to correlation misspecification
assumed_rho <- 0.5
true_rhos <- c(0, 0.3, 0.5, 0.7, 0.9)

# Calculate sample size assuming rho = 0.5
ss_assumed <- ss2Continuous(
  delta1 = 0.5, delta2 = 0.5,
  sd1 = 1, sd2 = 1,
  rho = assumed_rho, r = 1,
  alpha = 0.025, beta = 0.2,
  known_var = TRUE
)

# Calculate achieved power under different true correlations
sensitivity_results <- data.frame(
  Assumed_rho = assumed_rho,
  True_rho = true_rhos,
  n_per_group = ss_assumed$n2,
  Achieved_power = sapply(true_rhos, function(true_rho) {
    power2Continuous(
      n1 = ss_assumed$n1, n2 = ss_assumed$n2,
      delta1 = 0.5, delta2 = 0.5,
      sd1 = 1, sd2 = 1,
      rho = true_rho,
      alpha = 0.025,
      known_var = TRUE
    )$powerCoprimary
  })
)

kable(sensitivity_results,
      caption = "Sensitivity Analysis: Impact of Correlation Misspecification",
      digits = 3,
      col.names = c("Assumed rho", "True rho", "n per group", "Achieved Power"))