Equivalence & Non-Inferiority

# TOST Equivalence Test: Generic vs Brand-name Bioequivalence
# FDA bioequivalence testing using Two One-Sided Tests (TOST)

library(equivalence)  # For TOST
library(TOSTER)       # Alternative TOST package
library(ggplot2)

# Set seed
set.seed(2025)

# === Simulate crossover bioequivalence data ===
# True ratio ≈ 1.05 (5% higher generic, but within margin)
# Log-normal distribution for AUC

n <- 24  # Number of subjects in crossover design

# Each subject gets both treatments
# True geometric mean ratio = 1.05 (Generic 5% higher than Brand)
log_mean_brand <- 4.5    # log(AUC) for brand
log_mean_generic <- log(exp(log_mean_brand) * 1.05)  # 5% higher
sd_within <- 0.25        # Within-subject SD on log scale

# Subject-specific random effects
subject_effect <- rnorm(n, mean=0, sd=0.3)

# Generate paired data (crossover)
data <- data.frame(
  subject = rep(1:n, each=2),
  treatment = rep(c("Brand", "Generic"), times=n),
  log_AUC = c(
    log_mean_brand + subject_effect + rnorm(n, 0, sd_within),      # Brand
    log_mean_generic + subject_effect + rnorm(n, 0, sd_within)     # Generic
  )
)

# Convert to wide format for paired test
library(tidyr)
data_wide <- data %>%
  pivot_wider(id_cols=subject, names_from=treatment, values_from=log_AUC)

print("=== Sample Data(first 6 subjects) ===")
print(head(data_wide))

# === STEP 1: Check Assumptions ===
cat("\n=== Assumption Checks ===\n")

# 1. Normality of differences (log scale)
differences <- data_wide$Generic - data_wide$Brand
shapiro_result <- shapiro.test(differences)
cat("Shapiro-Wilk test for normality: W =", round(shapiro_result$statistic, 3), 
    ", p =", round(shapiro_result$p.value, 3), 
    ifelse(shapiro_result$p.value > 0.05, "✓", "✗"), "\n")

# Q-Q plot
qqnorm(differences, main="Q-Q Plot: Log(AUC) Differences")
qqline(differences)

cat("\nCrossover design: Each subject as own control ✓\n")
cat("Adequate washout period assumed ✓\n")

# === STEP 2: Define Equivalence Margin ===
# FDA bioequivalence: 80-125% ratio on original scale
# On log scale: log(0.80) to log(1.25)
lower_margin <- log(0.80)  # -0.223
upper_margin <- log(1.25)  # +0.223

cat("\n=== Equivalence Margins(FDA Standard) ===\n")
cat("Original scale: 80% to 125% ratio\n")
cat("Log scale: ", round(lower_margin, 3), " to ", round(upper_margin, 3), "\n")

# === STEP 3: Calculate Mean Difference and 90% CI ===
cat("\n=== Mean Difference(Generic - Brand) ===\n")

mean_diff <- mean(differences)
sd_diff <- sd(differences)
se_diff <- sd_diff / sqrt(n)

# 90% CI (for α = .05 TOST)
t_crit <- qt(0.95, df=n-1)  # One-sided 5%, two-sided 10%
ci_90_lower <- mean_diff - t_crit * se_diff
ci_90_upper <- mean_diff + t_crit * se_diff

cat("Mean log(AUC) difference:", round(mean_diff, 4), "\n")
cat("90% CI:", round(ci_90_lower, 4), "to", round(ci_90_upper, 4), "\n")

# Back-transform to ratio scale
ratio <- exp(mean_diff)
ratio_ci_lower <- exp(ci_90_lower)
ratio_ci_upper <- exp(ci_90_upper)

cat("\nGeometric mean ratio(Generic/Brand):", round(ratio, 3), "\n")
cat("90% CI for ratio:", round(ratio_ci_lower * 100, 1), "% to", 
    round(ratio_ci_upper * 100, 1), "%\n")

# === STEP 4: TOST (Two One-Sided Tests) ===
cat("\n=== TOST Equivalence Test ===\n")

# Test 1: H₀: μ_diff ≤ lower_margin vs Hₐ: μ_diff > lower_margin
t1 <- (mean_diff - lower_margin) / se_diff
p1 <- pt(t1, df=n-1, lower.tail=FALSE)

cat("Test 1 (lower bound): t =", round(t1, 3), ", p =", round(p1, 4), "\n")

# Test 2: H₀: μ_diff ≥ upper_margin vs Hₐ: μ_diff < upper_margin
t2 <- (mean_diff - upper_margin) / se_diff
p2 <- pt(t2, df=n-1, lower.tail=TRUE)

cat("Test 2 (upper bound): t =", round(t2, 3), ", p =", round(p2, 4), "\n")

# TOST p-value = max(p1, p2)
p_tost <- max(p1, p2)
cat("\nTOST p-value(max of two tests):", round(p_tost, 4), "\n")

# === STEP 5: Equivalence Decision ===
cat("\n=== Equivalence Decision ===\n")

if (ci_90_lower > lower_margin & ci_90_upper < upper_margin) {
  cat("CONCLUSION: BIOEQUIVALENT ✓\n")
  cat("90% CI [", round(ci_90_lower, 4), ",", round(ci_90_upper, 4), 
      "] entirely within margins [", round(lower_margin, 3), ",", 
      round(upper_margin, 3), "]\n")
  cat("On ratio scale: 90% CI [", round(ratio_ci_lower * 100, 1), "%, ", 
      round(ratio_ci_upper * 100, 1), "%] within [80%, 125%]\n")
  cat("TOST p-value =", round(p_tost, 4), "< .05: Reject H₀ of non-equivalence\n")
} else {
  cat("CONCLUSION: NOT BIOEQUIVALENT ✗\n")
  cat("90% CI crosses equivalence margin\n")
  cat("Cannot conclude bioequivalence at α = .05 level\n")
}

# === STEP 6: Visualizations ===

# 6.1 Forest plot with equivalence margins
forest_data <- data.frame(
  Study = "Generic vs Brand",
  Estimate = mean_diff,
  CI_lower = ci_90_lower,
  CI_upper = ci_90_upper
)

ggplot(forest_data, aes(y=Study, x=Estimate)) +
  geom_rect(aes(xmin=lower_margin, xmax=upper_margin, ymin=-Inf, ymax=Inf),
            fill="lightgreen", alpha=0.3) +
  geom_point(size=5, color="darkblue") +
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), height=0.2, linewidth=1.5) +
  geom_vline(xintercept=0, linetype="solid", color="black", linewidth=1) +
  geom_vline(xintercept=lower_margin, linetype="dashed", color="red", linewidth=1) +
  geom_vline(xintercept=upper_margin, linetype="dashed", color="red", linewidth=1) +
  annotate("text", x=lower_margin, y=1.3, label="Lower Margin", 
           size=3, hjust=1.1) +
  annotate("text", x=upper_margin, y=1.3, label="Upper Margin", 
           size=3, hjust=-0.1) +
  labs(title="Bioequivalence: Generic vs Brand-name Drug",
       subtitle="90% CI for log(AUC) difference with FDA equivalence margins",
       x="Log(Generic) - Log(Brand)\n[Equivalence Zone in Green]",
       y="") +
  theme_classic(base_size=14) +
  theme(plot.title=element_text(hjust=0.5, face="bold"),
        plot.subtitle=element_text(hjust=0.5))

ggsave("bioequivalence_forest.png", width=12, height=6)

# 6.2 Ratio scale visualization
ggplot(data.frame(x=1), aes(x=x)) +
  geom_rect(aes(xmin=0.80, xmax=1.25, ymin=-Inf, ymax=Inf),
            fill="lightgreen", alpha=0.3) +
  geom_point(aes(x=ratio, y=1), size=5, color="darkblue") +
  geom_errorbarh(aes(y=1, xmin=ratio_ci_lower, xmax=ratio_ci_upper), 
                 height=0.3, linewidth=1.5) +
  geom_vline(xintercept=1, linetype="solid", color="black", linewidth=1) +
  geom_vline(xintercept=0.80, linetype="dashed", color="red", linewidth=1) +
  geom_vline(xintercept=1.25, linetype="dashed", color="red", linewidth=1) +
  scale_x_continuous(limits=c(0.7, 1.4), breaks=seq(0.7, 1.4, 0.1)) +
  labs(title="Bioequivalence: Geometric Mean Ratio",
       subtitle="90% CI must be entirely within 80-125% for bioequivalence",
       x="Ratio: Generic AUC / Brand AUC(%)",
       y="") +
  theme_classic(base_size=14) +
  theme(plot.title=element_text(hjust=0.5, face="bold"),
        plot.subtitle=element_text(hjust=0.5),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

ggsave("bioequivalence_ratio.png", width=12, height=6)

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat(sprintf(
  "A crossover bioequivalence study(n=%d healthy volunteers) compared a generic\nformulation to the brand-name drug using AUC(area under concentration-time curve).\n\nUsing FDA's Two One-Sided Tests(TOST) procedure with equivalence margins of\n80%%-125%% (±20%% on ratio scale), the geometric mean ratio was %.2f (90%% CI [%.2f%%, %.2f%%]).\n\nThe 90%% confidence interval fell entirely within the pre-specified equivalence\nmargins of 80%%-125%% (TOST p = %.4f < .05), supporting bioequivalence of the\ngeneric to the brand-name formulation.\n\nOn the log scale, the mean difference was %.4f (90%% CI [%.4f, %.4f]), within\nthe margins of [%.3f, %.3f].\n\nConclusion: The generic formulation is bioequivalent to the brand-name drug,\nmeeting FDA criteria for therapeutic equivalence.",
  n, ratio, ratio_ci_lower * 100, ratio_ci_upper * 100, p_tost,
  mean_diff, ci_90_lower, ci_90_upper, lower_margin, upper_margin
))