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

## ----correlation_bounds-------------------------------------------------------
# Scenario: lambda = 1.25, nu = 0.8, mu = 0, sigma = 250
bounds1 <- corrbound2MixedCountContinuous(lambda = 1.25, nu = 0.8, mu = 0, sd = 250)
cat("Correlation bounds for NB(1.25, 0.8) and N(0, 250²):\n")
cat("Lower bound:", round(bounds1[1], 3), "\n")
cat("Upper bound:", round(bounds1[2], 3), "\n\n")

# Higher dispersion (smaller nu) typically restricts bounds
bounds2 <- corrbound2MixedCountContinuous(lambda = 1.25, nu = 0.5, mu = 0, sd = 250)
cat("Correlation bounds for NB(1.25, 0.5) and N(0, 250²):\n")
cat("Lower bound:", round(bounds2[1], 3), "\n")
cat("Upper bound:", round(bounds2[2], 3), "\n\n")

# Higher mean count
bounds3 <- corrbound2MixedCountContinuous(lambda = 2.0, nu = 0.8, mu = 0, sd = 250)
cat("Correlation bounds for NB(2.0, 0.8) and N(0, 250²):\n")
cat("Lower bound:", round(bounds3[1], 3), "\n")
cat("Upper bound:", round(bounds3[2], 3), "\n")

## ----table1_case_B------------------------------------------------------------
# Define scenarios for Table 1 Case B
scenarios_table1_B <- expand.grid(
  nu = c(3, 5),
  rho = c(0, 0.2, 0.4, 0.6, 0.8),
  stringsAsFactors = FALSE
)

# Calculate sample sizes for each scenario
results_table1_B <- lapply(1:nrow(scenarios_table1_B), function(i) {
  nu_val <- scenarios_table1_B$nu[i]
  rho_val <- scenarios_table1_B$rho[i]
  
  result <- ss2MixedCountContinuous(
    r1 = 1,
    r2 = 2,
    nu = nu_val,
    t = 1,
    mu1 = -50,
    mu2 = 0,
    sd = 75,
    r = 1,
    rho1 = rho_val,
    rho2 = rho_val,
    alpha = 0.025,
    beta = 0.1
  )
  
  data.frame(
    nu = nu_val,
    rho = rho_val,
    n2 = result$n2,
    n1 = result$n1,
    N = result$N
  )
})

table1_B_results <- bind_rows(results_table1_B)

# Display results grouped by nu
table1_B_nu3 <- table1_B_results %>%
  filter(nu == 3) %>%
  select(rho, n2, N)

table1_B_nu5 <- table1_B_results %>%
  filter(nu == 5) %>%
  select(rho, n2, N)

kable(table1_B_nu3, 
      caption = "Table 1 Case B: Sample Sizes (ν = 3, Balanced Design, α = 0.025, 1-β = 0.9)",
      digits = 1,
      col.names = c("ρ", "n per group", "N total"))

kable(table1_B_nu5, 
      caption = "Table 1 Case B: Sample Sizes (ν = 5, Balanced Design, α = 0.025, 1-β = 0.9)",
      digits = 1,
      col.names = c("ρ", "n per group", "N total"))

## ----example1-----------------------------------------------------------------
# Balanced design: n1 = n2 (i.e., r = 1)
result_balanced <- ss2MixedCountContinuous(
  r1 = 1.0,              # Count rate in treatment group
  r2 = 1.25,             # Count rate in control group
  nu = 0.8,              # Dispersion parameter
  t = 1,                 # Follow-up time
  mu1 = -50,             # Mean for treatment (negative = benefit)
  mu2 = 0,               # Mean for control
  sd = 250,              # Standard deviation
  r = 1,                 # Balanced allocation
  rho1 = 0.5,            # Correlation in treatment group
  rho2 = 0.5,            # Correlation in control group
  alpha = 0.025,
  beta = 0.2
)

print(result_balanced)

## ----example2-----------------------------------------------------------------
# Fixed effect sizes
r1 <- 1.0
r2 <- 1.25
nu <- 0.8
mu1 <- -50
mu2 <- 0
sd <- 250

# Range of correlations
rho_values <- c(0, 0.2, 0.4, 0.6, 0.8)

ss_by_rho <- sapply(rho_values, function(rho) {
  result <- ss2MixedCountContinuous(
    r1 = r1, r2 = r2, nu = nu, t = 1,
    mu1 = mu1, mu2 = mu2, sd = sd,
    r = 1,
    rho1 = rho, rho2 = rho,
    alpha = 0.025,
    beta = 0.2
  )
  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 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 = "Correlation (ρ)", 
     ylab = "Sample size per group (n)",
     main = "Effect of Correlation on Sample Size",
     ylim = c(600, max(ss_by_rho) * 1.1))
grid()

## ----example3-----------------------------------------------------------------
# Fixed design parameters
r1 <- 1.0
r2 <- 1.25
mu1 <- -50
mu2 <- 0
sd <- 250
rho <- 0.5

# Range of dispersion parameters
nu_values <- c(0.5, 0.8, 1.0, 2.0, 5.0)

ss_by_nu <- sapply(nu_values, function(nu) {
  result <- ss2MixedCountContinuous(
    r1 = r1, r2 = r2, nu = nu, t = 1,
    mu1 = mu1, mu2 = mu2, sd = sd,
    r = 1,
    rho1 = rho, rho2 = rho,
    alpha = 0.025,
    beta = 0.2
  )
  result$n2
})

result_df_nu <- data.frame(
  nu = nu_values,
  VMR = round(1 + 1.125/nu_values, 2),  # VMR at lambda = 1.125 (average)
  n_per_group = ss_by_nu,
  N_total = ss_by_nu * 2
)

kable(result_df_nu,
      caption = "Effect of Dispersion Parameter on Sample Size",
      digits = c(1, 2, 0, 0),
      col.names = c("ν", "VMR", "n per group", "N total"))

## ----example4-----------------------------------------------------------------
# Balanced design (r = 1)
result_balanced <- ss2MixedCountContinuous(
  r1 = 1.0, r2 = 1.25, nu = 0.8, t = 1,
  mu1 = -50, mu2 = 0, sd = 250,
  r = 1,
  rho1 = 0.5, rho2 = 0.5,
  alpha = 0.025,
  beta = 0.2
)

# Unbalanced design (r = 2, i.e., n1 = 2*n2)
result_unbalanced <- ss2MixedCountContinuous(
  r1 = 1.0, r2 = 1.25, nu = 0.8, t = 1,
  mu1 = -50, mu2 = 0, sd = 250,
  r = 2,
  rho1 = 0.5, rho2 = 0.5,
  alpha = 0.025,
  beta = 0.2
)

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

kable(comparison_allocation,
      caption = "Comparison: Balanced vs Unbalanced Allocation",
      col.names = c("Design", "n₁", "n₂", "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 <- power2MixedCountContinuous(
  n1 = result_balanced$n1,
  n2 = result_balanced$n2,
  r1 = 1.0,
  r2 = 1.25,
  nu = 0.8,
  t = 1,
  mu1 = -50,
  mu2 = 0,
  sd = 250,
  rho1 = 0.5,
  rho2 = 0.5,
  alpha = 0.025
)

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