BE Simulation Report
2025-10-27
Summary of Simulation Workflow
This R script implements a reproducible framework for virtual bioequivalence (BE) assessment using physiologically based pharmacokinetic (PBPK)–inspired concentration–time profiles and Monte Carlo methods.
Step 1 – Define Reference PK Profile
A mean reference concentration–time profile is provided (e.g.,
ref <- data.frame(time_hr, mean_conc_ng_per_mL)). This
serves as the baseline for simulating both Test and
Reference arms under user-specified variability
assumptions.
Step 2 – Set Simulation Parameters
Key parameters include: - n_per_arm – sample size per
treatment arm
- CV_between, CV_resid – between-subject and
residual variability (coefficient of variation)
- true_GMR – assumed geometric mean ratio
(Test/Reference)
- BE acceptance limits (default 0.80 – 1.25) and type I error
alpha
Step 3 – Monte Carlo Trial Generation
For each virtual trial: 1. Individual concentration profiles are
generated with log-normal random effects (CV_between) and
residual error (CV_resid). 2. Individual AUC and Cmax
values are calculated using trapezoidal integration
(trap_auc()). 3. Bioequivalence statistics (90 % CI of
log-transformed ratios) are computed using a standard parallel-design
t-test (be_stats_parallel()). 4. Pass/fail outcomes for
AUC, Cmax, and joint BE are recorded.
Step 4 – Decision Probability Estimation
Multiple (e.g., > 1000) virtual trials are simulated to estimate the empirical probability of meeting BE criteria for the specified design and parameters: \[ P_\text{BE} = \Pr(90\%\,CI \subseteq [0.80,\,1.25]). \]
Step 5 – Decision Heatmap
A two-dimensional grid of true GMR values and sample sizes is explored. For each grid cell, the Monte Carlo procedure is repeated, yielding a heatmap of \(P_\text{BE}\) for AUC and Cmax. This allows visualization of power and design sensitivity across plausible conditions.
Step 6 – Parameter-Uncertainty Propagation
To evaluate robustness, uncertainty in variability parameters is
propagated by repeatedly sampling plausible CV_between and
CV_resid values (log-normal distribution around
nominal).
Each sampled pair triggers a short Monte Carlo run, generating a
distribution of \(P_\text{BE}\)
estimates.
Step 7 – Visualization and Reporting
Plots include: - Histograms of simulated Cmax and AUC ratios with BE
limits
- Heatmaps of \(P_\text{BE}\) vs GMR
and sample size
- Sensitivity plots showing impact of variability on BE probability
BE Simulation Report
# virtual BE Monte Carlo with:
# ------------------------------------------------------------------------------
# ---------------------------
# PACKAGES
# ---------------------------
library(ggplot2)
# ---------------------------
# 0. EXAMPLE REFERENCE PROFILE
# ---------------------------
ref <- data.frame(
time_hr = c(0, 0.25, 0.5, 1, 2, 4, 6, 8, 12, 24),
mean_conc_ng_per_mL = c(0, 12.5, 28.3, 42.8, 36.0, 20.1, 12.6, 7.8, 2.1, 0.2)
)
# ---------------------------
# 1. USER PARAMETERS (tweak for speed/precision)
# ---------------------------
set.seed(12345)
# Basic simulation parameters
N_trials_main <- 1200 # Monte Carlo trials for main simulation (per baseline scenario)
n_per_arm <- 24
CV_between_nom <- 0.25
CV_resid_nom <- 0.20
true_GMR_nom <- 0.95
alpha <- 0.10
BE_lower <- 0.80
BE_upper <- 1.25
# Extended analyses parameters
# 1) Heatmap grid for decision visualization
gmr_grid <- seq(0.88, 1.04, length.out = 9) # range of true GMRs to explore
n_grid <- seq(12, 60, by = 6) # sample sizes per arm to evaluate
reps_per_cell <- 250 # Monte Carlo replicates per (GMR, n) cell
# 2) Parameter uncertainty propagation
param_sims <- 60 # number of different (CV_between, CV_resid) draws
trials_per_param_sim <- 600 # MC trials within each param simulation
# Files
# results_dir <- "."
# results_csv_main <- file.path(results_dir, "results_summary_virtual_BE_main.csv")
# results_csv_heatmap <- file.path(results_dir, "results_heatmap_BE.csv")
# results_csv_param <- file.path(results_dir, "results_param_uncertainty.csv")
# report_rmd <- file.path(results_dir, "virtual_BE_report.Rmd")
# report_html <- file.path(results_dir, "virtual_BE_report.html")
# ---------------------------
# 2. HELPER FUNCTIONS
# ---------------------------
trap_auc <- function(t, conc){
sum(0.5 * (conc[-1] + conc[-length(conc)]) * diff(t))
}
cv_to_lnsd <- function(cv) sqrt(log(1 + cv^2))
be_stats_parallel <- function(yT, yR, alpha = 0.10){
# returns ratio, lower, upper (on original scale)
lnT <- log(yT); lnR <- log(yR)
nT <- length(lnT); nR <- length(lnR)
meanT <- mean(lnT); meanR <- mean(lnR)
s2T <- var(lnT); s2R <- var(lnR)
df <- nT + nR - 2
if (df <= 0) return(list(ratio = NA, lower = NA, upper = NA))
s_pooled <- sqrt(((nT - 1) * s2T + (nR - 1) * s2R) / df)
se_diff <- s_pooled * sqrt(1/nT + 1/nR)
t_crit <- qt(1 - alpha/2, df)
diff <- meanT - meanR
lower_ln <- diff - t_crit * se_diff
upper_ln <- diff + t_crit * se_diff
list(ratio = exp(diff), lower = exp(lower_ln), upper = exp(upper_ln))
}
simulate_trial_parallel <- function(mean_profile, n_per_arm, CV_between, CV_resid, true_GMR, alpha){
# Simulate one virtual two-arm (parallel) trial and return BE stats for AUC and Cmax
times <- mean_profile$time_hr
mean_conc <- mean_profile$mean_conc_ng_per_mL
sd_between_ln <- cv_to_lnsd(CV_between)
sd_resid_ln <- cv_to_lnsd(CV_resid)
# subject-level scaling and residuals
etaR <- exp(rnorm(n_per_arm, 0, sd_between_ln))
etaT <- exp(rnorm(n_per_arm, 0, sd_between_ln))
aucR <- aucT <- cmaxR <- cmaxT <- numeric(n_per_arm)
for (i in seq_len(n_per_arm)) {
residR <- exp(rnorm(length(mean_conc), 0, sd_resid_ln))
residT <- exp(rnorm(length(mean_conc), 0, sd_resid_ln))
concR <- mean_conc * etaR[i] * residR
concT <- mean_conc * true_GMR * etaT[i] * residT
aucR[i] <- trap_auc(times, concR); aucT[i] <- trap_auc(times, concT)
cmaxR[i] <- max(concR); cmaxT[i] <- max(concT)
}
bs_auc <- be_stats_parallel(aucT, aucR, alpha)
bs_cmax <- be_stats_parallel(cmaxT, cmaxR, alpha)
pass_auc <- (bs_auc$lower >= BE_lower) & (bs_auc$upper <= BE_upper)
pass_cmax <- (bs_cmax$lower >= BE_lower) & (bs_cmax$upper <= BE_upper)
list(
ratio_auc = bs_auc$ratio, lower_auc = bs_auc$lower, upper_auc = bs_auc$upper,
ratio_cmax = bs_cmax$ratio, lower_cmax = bs_cmax$lower, upper_cmax = bs_cmax$upper,
pass_auc = pass_auc, pass_cmax = pass_cmax, pass_both = (pass_auc & pass_cmax)
)
}
# ---------------------------
# 3. MAIN / baseline simulation (single scenario)
# ---------------------------
cat("Running main baseline Monte Carlo simulation...\n")## Running main baseline Monte Carlo simulation...main_results <- data.frame(
trial = integer(), ratio_auc = numeric(), lower_auc = numeric(), upper_auc = numeric(),
ratio_cmax = numeric(), lower_cmax = numeric(), upper_cmax = numeric(),
pass_auc = logical(), pass_cmax = logical(), pass_both = logical(),
stringsAsFactors = FALSE
)
for (tr in seq_len(N_trials_main)){
res <- simulate_trial_parallel(ref, n_per_arm, CV_between_nom, CV_resid_nom, true_GMR_nom, alpha)
main_results[tr, ] <- c(tr, res$ratio_auc, res$lower_auc, res$upper_auc,
res$ratio_cmax, res$lower_cmax, res$upper_cmax,
res$pass_auc, res$pass_cmax, res$pass_both)
}
main_results$pass_auc <- as.logical(main_results$pass_auc)
main_results$pass_cmax <- as.logical(main_results$pass_cmax)
main_results$pass_both <- as.logical(main_results$pass_both)
print(head(main_results,100 ))## trial ratio_auc lower_auc upper_auc ratio_cmax lower_cmax upper_cmax
## 1 1 0.9996400 0.8721453 1.1457726 0.9484519 0.8221302 1.0941830
## 2 2 0.9423437 0.8230330 1.0789501 0.8455979 0.7300059 0.9794932
## 3 3 0.9810255 0.8558064 1.1245663 0.9902171 0.8436694 1.1622206
## 4 4 1.0048606 0.8982793 1.1240879 1.0270765 0.9248383 1.1406170
## 5 5 0.8421703 0.7467653 0.9497640 0.8353708 0.7240421 0.9638174
## 6 6 0.9510506 0.8475723 1.0671623 1.0019045 0.8693811 1.1546290
## 7 7 1.0765091 0.9490575 1.2210766 1.0488692 0.9055793 1.2148319
## 8 8 0.8930954 0.7918531 1.0072819 0.9103803 0.7885956 1.0509724
## 9 9 0.9080836 0.7922595 1.0408405 0.8695043 0.7580529 0.9973417
## 10 10 1.0223112 0.8969188 1.1652339 0.9938903 0.8470271 1.1662176
## 11 11 1.0136316 0.9045262 1.1358975 1.0124893 0.8758517 1.1704432
## 12 12 0.9219964 0.8042525 1.0569782 0.9642832 0.8287602 1.1219676
## 13 13 0.9695222 0.8389038 1.1204780 1.0152967 0.8713500 1.1830233
## 14 14 0.9911238 0.8685074 1.1310513 0.9675010 0.8321372 1.1248844
## 15 15 1.1001309 0.9822153 1.2322024 1.1206050 0.9828989 1.2776040
## 16 16 1.0106199 0.8678090 1.1769323 0.9991199 0.8493050 1.1753618
## 17 17 1.0562776 0.9255271 1.2054995 1.1197668 0.9629822 1.3020777
## 18 18 0.9342779 0.8083973 1.0797602 0.8929516 0.7503722 1.0626226
## 19 19 0.9715708 0.8579193 1.1002782 0.9644112 0.8340950 1.1150876
## 20 20 0.9144557 0.8319846 1.0051018 0.8798880 0.7762186 0.9974031
## 21 21 0.9413700 0.8335212 1.0631734 0.9303071 0.8103606 1.0680077
## 22 22 0.9795409 0.8672333 1.1063923 1.0142430 0.8772442 1.1726370
## 23 23 1.0546470 0.9504485 1.1702689 1.0836880 0.9673398 1.2140301
## 24 24 0.8984462 0.7979952 1.0115419 0.8768798 0.7676982 1.0015891
## 25 25 0.9123308 0.8057517 1.0330076 0.8926120 0.7638807 1.0430374
## 26 26 0.9485215 0.8350721 1.0773837 0.9167786 0.7859203 1.0694253
## 27 27 0.8679907 0.7606143 0.9905254 0.9057961 0.7809138 1.0506493
## 28 28 0.8948654 0.7866855 1.0179214 0.9108773 0.7924102 1.0470554
## 29 29 0.9460754 0.8500163 1.0529900 0.9357989 0.8232247 1.0637673
## 30 30 0.9560771 0.8447407 1.0820875 0.9270054 0.8021606 1.0712806
## 31 31 0.9676872 0.8600429 1.0888045 1.0388654 0.8951195 1.2056951
## 32 32 1.0069046 0.8939354 1.1341500 0.9805534 0.8571904 1.1216703
## 33 33 1.0059672 0.8789686 1.1513152 0.9876074 0.8425045 1.1577010
## 34 34 1.0069638 0.8785046 1.1542070 0.9520043 0.8221502 1.1023683
## 35 35 1.0129895 0.8927303 1.1494489 1.0395784 0.8941883 1.2086081
## 36 36 0.9785023 0.8591207 1.1144729 0.9995708 0.8547504 1.1689282
## 37 37 0.8216309 0.7462151 0.9046687 0.8129649 0.7261257 0.9101894
## 38 38 1.0189929 0.9115099 1.1391501 1.0097123 0.8985741 1.1345963
## 39 39 0.8167680 0.7388491 0.9029043 0.8176754 0.7241564 0.9232716
## 40 40 0.9728685 0.8661685 1.0927124 1.0385144 0.8920392 1.2090413
## 41 41 0.9765901 0.8460085 1.1273270 1.0326371 0.8669368 1.2300083
## 42 42 1.0180999 0.8962642 1.1564976 1.0712887 0.9165557 1.2521438
## 43 43 0.8798687 0.7922633 0.9771611 0.8223522 0.7332199 0.9223196
## 44 44 0.8882160 0.7961381 0.9909433 0.9133030 0.8015471 1.0406405
## 45 45 0.8745648 0.7761047 0.9855160 0.8705844 0.7722481 0.9814425
## 46 46 0.9150932 0.8152171 1.0272055 0.9595930 0.8486329 1.0850612
## 47 47 0.8955942 0.7799648 1.0283656 0.8404682 0.7163910 0.9860353
## 48 48 1.1716428 1.0352536 1.3260005 1.1704581 1.0136686 1.3514989
## 49 49 0.8741434 0.7920879 0.9646994 0.8389449 0.7501178 0.9382908
## 50 50 0.9242248 0.8106663 1.0536905 0.9518864 0.8124867 1.1152031
## 51 51 0.9201792 0.8172200 1.0361100 0.9353184 0.8116325 1.0778529
## 52 52 0.8991386 0.7933062 1.0190897 0.9242726 0.8106438 1.0538288
## 53 53 0.8249262 0.7236626 0.9403599 0.8226404 0.7068919 0.9573418
## 54 54 0.9251058 0.8366737 1.0228848 0.8902833 0.8018308 0.9884933
## 55 55 1.0024678 0.8933662 1.1248932 0.9909191 0.8801020 1.1156895
## 56 56 0.9764114 0.8594218 1.1093263 0.9289047 0.8145358 1.0593320
## 57 57 0.9683782 0.8458149 1.1087016 0.9299558 0.8037335 1.0760007
## 58 58 0.9599863 0.8317893 1.1079413 0.9348435 0.8006681 1.0915040
## 59 59 1.0185551 0.9069054 1.1439501 1.0438823 0.9037069 1.2058004
## 60 60 1.0168050 0.8944918 1.1558435 1.0392360 0.8964410 1.2047769
## 61 61 0.9253665 0.8176214 1.0473102 0.8964414 0.7720881 1.0408233
## 62 62 0.9332835 0.8411144 1.0355523 0.8880971 0.7817532 1.0089072
## 63 63 0.9032999 0.8200066 0.9950538 0.8432147 0.7440659 0.9555753
## 64 64 1.0631530 0.9209598 1.2273002 1.0469863 0.9012850 1.2162417
## 65 65 1.0080553 0.8875280 1.1449504 1.0367619 0.9100888 1.1810663
## 66 66 0.9148773 0.8101154 1.0331866 0.9158615 0.8058761 1.0408575
## 67 67 0.9225582 0.8135733 1.0461424 0.9380631 0.7953725 1.1063524
## 68 68 0.9990127 0.8752863 1.1402285 1.0818986 0.9394832 1.2459025
## 69 69 0.9750641 0.8627198 1.1020380 0.9114141 0.7955294 1.0441797
## 70 70 0.8787362 0.7686590 1.0045774 0.8618121 0.7493083 0.9912077
## 71 71 0.8266504 0.7416535 0.9213883 0.7913593 0.6921477 0.9047918
## 72 72 1.0715989 0.9428104 1.2179800 1.0993930 0.9399805 1.2858405
## 73 73 0.9799456 0.8757144 1.0965829 0.9369447 0.8220391 1.0679120
## 74 74 0.8851030 0.7780509 1.0068843 0.9002512 0.7744694 1.0464612
## 75 75 0.8413288 0.7468684 0.9477360 0.7836862 0.6861794 0.8950489
## 76 76 0.9492253 0.8398251 1.0728765 0.9622272 0.8330251 1.1114685
## 77 77 0.8662413 0.7489408 1.0019135 0.8854094 0.7458647 1.0510617
## 78 78 1.0106588 0.8796514 1.1611773 1.0215997 0.8845988 1.1798184
## 79 79 1.0420038 0.9328404 1.1639419 1.0865821 0.9558506 1.2351937
## 80 80 1.0760951 0.9468704 1.2229559 1.1324566 0.9774782 1.3120067
## 81 81 0.9816567 0.8795259 1.0956469 1.0450151 0.9134710 1.1955022
## 82 82 0.8386049 0.7364805 0.9548904 0.8578429 0.7553655 0.9742231
## 83 83 0.8985804 0.7905647 1.0213544 0.9135978 0.8084490 1.0324226
## 84 84 0.9484188 0.8436235 1.0662318 0.9924687 0.8562457 1.1503638
## 85 85 0.8635584 0.7543963 0.9885164 0.8752859 0.7573886 1.0115355
## 86 86 1.0047517 0.8993557 1.1224992 1.0448644 0.9146610 1.1936024
## 87 87 0.9208991 0.7986010 1.0619259 0.8696044 0.7368383 1.0262927
## 88 88 0.9173720 0.8002384 1.0516509 1.0125480 0.8673811 1.1820102
## 89 89 0.8889553 0.7968673 0.9916853 0.9025403 0.7842374 1.0386893
## 90 90 0.8657428 0.7601392 0.9860177 0.8692651 0.7468528 1.0117413
## 91 91 0.9980374 0.8688702 1.1464068 0.9481586 0.8143984 1.1038882
## 92 92 1.0269357 0.9196189 1.1467759 0.9913264 0.8740126 1.1243865
## 93 93 1.0273328 0.9070115 1.1636157 0.9861100 0.8636622 1.1259182
## 94 94 1.0164660 0.8871071 1.1646882 1.0216166 0.8692655 1.2006693
## 95 95 1.0435369 0.9092126 1.1977058 1.0988241 0.9493142 1.2718807
## 96 96 0.8858391 0.8050860 0.9746920 0.9335381 0.8331965 1.0459637
## 97 97 0.9519321 0.8467025 1.0702399 0.9431117 0.8233805 1.0802535
## 98 98 0.8972216 0.7875678 1.0221425 0.8970517 0.7573655 1.0625011
## 99 99 0.8613358 0.7475197 0.9924814 0.8699290 0.7448589 1.0159997
## 100 100 0.9916079 0.8667607 1.1344379 0.9948726 0.8558135 1.1565270
## pass_auc pass_cmax pass_both
## 1 TRUE TRUE TRUE
## 2 TRUE FALSE FALSE
## 3 TRUE TRUE TRUE
## 4 TRUE TRUE TRUE
## 5 FALSE FALSE FALSE
## 6 TRUE TRUE TRUE
## 7 TRUE TRUE TRUE
## 8 FALSE FALSE FALSE
## 9 FALSE FALSE FALSE
## 10 TRUE TRUE TRUE
## 11 TRUE TRUE TRUE
## 12 TRUE TRUE TRUE
## 13 TRUE TRUE TRUE
## 14 TRUE TRUE TRUE
## 15 TRUE FALSE FALSE
## 16 TRUE TRUE TRUE
## 17 TRUE FALSE FALSE
## 18 TRUE FALSE FALSE
## 19 TRUE TRUE TRUE
## 20 TRUE FALSE FALSE
## 21 TRUE TRUE TRUE
## 22 TRUE TRUE TRUE
## 23 TRUE TRUE TRUE
## 24 FALSE FALSE FALSE
## 25 TRUE FALSE FALSE
## 26 TRUE FALSE FALSE
## 27 FALSE FALSE FALSE
## 28 FALSE FALSE FALSE
## 29 TRUE TRUE TRUE
## 30 TRUE TRUE TRUE
## 31 TRUE TRUE TRUE
## 32 TRUE TRUE TRUE
## 33 TRUE TRUE TRUE
## 34 TRUE TRUE TRUE
## 35 TRUE TRUE TRUE
## 36 TRUE TRUE TRUE
## 37 FALSE FALSE FALSE
## 38 TRUE TRUE TRUE
## 39 FALSE FALSE FALSE
## 40 TRUE TRUE TRUE
## 41 TRUE TRUE TRUE
## 42 TRUE FALSE FALSE
## 43 FALSE FALSE FALSE
## 44 FALSE TRUE FALSE
## 45 FALSE FALSE FALSE
## 46 TRUE TRUE TRUE
## 47 FALSE FALSE FALSE
## 48 FALSE FALSE FALSE
## 49 FALSE FALSE FALSE
## 50 TRUE TRUE TRUE
## 51 TRUE TRUE TRUE
## 52 FALSE TRUE FALSE
## 53 FALSE FALSE FALSE
## 54 TRUE TRUE TRUE
## 55 TRUE TRUE TRUE
## 56 TRUE TRUE TRUE
## 57 TRUE TRUE TRUE
## 58 TRUE TRUE TRUE
## 59 TRUE TRUE TRUE
## 60 TRUE TRUE TRUE
## 61 TRUE FALSE FALSE
## 62 TRUE FALSE FALSE
## 63 TRUE FALSE FALSE
## 64 TRUE TRUE TRUE
## 65 TRUE TRUE TRUE
## 66 TRUE TRUE TRUE
## 67 TRUE FALSE FALSE
## 68 TRUE TRUE TRUE
## 69 TRUE FALSE FALSE
## 70 FALSE FALSE FALSE
## 71 FALSE FALSE FALSE
## 72 TRUE FALSE FALSE
## 73 TRUE TRUE TRUE
## 74 FALSE FALSE FALSE
## 75 FALSE FALSE FALSE
## 76 TRUE TRUE TRUE
## 77 FALSE FALSE FALSE
## 78 TRUE TRUE TRUE
## 79 TRUE TRUE TRUE
## 80 TRUE FALSE FALSE
## 81 TRUE TRUE TRUE
## 82 FALSE FALSE FALSE
## 83 FALSE TRUE FALSE
## 84 TRUE TRUE TRUE
## 85 FALSE FALSE FALSE
## 86 TRUE TRUE TRUE
## 87 FALSE FALSE FALSE
## 88 TRUE TRUE TRUE
## 89 FALSE FALSE FALSE
## 90 FALSE FALSE FALSE
## 91 TRUE TRUE TRUE
## 92 TRUE TRUE TRUE
## 93 TRUE TRUE TRUE
## 94 TRUE TRUE TRUE
## 95 TRUE FALSE FALSE
## 96 TRUE TRUE TRUE
## 97 TRUE TRUE TRUE
## 98 FALSE FALSE FALSE
## 99 FALSE FALSE FALSE
## 100 TRUE TRUE TRUE# Summary probabilities
prob_pass_auc_main <- mean(main_results$pass_auc, na.rm = TRUE)
prob_pass_cmax_main <- mean(main_results$pass_cmax, na.rm = TRUE)
prob_pass_both_main <- mean(main_results$pass_both, na.rm = TRUE)
cat(sprintf("Baseline P(pass AUC) = %.3f | P(pass Cmax) = %.3f | P(pass both) = %.3f\n",
prob_pass_auc_main, prob_pass_cmax_main, prob_pass_both_main))## Baseline P(pass AUC) = 0.682 | P(pass Cmax) = 0.557 | P(pass both) = 0.522# ---------------------------
# 4. Heatmap: P(BE) vs true_GMR and sample size (Cmax & AUC)
# ---------------------------
cat("Computing heatmap grid (this may take a while)...\n")## Computing heatmap grid (this may take a while)...heatmap_rows <- expand.grid(true_GMR = gmr_grid, n = n_grid)
heatmap_rows$prob_cmax <- NA
heatmap_rows$prob_auc <- NA
# iterate cells
cell_idx <- 0
total_cells <- nrow(heatmap_rows)
for (r in seq_len(total_cells)){
cell_idx <- cell_idx + 1
gmr_val <- heatmap_rows$true_GMR[r]
n_val <- heatmap_rows$n[r]
pass_cmax_cnt <- 0
pass_auc_cnt <- 0
for (rep in seq_len(reps_per_cell)){
res <- simulate_trial_parallel(ref, n_val, CV_between_nom, CV_resid_nom, gmr_val, alpha)
pass_cmax_cnt <- pass_cmax_cnt + as.integer(res$pass_cmax)
pass_auc_cnt <- pass_auc_cnt + as.integer(res$pass_auc)
}
heatmap_rows$prob_cmax[r] <- pass_cmax_cnt / reps_per_cell
heatmap_rows$prob_auc[r] <- pass_auc_cnt / reps_per_cell
if (cell_idx %% 8 == 0) cat(sprintf("Completed %d / %d cells...\n", cell_idx, total_cells))
}## Completed 8 / 81 cells...
## Completed 16 / 81 cells...
## Completed 24 / 81 cells...
## Completed 32 / 81 cells...
## Completed 40 / 81 cells...
## Completed 48 / 81 cells...
## Completed 56 / 81 cells...
## Completed 64 / 81 cells...
## Completed 72 / 81 cells...
## Completed 80 / 81 cells...print(heatmap_rows )## true_GMR n prob_cmax prob_auc
## 1 0.88 12 0.080 0.124
## 2 0.90 12 0.108 0.176
## 3 0.92 12 0.104 0.212
## 4 0.94 12 0.088 0.276
## 5 0.96 12 0.116 0.260
## 6 0.98 12 0.132 0.352
## 7 1.00 12 0.172 0.336
## 8 1.02 12 0.156 0.320
## 9 1.04 12 0.132 0.328
## 10 0.88 18 0.236 0.300
## 11 0.90 18 0.252 0.344
## 12 0.92 18 0.296 0.420
## 13 0.94 18 0.384 0.504
## 14 0.96 18 0.428 0.608
## 15 0.98 18 0.424 0.608
## 16 1.00 18 0.424 0.636
## 17 1.02 18 0.376 0.632
## 18 1.04 18 0.440 0.596
## 19 0.88 24 0.284 0.320
## 20 0.90 24 0.364 0.472
## 21 0.92 24 0.404 0.488
## 22 0.94 24 0.536 0.696
## 23 0.96 24 0.576 0.768
## 24 0.98 24 0.656 0.832
## 25 1.00 24 0.668 0.804
## 26 1.02 24 0.580 0.744
## 27 1.04 24 0.596 0.780
## 28 0.88 30 0.300 0.384
## 29 0.90 30 0.492 0.584
## 30 0.92 30 0.548 0.668
## 31 0.94 30 0.680 0.776
## 32 0.96 30 0.736 0.860
## 33 0.98 30 0.748 0.868
## 34 1.00 30 0.768 0.920
## 35 1.02 30 0.764 0.892
## 36 1.04 30 0.684 0.816
## 37 0.88 36 0.400 0.472
## 38 0.90 36 0.560 0.656
## 39 0.92 36 0.660 0.728
## 40 0.94 36 0.648 0.788
## 41 0.96 36 0.820 0.900
## 42 0.98 36 0.844 0.928
## 43 1.00 36 0.872 0.952
## 44 1.02 36 0.844 0.952
## 45 1.04 36 0.812 0.920
## 46 0.88 42 0.428 0.508
## 47 0.90 42 0.520 0.616
## 48 0.92 42 0.628 0.744
## 49 0.94 42 0.812 0.908
## 50 0.96 42 0.876 0.932
## 51 0.98 42 0.908 0.980
## 52 1.00 42 0.932 0.996
## 53 1.02 42 0.908 0.956
## 54 1.04 42 0.900 0.960
## 55 0.88 48 0.508 0.580
## 56 0.90 48 0.604 0.712
## 57 0.92 48 0.696 0.808
## 58 0.94 48 0.876 0.948
## 59 0.96 48 0.932 0.964
## 60 0.98 48 0.948 0.980
## 61 1.00 48 0.960 0.988
## 62 1.02 48 0.952 0.988
## 63 1.04 48 0.904 0.972
## 64 0.88 54 0.516 0.604
## 65 0.90 54 0.660 0.760
## 66 0.92 54 0.764 0.880
## 67 0.94 54 0.860 0.908
## 68 0.96 54 0.936 0.968
## 69 0.98 54 0.960 0.996
## 70 1.00 54 0.988 1.000
## 71 1.02 54 0.964 0.992
## 72 1.04 54 0.928 0.984
## 73 0.88 60 0.536 0.632
## 74 0.90 60 0.712 0.824
## 75 0.92 60 0.840 0.928
## 76 0.94 60 0.920 0.972
## 77 0.96 60 0.964 0.976
## 78 0.98 60 0.960 0.988
## 79 1.00 60 0.988 1.000
## 80 1.02 60 0.984 0.992
## 81 1.04 60 0.972 0.996# ---------------------------
# 5. Parameter uncertainty propagation
# ---------------------------
cat("Running parameter uncertainty propagation...\n")## Running parameter uncertainty propagation...# Assume uncertainty on CVs ~ lognormal around nominal with 20% CV on each CV
param_results <- data.frame(sim = integer(), CV_between = numeric(), CV_resid = numeric(),
prob_pass_auc = numeric(), prob_pass_cmax = numeric(), prob_pass_both = numeric(),
stringsAsFactors = FALSE)
for (s in seq_len(param_sims)){
# sample CVs in a plausible range (lognormal multiplicative uncertainty)
CV_between_s <- CV_between_nom * exp(rnorm(1, 0, 0.2)) # ~20% multiplicative uncertainty
CV_resid_s <- CV_resid_nom * exp(rnorm(1, 0, 0.2))
pass_auc_cnt <- pass_cmax_cnt <- pass_both_cnt <- 0
for (tr in seq_len(trials_per_param_sim)){
res <- simulate_trial_parallel(ref, n_per_arm, CV_between_s, CV_resid_s, true_GMR_nom, alpha)
pass_auc_cnt <- pass_auc_cnt + as.integer(res$pass_auc)
pass_cmax_cnt <- pass_cmax_cnt + as.integer(res$pass_cmax)
pass_both_cnt <- pass_both_cnt + as.integer(res$pass_both)
}
param_results[s, ] <- c(s, CV_between_s, CV_resid_s,
pass_auc_cnt / trials_per_param_sim,
pass_cmax_cnt / trials_per_param_sim,
pass_both_cnt / trials_per_param_sim)
if (s %% 10 == 0) cat(sprintf("Completed %d / %d parameter sims...\n", s, param_sims))
}## Completed 10 / 60 parameter sims...
## Completed 20 / 60 parameter sims...
## Completed 30 / 60 parameter sims...
## Completed 40 / 60 parameter sims...
## Completed 50 / 60 parameter sims...
## Completed 60 / 60 parameter sims...names(param_results) <- c("sim", "CV_between", "CV_resid", "prob_pass_auc", "prob_pass_cmax", "prob_pass_both")
print(param_results )## sim CV_between CV_resid prob_pass_auc prob_pass_cmax prob_pass_both
## 1 1 0.3293446 0.1762974 0.4716667 0.34166667 0.29166667
## 2 2 0.2364205 0.2140366 0.7733333 0.59166667 0.56333333
## 3 3 0.2411367 0.1501974 0.7700000 0.66666667 0.63333333
## 4 4 0.2630260 0.2217387 0.6433333 0.49166667 0.44333333
## 5 5 0.2249008 0.2208061 0.8000000 0.60833333 0.59000000
## 6 6 0.2220461 0.2585150 0.7683333 0.55500000 0.53166667
## 7 7 0.2354917 0.2624163 0.7200000 0.50166667 0.47500000
## 8 8 0.1880698 0.2399441 0.8833333 0.70833333 0.69666667
## 9 9 0.2275708 0.1370904 0.8100000 0.72000000 0.69833333
## 10 10 0.2545568 0.2210631 0.7216667 0.53666667 0.51000000
## 11 11 0.1737587 0.2593253 0.8766667 0.64666667 0.62500000
## 12 12 0.2482597 0.1982223 0.7150000 0.57333333 0.52833333
## 13 13 0.3553104 0.1955019 0.3416667 0.23333333 0.18666667
## 14 14 0.3542465 0.2224868 0.3516667 0.24166667 0.19833333
## 15 15 0.2544213 0.1928533 0.7100000 0.56833333 0.54166667
## 16 16 0.4279056 0.1847380 0.1450000 0.07166667 0.04833333
## 17 17 0.2477182 0.1841773 0.7233333 0.58000000 0.55000000
## 18 18 0.2723122 0.2084648 0.6266667 0.47166667 0.42833333
## 19 19 0.2570297 0.3396544 0.6266667 0.30333333 0.27500000
## 20 20 0.2906534 0.1566050 0.6266667 0.53333333 0.49166667
## 21 21 0.2767014 0.1831162 0.6200000 0.46500000 0.43166667
## 22 22 0.2878676 0.1522014 0.6033333 0.53333333 0.50000000
## 23 23 0.2035561 0.2575204 0.8466667 0.63666667 0.61000000
## 24 24 0.3648930 0.1663327 0.3550000 0.28000000 0.21500000
## 25 25 0.2176315 0.1453647 0.8516667 0.73666667 0.71833333
## 26 26 0.2665153 0.2184781 0.6450000 0.48166667 0.45333333
## 27 27 0.2986238 0.2408465 0.5350000 0.36500000 0.31333333
## 28 28 0.2856638 0.1891054 0.5583333 0.44500000 0.39666667
## 29 29 0.2296778 0.2223895 0.7833333 0.56333333 0.53833333
## 30 30 0.2745071 0.2059293 0.6333333 0.51166667 0.47833333
## 31 31 0.2509933 0.1931794 0.6933333 0.59166667 0.54833333
## 32 32 0.2309428 0.1866149 0.7700000 0.67666667 0.64000000
## 33 33 0.2886597 0.2208823 0.5783333 0.44333333 0.39500000
## 34 34 0.1964760 0.2058530 0.8916667 0.73666667 0.71666667
## 35 35 0.3193155 0.1811458 0.4566667 0.37500000 0.32166667
## 36 36 0.2432041 0.2096628 0.7000000 0.54000000 0.51000000
## 37 37 0.1995232 0.2202051 0.8600000 0.67500000 0.66333333
## 38 38 0.1735615 0.1895485 0.9266667 0.78833333 0.78333333
## 39 39 0.2900513 0.1552624 0.5733333 0.48166667 0.43666667
## 40 40 0.2965755 0.2136521 0.5416667 0.42833333 0.37666667
## 41 41 0.2952624 0.2586338 0.5033333 0.30666667 0.27000000
## 42 42 0.2053752 0.2308131 0.8433333 0.63000000 0.61333333
## 43 43 0.3054573 0.1842097 0.5266667 0.46166667 0.36333333
## 44 44 0.2521400 0.1999093 0.6866667 0.52166667 0.49500000
## 45 45 0.2126741 0.2454379 0.7750000 0.56000000 0.52000000
## 46 46 0.2494217 0.3142493 0.6466667 0.38333333 0.34833333
## 47 47 0.2138837 0.1787863 0.8183333 0.69500000 0.67500000
## 48 48 0.2168240 0.2230066 0.7900000 0.61333333 0.59000000
## 49 49 0.1664721 0.2432965 0.9016667 0.67000000 0.65833333
## 50 50 0.2484772 0.1876833 0.7366667 0.58666667 0.55000000
## 51 51 0.2696719 0.1683350 0.6400000 0.53333333 0.47500000
## 52 52 0.1477986 0.1721310 0.9716667 0.87666667 0.87166667
## 53 53 0.3585870 0.2428305 0.3116667 0.19500000 0.14500000
## 54 54 0.2791049 0.2035791 0.5700000 0.43500000 0.38666667
## 55 55 0.2358034 0.2225726 0.7366667 0.56500000 0.52833333
## 56 56 0.2371105 0.1704332 0.7683333 0.64666667 0.62166667
## 57 57 0.3090121 0.2853253 0.4633333 0.29000000 0.22666667
## 58 58 0.2552678 0.2170866 0.6450000 0.51666667 0.47333333
## 59 59 0.2116554 0.2069662 0.8150000 0.66166667 0.63666667
## 60 60 0.2392023 0.1600761 0.7516667 0.66166667 0.63500000# ---------------------------
# 6. PLOTS (save to files & keep in environment)
# ---------------------------
cat("Creating plots... (files will be saved in working directory)\n")## Creating plots... (files will be saved in working directory)# Main: distributions of ratios (AUC & Cmax)
p1 <- ggplot(main_results, aes(x = ratio_cmax)) +
geom_histogram(binwidth = 0.02, fill = "skyblue", color = "black") +
geom_vline(xintercept = c(BE_lower, BE_upper), linetype = "dashed", color = "red") +
labs(title = "Cmax Ratio Distribution (Main Scenario)", x = "Cmax Ratio", y = "Count")
p2 <- ggplot(main_results, aes(x = ratio_auc)) +
geom_histogram(binwidth = 0.02, fill = "lightcoral", color = "black") +
geom_vline(xintercept = c(BE_lower, BE_upper), linetype = "dashed", color = "red") +
labs(title = "AUC Ratio Distribution (Main Scenario)", x = "AUC Ratio", y = "Count")
p1
p2
# Heatmap (Cmax)
library(reshape2)
heat_cmax <- dcast(heatmap_rows, n ~ true_GMR, value.var = "prob_cmax")
# convert to plotting-friendly long form
heat_long_cmax <- heatmap_rows[, c("true_GMR", "n", "prob_cmax")]
p_heat_cmax <- ggplot(heat_long_cmax, aes(x = factor(true_GMR), y = factor(n), fill = prob_cmax)) +
geom_tile() + scale_fill_gradient(low = "white", high = "steelblue") +
labs(x = "True GMR", y = "Sample Size per Arm", fill = "P(BE) (Cmax)",
title = "Heatmap: P(BE) for Cmax vs True GMR & Sample Size") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_heat_cmax
# Heatmap (AUC)
heat_long_auc <- heatmap_rows[, c("true_GMR", "n", "prob_auc")]
p_heat_auc <- ggplot(heat_long_auc, aes(x = factor(true_GMR), y = factor(n), fill = prob_auc)) +
geom_tile() + scale_fill_gradient(low = "white", high = "tomato") +
labs(x = "True GMR", y = "Sample Size per Arm", fill = "P(BE) (AUC)",
title = "Heatmap: P(BE) for AUC vs True GMR & Sample Size") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_heat_auc
# Parameter uncertainty plots
p_param_cv_between <- ggplot(param_results, aes(x = CV_between, y = prob_pass_both)) +
geom_point() + geom_smooth(method = "loess", se = TRUE) +
labs(title = "Effect of CV_between on P(BE both)", x = "CV_between", y = "P(pass both)")
p_param_cv_resid <- ggplot(param_results, aes(x = CV_resid, y = prob_pass_both)) +
geom_point() + geom_smooth(method = "loess", se = TRUE) +
labs(title = "Effect of CV_resid on P(BE both)", x = "CV_resid", y = "P(pass both)")
p_param_cv_between## `geom_smooth()` using formula = 'y ~ x'
p_param_cv_resid## `geom_smooth()` using formula = 'y ~ x'
# ggsave("param_cv_between_vs_prob.png", p_param_cv_between, width = 7, height = 4)
# ggsave("param_cv_resid_vs_prob.png", p_param_cv_resid, width = 7, height = 4)