Monte Carlo Study: Power of the Heterogeneity Test
The worked examples show the cascade on individual datasets. We now assess the statistical power of the shape heterogeneity test (step 1) as a function of how different the component shapes actually are.
Design
We parameterize shape heterogeneity by the coefficient of variation (CV) of the shape parameters. All \(m = 5\) shapes are set to \(k_j = \bar{k}(1 + \text{cv} \cdot z_j)\) where \(\bar{k} = 1.5\) and \(z_1, \ldots, z_m\) are fixed offsets with mean 0 and standard deviation 1.
| CV | Meaning |
|---|---|
| 0% | All shapes equal — homogeneous model is correct |
| 5–10% | Mild heterogeneity — hard to detect |
| 20–30% | Moderate heterogeneity — detectable with adequate \(n\) |
| 50% | Strong heterogeneity — reliably detected |
Other settings: \(n = 500\) observations, \(p = 0.2\) (masking probability), right-censoring at the 80th percentile of the system lifetime under the homogeneous baseline, \(\alpha = 0.05\).
set.seed(42)
m <- 5
base_k <- 1.5
base_scales <- c(300, 400, 500, 600, 700)
cv_levels <- c(0, 0.05, 0.10, 0.20, 0.30, 0.50)
n_reps <- 200
n_obs <- 500
p_mask <- 0.2
alpha <- 0.05
# Fixed offsets: mean 0, sd 1
offsets <- c(-2, -1, 0, 1, 2)
offsets <- offsets / sd(offsets)
# Right-censoring time (based on homogeneous baseline at CV = 0)
beta_sys <- wei_series_system_scale(base_k, base_scales)
tau <- qweibull(0.8, shape = base_k, scale = beta_sys)
# Models and solvers
model_exp <- exp_series_md_c1_c2_c3()
model_hom <- wei_series_homogeneous_md_c1_c2_c3()
model_het <- wei_series_md_c1_c2_c3()
solver_exp <- fit(model_exp)
solver_hom <- fit(model_hom)
solver_het <- fit(model_het)
gen_het <- rdata(model_het)
np_exp <- m
np_hom <- m + 1
np_het <- 2 * mresults_list <- vector("list", length(cv_levels) * n_reps)
idx <- 1
for (cv in cv_levels) {
shapes <- base_k * (1 + cv * offsets)
theta <- as.numeric(rbind(shapes, base_scales))
for (r in seq_len(n_reps)) {
df_r <- gen_het(theta, n = n_obs, p = p_mask,
observe = observe_right_censor(tau = tau))
row <- list(cv = cv, rep = r,
reject_het_vs_hom = NA, reject_hom_vs_exp = NA,
ll_exp = NA, ll_hom = NA, ll_het = NA,
aic_exp = NA, aic_hom = NA, aic_het = NA,
bic_exp = NA, bic_hom = NA, bic_het = NA,
converged_all = FALSE)
tryCatch({
fe <- solver_exp(df_r, par = rep(0.002, m), method = "Nelder-Mead")
# Two-phase for hom: Nelder-Mead warm-up -> L-BFGS-B refinement
fh_nm <- solver_hom(df_r, par = c(1, rep(500, m)),
method = "Nelder-Mead",
control = list(maxit = 5000))
fh <- solver_hom(df_r, par = fh_nm$par,
method = "L-BFGS-B", lower = rep(1e-6, np_hom))
# Two-phase for het: warm-start from hom, Nelder-Mead -> L-BFGS-B
het_start <- as.numeric(rbind(rep(fh$par[1], m), fh$par[-1]))
fw_nm <- solver_het(df_r, par = het_start,
method = "Nelder-Mead",
control = list(maxit = 10000))
fw <- solver_het(df_r, par = fw_nm$par,
method = "L-BFGS-B", lower = rep(1e-6, np_het))
if (fe$converged && fh$converged && fw$converged) {
row$ll_exp <- fe$loglik
row$ll_hom <- fh$loglik
row$ll_het <- fw$loglik
row$aic_exp <- -2 * fe$loglik + 2 * np_exp
row$aic_hom <- -2 * fh$loglik + 2 * np_hom
row$aic_het <- -2 * fw$loglik + 2 * np_het
row$bic_exp <- -2 * fe$loglik + log(n_obs) * np_exp
row$bic_hom <- -2 * fh$loglik + log(n_obs) * np_hom
row$bic_het <- -2 * fw$loglik + log(n_obs) * np_het
LRT_het <- -2 * (fh$loglik - fw$loglik)
LRT_hom <- -2 * (fe$loglik - fh$loglik)
row$reject_het_vs_hom <- pchisq(LRT_het, df = m - 1,
lower.tail = FALSE) < alpha
row$reject_hom_vs_exp <- pchisq(LRT_hom, df = 1,
lower.tail = FALSE) < alpha
row$converged_all <- TRUE
}
}, error = function(e) NULL)
results_list[[idx]] <- row
idx <- idx + 1
}
cat(sprintf("CV = %.0f%% complete\n", cv * 100))
}
power_study <- do.call(rbind, lapply(results_list, as.data.frame))Rejection rates
ps <- mc_data$power_study
ps_conv <- ps[ps$converged_all, ]
rej_summary <- do.call(rbind, lapply(split(ps_conv, ps_conv$cv), function(d) {
data.frame(
CV_pct = d$cv[1] * 100,
n_converged = nrow(d),
reject_het_vs_hom = mean(d$reject_het_vs_hom),
reject_hom_vs_exp = mean(d$reject_hom_vs_exp)
)
}))
knitr::kable(rej_summary, digits = 3, row.names = FALSE,
caption = sprintf(
"Rejection rates at alpha = %.2f (n = %d, m = %d)",
mc_data$config$alpha, mc_data$config$n, mc_data$config$m),
col.names = c("Shape CV (%)", "Converged",
"Reject Hom (Het vs Hom)",
"Reject Exp (Hom vs Exp)"))| Shape CV (%) | Converged | Reject Hom (Het vs Hom) | Reject Exp (Hom vs Exp) |
|---|---|---|---|
| 0 | 200 | 0.045 | 1.00 |
| 5 | 200 | 0.050 | 1.00 |
| 10 | 200 | 0.145 | 1.00 |
| 20 | 200 | 0.350 | 1.00 |
| 30 | 200 | 0.860 | 0.78 |
| 50 | 200 | 1.000 | 0.99 |
Power curve
plot(rej_summary$CV_pct, rej_summary$reject_het_vs_hom,
type = "b", pch = 19, col = "steelblue", lwd = 2,
xlab = "Coefficient of Variation of Shapes (%)",
ylab = "Rejection Rate",
ylim = c(0, 1),
main = "Power of the LRT Cascade")
lines(rej_summary$CV_pct, rej_summary$reject_hom_vs_exp,
type = "b", pch = 17, col = "firebrick", lwd = 2, lty = 2)
abline(h = mc_data$config$alpha, lty = 3, col = "gray50")
legend("right",
c("Het vs Hom (shape heterogeneity)",
"Hom vs Exp (aging)",
expression(alpha == 0.05)),
col = c("steelblue", "firebrick", "gray50"),
lty = c(1, 2, 3), pch = c(19, 17, NA), lwd = 2, cex = 0.85)
grid()
At CV = 0% (no heterogeneity), the rejection rate of the shape heterogeneity test should be near the nominal level \(\alpha = 0.05\) — the Type I error rate. As CV increases, the power rises, eventually reaching near 1 at CV = 50%.
The aging test (hom vs exp) maintains high rejection rates across all CV levels, because the baseline shape \(\bar{k} = 1.5\) is far from 1 regardless of how much the individual shapes vary around it.
AIC/BIC model selection accuracy
We also evaluate how often each information criterion selects the correct model. At CV = 0, the correct model is homogeneous Weibull. At CV > 0, it is heterogeneous Weibull.
aic_correct <- do.call(rbind, lapply(split(ps_conv, ps_conv$cv), function(d) {
cv <- d$cv[1]
aic_choice <- ifelse(d$aic_het < d$aic_hom & d$aic_het < d$aic_exp, "het",
ifelse(d$aic_hom < d$aic_exp, "hom", "exp"))
bic_choice <- ifelse(d$bic_het < d$bic_hom & d$bic_het < d$bic_exp, "het",
ifelse(d$bic_hom < d$bic_exp, "hom", "exp"))
correct_model <- if (cv == 0) "hom" else "het"
data.frame(
CV_pct = cv * 100,
AIC_correct = mean(aic_choice == correct_model),
BIC_correct = mean(bic_choice == correct_model)
)
}))
knitr::kable(aic_correct, digits = 3, row.names = FALSE,
caption = "Proportion selecting the correct model",
col.names = c("Shape CV (%)", "AIC Correct", "BIC Correct"))| Shape CV (%) | AIC Correct | BIC Correct |
|---|---|---|
| 0 | 0.925 | 1.000 |
| 5 | 0.080 | 0.000 |
| 10 | 0.210 | 0.000 |
| 20 | 0.445 | 0.025 |
| 30 | 0.920 | 0.405 |
| 50 | 1.000 | 1.000 |