## ----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)

## ----table2_replication-------------------------------------------------------
# Recreate Sozu et al. (2012) Table 2
library(dplyr)
library(tidyr)

param_grid_mixed_cb_ss <- expand.grid(
  delta = 4.4,
  sd = c(19, 20, 21, 22),
  p1 = 0.59,
  p2 = 0.46
)

result_mixed_cb_ss <- design_table(
  param_grid = param_grid_mixed_cb_ss,
  rho_values = c(0, 0.3, 0.5, 0.8),
  r = 1,
  alpha = 0.025,
  beta = 0.2,
  endpoint_type = "mixed_cont_binary",
  Test = "AN"
) %>% 
  mutate_at(vars(starts_with("rho_")), ~ . / 2)

kable(result_mixed_cb_ss,
      caption = "Table 2: Sample Size per Group (n) for PREMIER Study Scenario (delta1 = 4.4, p1 = 0.59, p2 = 0.46, α = 0.025, 1-β = 0.80)",
      digits = 2,
      col.names = c("delta", "σ", "p1", "p2", "ρ=0.0", "ρ=0.3", "ρ=0.5", "ρ=0.8"))

## ----table5_replication-------------------------------------------------------
# Recreate Supporting Information Table 5
param_grid_mixed_cb_ss2 <- tibble(
  delta = c(0.235, 0.397, 0.521, 0.190, 0.335, 0.457),
  sd = 1,
  p1 = c(rep(0.99, 3), rep(0.95, 3)),
  p2 = c(seq(0.95, 0.85, length.out = 3), seq(0.90, 0.80, length.out = 3))
)

result_mixed_cb_ss2 <- do.call(
  bind_rows,
  lapply(c("ANc", "ASc"), function(test) {
    design_table(
      param_grid = param_grid_mixed_cb_ss2,
      rho_values = c(0, 0.3, 0.5, 0.8),
      r = 1,
      alpha = 0.025,
      beta = 0.2,
      endpoint_type = "mixed_cont_binary",
      Test = test
    ) %>% 
      mutate_at(vars(starts_with("rho_")), ~ . / 2) %>% 
      mutate(Test = test)
  })
) %>% 
  arrange(desc(p1), delta) %>% 
  select(delta, sd, p1, p2, Test, everything())

# Display for ANc
result_anc <- result_mixed_cb_ss2 %>% 
  filter(Test == "ANc") %>% 
  select(-Test)

kable(result_anc,
      caption = "Table 5 (Part A): Sample Size per Group (n) with Continuity Correction (ANc) (σ = 1, α = 0.025, 1-β = 0.80)^a^",
      digits = 3,
      col.names = c("delta", "σ", "p1", "p2", "ρ=0.0", "ρ=0.3", "ρ=0.5", "ρ=0.8"))

# Display for ASc
result_asc <- result_mixed_cb_ss2 %>% 
  filter(Test == "ASc") %>% 
  select(-Test)

kable(result_asc,
      caption = "Table 5 (Part B): Sample Size per Group (n) with Arcsine and Continuity Correction (ASc) (σ = 1, α = 0.025, 1-β = 0.80)^a^",
      digits = 3,
      col.names = c("delta", "σ", "p1", "p2", "ρ=0.0", "ρ=0.3", "ρ=0.5", "ρ=0.8"))

## ----example1-----------------------------------------------------------------
# Balanced design: nT = nC (i.e., r = 1, which corresponds to kappa = 1)
result_balanced <- ss2MixedContinuousBinary(
  delta = 0.5,           # Standardized effect for continuous endpoint
  sd = 1,                # Standard deviation
  p1 = 0.7,              # Success prob in treatment group
  p2 = 0.5,              # Success prob in control group
  rho = 0.5,             # Biserial correlation
  r = 1,                 # Balanced allocation (r = nT/nC = 1)
  alpha = 0.025,
  beta = 0.2,
  Test = "AN"
)

print(result_balanced)

## ----example2-----------------------------------------------------------------
# Fixed effect sizes
delta <- 0.5
p1 <- 0.7
p2 <- 0.5

# Range of correlations
rho_values <- c(0, 0.3, 0.5, 0.8)

ss_by_rho <- sapply(rho_values, function(rho) {
  result <- ss2MixedContinuousBinary(
    delta = delta,
    sd = 1,
    p1 = p1,
    p2 = p2,
    rho = rho,
    r = 1,
    alpha = 0.025,
    beta = 0.2,
    Test = "AN"
  )
  result$n2
})

result_df <- data.frame(
  rho = rho_values,
  n_per_group = ss_by_rho,
  N_total = ss_by_rho * 2,
  reduction_pct = round((1 - ss_by_rho / ss_by_rho[1]) * 100, 1)
)

kable(result_df,
      caption = "Effect of Biserial Correlation on Sample Size",
      col.names = c("ρ", "n per group", "N total", "Reduction (%)"))

# Plot
plot(rho_values, ss_by_rho, 
     type = "b", pch = 19,
     xlab = "Biserial Correlation (ρ)", 
     ylab = "Sample size per group (n)",
     main = "Effect of Correlation on Sample Size",
     ylim = c(90, max(ss_by_rho) * 1.1))
grid()

## ----example3-----------------------------------------------------------------
# Fixed design parameters
delta <- 0.5
p1 <- 0.7
p2 <- 0.5
rho <- 0.5

test_methods <- c("AN", "ANc", "AS", "ASc")

test_comparison <- lapply(test_methods, function(test_method) {
  result <- ss2MixedContinuousBinary(
    delta = delta,
    sd = 1,
    p1 = p1,
    p2 = p2,
    rho = rho,
    r = 1,
    alpha = 0.025,
    beta = 0.2,
    Test = test_method
  )
  
  data.frame(
    Test_method = test_method,
    n_per_group = result$n2,
    N_total = result$N
  )
})

test_comparison_table <- bind_rows(test_comparison)

kable(test_comparison_table,
      caption = "Comparison of Test Methods for Binary Endpoint",
      digits = 0,
      col.names = c("Test Method", "n per group", "N total"))

## ----example4-----------------------------------------------------------------
# Balanced design (r = 1, equivalent to kappa = 1)
result_balanced <- ss2MixedContinuousBinary(
  delta = 0.5,
  sd = 1,
  p1 = 0.7,
  p2 = 0.5,
  rho = 0.5,
  r = 1,
  alpha = 0.025,
  beta = 0.2,
  Test = "AN"
)

# Unbalanced design (r = 2, i.e., nT = 2*nC, equivalent to kappa = 0.5)
result_unbalanced <- ss2MixedContinuousBinary(
  delta = 0.5,
  sd = 1,
  p1 = 0.7,
  p2 = 0.5,
  rho = 0.5,
  r = 2,
  alpha = 0.025,
  beta = 0.2,
  Test = "AN"
)

comparison_allocation <- data.frame(
  Design = c("Balanced (1:1)", "Unbalanced (2:1)"),
  n_treatment = c(result_balanced$n1, result_unbalanced$n1),
  n_control = c(result_balanced$n2, result_unbalanced$n2),
  N_total = c(result_balanced$N, result_unbalanced$N),
  kappa = c(1, 0.5)
)

kable(comparison_allocation,
      caption = "Comparison: Balanced vs Unbalanced Allocation",
      col.names = c("Design", "nT", "nC", "N total", "κ"))

cat("\nIncrease in total sample size:", 
    round((result_unbalanced$N - result_balanced$N) / result_balanced$N * 100, 1), "%\n")

## ----power_verification-------------------------------------------------------
# Use result from Example 1
power_result <- power2MixedContinuousBinary(
  n1 = result_balanced$n1,
  n2 = result_balanced$n2,
  delta = 0.5,
  sd = 1,
  p1 = 0.7,
  p2 = 0.5,
  rho = 0.5,
  alpha = 0.025,
  Test = "AN"
)

cat("Target power: 0.80\n")
cat("Achieved power (Continuous endpoint):", round(as.numeric(power_result$power1), 4), "\n")
cat("Achieved power (Binary endpoint):", round(as.numeric(power_result$power2), 4), "\n")
cat("Achieved power (Co-primary):", round(as.numeric(power_result$powerCoprimary), 4), "\n")