## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.align = "center",
  comment = "",
  fig.width = 10,
  fig.height = 5,
  dpi = 150
)

# SEED STRATEGY:
# - Global seed (2024) for reproducibility
# - Local seeds (123, 456, etc.) used for specific data visualizations
# - The Monte Carlo study is precomputed and seeded in
#   data-raw/vignette-simulations.R
set.seed(2024)

# Standardized color palette for all comparisons
MODEL_COLORS <- c(
  "Beta (betareg)" = "#D32F2F", # Red
  "Beta (gkwreg)" = "#1976D2", # Blue
  "Kumaraswamy" = "#388E3C", # Green
  "Exp. Kumaraswamy" = "#7B1FA2" # Purple
)

MODEL_NAMES <- c(
  betareg = "Beta (betareg)",
  gkw_beta = "Beta (gkwreg)",
  gkw_kw = "Kumaraswamy",
  gkw_ekw = "Exp. Kumaraswamy"
)

## ----load_precomputed, include=FALSE------------------------------------------
# The Monte Carlo study below (3 scenarios x 200 replications x 4 models) is
# precomputed. Running it live made the vignette take about 11 minutes to
# rebuild on the CRAN check machines, which is more than CRAN can afford to run
# regularly. The study itself is unchanged -- still 200 replications -- but the
# summaries are read from a stored object instead of being recomputed here.
#
# Regenerate with: Rscript data-raw/vignette-simulations.R
sim_results <- readRDS(
  system.file("extdata", "vignette-simulations.rds", package = "gkwreg")
)

# Quantities quoted in the prose below. They are derived from the tables rather
# than written out by hand, so the narrative cannot drift away from the results
# when the study is regenerated.
.pick <- function(tab, model, col) tab[[col]][tab$Model == model]
.speedup <- function(tab) {
  .pick(tab, "Beta (betareg)", "Time") / .pick(tab, "Kumaraswamy", "Time")
}

speedup_s1 <- .speedup(sim_results$comp_s1)
speedup_s2 <- .speedup(sim_results$comp_s2)
speedup_s3 <- .speedup(sim_results$comp_s3)
speedup_lo <- min(speedup_s1, speedup_s2, speedup_s3)
speedup_hi <- max(speedup_s1, speedup_s2, speedup_s3)

conv_beta_s3 <- .pick(sim_results$comp_s3, "Beta (betareg)", "Conv_Rate")
conv_kw_s3 <- .pick(sim_results$comp_s3, "Kumaraswamy", "Conv_Rate")
n_conv_beta_s3 <- round(conv_beta_s3 / 100 * sim_results$n_sim)

aic_beta_s1 <- .pick(sim_results$comp_s1, "Beta (betareg)", "AIC")
aic_kw_s1 <- .pick(sim_results$comp_s1, "Kumaraswamy", "AIC")

fmt_x <- function(x) sprintf("%.1f", x)

## ----packages-----------------------------------------------------------------
# betareg is only in Suggests: skip the whole vignette when unavailable
# (CRAN policy on the conditional use of suggested packages).
if (!requireNamespace("betareg", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
  knitr::asis_output(
    "> **Note**: this vignette requires the **betareg** package, which is not\n> installed. All code below is shown but not evaluated.\n"
  )
} else {
  library(betareg)
}
library(gkwreg)
library(ggplot2)

## ----scenario1_dgp------------------------------------------------------------
dgp_beta <- function(n, params) {
  x1 <- rnorm(n, 0, 1)
  x2 <- runif(n, -1, 1)

  eta_mu <- params$beta_mu[1] + params$beta_mu[2] * x1 + params$beta_mu[3] * x2
  eta_phi <- params$beta_phi[1] + params$beta_phi[2] * x1

  mu <- plogis(eta_mu)
  phi <- exp(eta_phi)

  y <- rbeta(n, mu * phi, (1 - mu) * phi)
  y <- pmax(pmin(y, 0.999), 0.001)

  data.frame(y = y, x1 = x1, x2 = x2, mu = mu, phi = phi)
}

## ----scenario1_visualization, fig.height=6, fig.cap="Distributional Characteristics - Scenario 1 (Well-Specified Beta)"----
set.seed(123)
vis_data_s1 <- dgp_beta(1000, list(
  beta_mu = c(0.5, -0.8, 0.6),
  beta_phi = c(1.5, 0.4)
))

par(mfrow = c(2, 2), mar = c(4, 4, 2, 2))

# Response distribution
hist(vis_data_s1$y,
  breaks = 40, col = "#388E3C", border = "white",
  main = "Response Distribution", xlab = "Y", prob = TRUE,
  cex.main = 1.2, cex.lab = 1.1
)
lines(density(vis_data_s1$y), col = "#D32F2F", lwd = 2)
legend("topright", "Kernel density", col = "#D32F2F", lwd = 2, bty = "n")

# Effect of x1
plot(vis_data_s1$x1, vis_data_s1$y,
  pch = 16, col = rgb(0, 0, 0, 0.3),
  main = "Covariate Effect: x1", xlab = "x1", ylab = "Y",
  cex.main = 1.2, cex.lab = 1.1
)
lines(lowess(vis_data_s1$x1, vis_data_s1$y), col = "#D32F2F", lwd = 3)
abline(h = 0.5, lty = 2, col = "gray50")

# Effect of x2
plot(vis_data_s1$x2, vis_data_s1$y,
  pch = 16, col = rgb(0, 0, 0, 0.3),
  main = "Covariate Effect: x2", xlab = "x2", ylab = "Y",
  cex.main = 1.2, cex.lab = 1.1
)
lines(lowess(vis_data_s1$x2, vis_data_s1$y), col = "#388E3C", lwd = 3)
abline(h = 0.5, lty = 2, col = "gray50")

# Conditional variance
plot(abs(vis_data_s1$x1), (vis_data_s1$y - vis_data_s1$mu)^2,
  pch = 16, col = rgb(0, 0, 0, 0.2),
  main = "Heteroscedasticity Structure",
  xlab = "|x1|", ylab = "Squared Residuals",
  cex.main = 1.2, cex.lab = 1.1
)
lines(lowess(abs(vis_data_s1$x1), (vis_data_s1$y - vis_data_s1$mu)^2),
  col = "#1976D2", lwd = 3
)

## ----scenario1_sim------------------------------------------------------------
# Precomputed; see data-raw/vignette-simulations.R
comp_s1 <- sim_results$comp_s1

## ----scenario1_coef_comparison------------------------------------------------
# Fit models to one dataset for coefficient comparison
set.seed(456)
data_coef_s1 <- dgp_beta(300, list(
  beta_mu = c(0.5, -0.8, 0.6),
  beta_phi = c(1.5, 0.4)
))

fit_beta_s1 <- betareg(y ~ x1 + x2 | x1, data = data_coef_s1)
fit_kw_s1 <- gkwreg(y ~ x1 + x2 | x1, data = data_coef_s1, family = "kw")

# Extract coefficients
coef_beta_s1 <- coef(fit_beta_s1)
se_beta_s1 <- sqrt(diag(vcov(fit_beta_s1)))

coef_kw_s1 <- coef(fit_kw_s1)
se_kw_s1 <- fit_kw_s1$se

# Build comparison table
coef_names <- names(coef_beta_s1)
true_values <- c(0.5, -0.8, 0.6, 1.5, 0.4)

coef_comparison_s1 <- data.frame(
  Parameter = coef_names,
  True_Value = true_values,
  Beta_Est = round(coef_beta_s1, 4),
  Beta_SE = round(se_beta_s1, 4),
  Beta_Z = round(coef_beta_s1 / se_beta_s1, 2),
  Kw_Est = round(coef_kw_s1, 4),
  Kw_SE = round(se_kw_s1, 4),
  Kw_Z = round(coef_kw_s1 / se_kw_s1, 2)
)

# Calculate differences
coef_comparison_s1$Bias_Beta <- round(coef_comparison_s1$Beta_Est - true_values, 4)
coef_comparison_s1$Bias_Kw <- round(coef_comparison_s1$Kw_Est - true_values, 4)
coef_comparison_s1$SE_Ratio <- round(coef_comparison_s1$Kw_SE / coef_comparison_s1$Beta_SE, 3)

knitr::kable(
  coef_comparison_s1,
  row.names = FALSE,
  caption = "Table 1: Parameter Estimates - Scenario 1 (Well-Specified Beta, n=300)",
  col.names = c(
    "Parameter", "True", "Est.", "SE", "z", "Est.", "SE", "z",
    "Bias", "Bias", "SE Ratio"
  ),
  align = c("l", rep("r", 10)),
  digits = 4
)

## ----scenario2_dgp------------------------------------------------------------
dgp_heavy_tails <- function(n, params) {
  x1 <- rnorm(n, 0, 1)
  x2 <- rbinom(n, 1, 0.5)

  eta_alpha <- params$beta_alpha[1] + params$beta_alpha[2] * x1
  eta_beta <- params$beta_beta[1] + params$beta_beta[2] * x2
  eta_lambda <- params$beta_lambda[1]

  alpha <- exp(eta_alpha)
  beta <- exp(eta_beta)
  lambda <- exp(eta_lambda)

  u <- runif(n)
  y <- (1 - (1 - u^(1 / beta))^(1 / alpha))^(1 / lambda)
  y <- pmax(pmin(y, 0.9999), 0.0001)

  data.frame(y = y, x1 = x1, x2 = factor(x2))
}

## ----scenario2_visualization, fig.height=6, fig.cap="Distributional Characteristics - Scenario 2 (Heavy Tails)"----
set.seed(789)
vis_data_s2 <- dgp_heavy_tails(1500, list(
  beta_alpha = c(0.8, -0.5),
  beta_beta = c(0.3, 0.4),
  beta_lambda = c(0.6)
))

par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))

# Response distribution
hist(vis_data_s2$y,
  breaks = 60, col = "#7B1FA2", border = "white",
  main = "Heavy-Tailed Distribution", xlab = "Y", prob = TRUE,
  cex.main = 1.2, cex.lab = 1.1
)
lines(density(vis_data_s2$y), col = "#D32F2F", lwd = 2)
# Add Beta approximation for comparison
beta_approx <- rbeta(1500, 2, 2)
lines(density(beta_approx), col = "#1976D2", lwd = 2, lty = 2)
legend("topright", c("True density", "Beta approximation"),
  col = c("#D32F2F", "#1976D2"), lwd = 2, lty = c(1, 2), bty = "n", cex = 0.9
)

# QQ plot against Beta
qqplot(rbeta(1500, 2, 2), vis_data_s2$y,
  main = "Q-Q Plot: EKw vs Beta(2,2)",
  xlab = "Beta(2,2) Quantiles", ylab = "Observed Quantiles",
  pch = 16, col = rgb(0, 0, 0, 0.3), cex.main = 1.2, cex.lab = 1.1
)
abline(0, 1, col = "#D32F2F", lwd = 2)
text(0.2, 0.8, "Heavier tails\n(departures)", col = "#D32F2F", cex = 0.9)

# Effect by group
boxplot(y ~ x2,
  data = vis_data_s2, col = c("#388E3C", "#7B1FA2"),
  main = "Distribution by Group", xlab = "x2", ylab = "Y",
  names = c("Group 0", "Group 1"), cex.main = 1.2, cex.lab = 1.1
)

# Tail behavior
plot(sort(vis_data_s2$y)[1:50],
  ylab = "Y", xlab = "Order (Lower Tail)",
  main = "Lower Tail Concentration", pch = 16, col = "#7B1FA2",
  cex.main = 1.2, cex.lab = 1.1
)
beta_tail <- sort(rbeta(1500, 2, 2))[1:50]
points(beta_tail, pch = 1, col = "#1976D2")
legend("topleft", c("EKw", "Beta(2,2)"),
  col = c("#7B1FA2", "#1976D2"),
  pch = c(16, 1), bty = "n"
)

## ----scenario2_sim------------------------------------------------------------
# Precomputed; see data-raw/vignette-simulations.R
comp_s2 <- sim_results$comp_s2

## ----scenario2_coef_comparison------------------------------------------------
# Fit models to one dataset
set.seed(101112)
data_coef_s2 <- dgp_heavy_tails(300, list(
  beta_alpha = c(0.8, -0.5),
  beta_beta = c(0.3, 0.4),
  beta_lambda = c(0.6)
))

fit_beta_s2 <- betareg(y ~ x1 | x2, data = data_coef_s2)
fit_kw_s2 <- gkwreg(y ~ x1 | x2, data = data_coef_s2, family = "kw")

coef_beta_s2 <- coef(fit_beta_s2)
se_beta_s2 <- sqrt(diag(vcov(fit_beta_s2)))

coef_kw_s2 <- coef(fit_kw_s2)
se_kw_s2 <- fit_kw_s2$se

# Note: True parameters are on different scale (α, β vs μ, φ)
# Focus on relative comparisons
coef_comparison_s2 <- data.frame(
  Parameter = names(coef_beta_s2),
  Beta_Est = round(coef_beta_s2, 4),
  Beta_SE = round(se_beta_s2, 4),
  Beta_Z = round(coef_beta_s2 / se_beta_s2, 2),
  Kw_Est = round(coef_kw_s2, 4),
  Kw_SE = round(se_kw_s2, 4),
  Kw_Z = round(coef_kw_s2 / se_kw_s2, 2),
  SE_Ratio = round(se_kw_s2 / se_beta_s2, 3)
)

knitr::kable(
  coef_comparison_s2,
  row.names = FALSE,
  caption = "Table 2: Parameter Estimates - Scenario 2 (Heavy Tails, n=300)",
  col.names = c("Parameter", "Est.", "SE", "z", "Est.", "SE", "z", "SE Ratio"),
  align = c("l", rep("r", 7)),
  digits = 4
)

## ----scenario3_dgp------------------------------------------------------------
dgp_extreme <- function(n, params) {
  x1 <- rnorm(n, 0, 1)
  group <- sample(c("J", "U"), n, replace = TRUE)

  alpha <- ifelse(
    group == "J",
    exp(params$alpha_J[1] + params$alpha_J[2] * x1),
    exp(params$alpha_U[1] + params$alpha_U[2] * x1)
  )

  beta <- ifelse(group == "J", exp(params$beta_J), exp(params$beta_U))

  u <- runif(n)
  y <- (1 - (1 - u)^(1 / beta))^(1 / alpha)
  y <- pmax(pmin(y, 0.9999), 0.0001)

  data.frame(y = y, x1 = x1, group = factor(group))
}

## ----scenario3_visualization, fig.height=7, fig.cap="Extreme Distributional Shapes - Scenario 3 (Boundary Concentration)"----
set.seed(131415)
vis_data_s3 <- dgp_extreme(2000, list(
  alpha_J = c(-1.5, 0.2),
  beta_J = 2.0,
  alpha_U = c(-1.8, 0.1),
  beta_U = -0.8
))

par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))

# J-shaped distribution
hist(vis_data_s3$y[vis_data_s3$group == "J"],
  breaks = 50,
  col = "#FF6F00", border = "white",
  main = "J-Shaped Distribution", xlab = "Y", prob = TRUE,
  xlim = c(0, 1), cex.main = 1.2, cex.lab = 1.1
)
lines(density(vis_data_s3$y[vis_data_s3$group == "J"]),
  col = "#D32F2F", lwd = 2
)
text(0.5, 15, sprintf(
  "80%% of mass\n< %.2f",
  quantile(vis_data_s3$y[vis_data_s3$group == "J"], 0.8)
),
col = "#D32F2F", cex = 0.9
)

# U-shaped distribution
hist(vis_data_s3$y[vis_data_s3$group == "U"],
  breaks = 50,
  col = "#0277BD", border = "white",
  main = "U-Shaped Distribution", xlab = "Y", prob = TRUE,
  xlim = c(0, 1), cex.main = 1.2, cex.lab = 1.1
)
lines(density(vis_data_s3$y[vis_data_s3$group == "U"]),
  col = "#D32F2F", lwd = 2
)
text(0.5, 3, "Bimodal:\nconcentration\nat boundaries",
  col = "#D32F2F", cex = 0.9
)

# Combined distribution
plot(density(vis_data_s3$y),
  main = "Mixture Distribution",
  xlab = "Y", ylab = "Density", lwd = 2, col = "#7B1FA2",
  cex.main = 1.2, cex.lab = 1.1
)
polygon(density(vis_data_s3$y), col = rgb(0.48, 0.12, 0.63, 0.3), border = NA)

# Empirical CDF showing boundary concentration
plot(ecdf(vis_data_s3$y),
  main = "Empirical CDF",
  xlab = "Y", ylab = "F(Y)", lwd = 2, col = "#388E3C",
  cex.main = 1.2, cex.lab = 1.1
)
abline(v = c(0.1, 0.9), lty = 2, col = "gray40")
text(0.05, 0.5, sprintf(
  "%.0f%%\n< 0.1",
  100 * mean(vis_data_s3$y < 0.1)
),
col = "#D32F2F", cex = 0.9
)
text(0.95, 0.5, sprintf(
  "%.0f%%\n> 0.9",
  100 * mean(vis_data_s3$y > 0.9)
),
col = "#D32F2F", cex = 0.9
)

## ----scenario3_sim------------------------------------------------------------
# Precomputed; see data-raw/vignette-simulations.R
comp_s3 <- sim_results$comp_s3

## ----scenario3_coef_comparison------------------------------------------------
# For this scenario, we'll fit to a single J-shaped dataset where Beta converges
set.seed(161718)
# Generate J-shaped only for better Beta convergence probability
data_coef_s3 <- dgp_extreme(400, list(
  alpha_J = c(-1.2, 0.15), # Slightly less extreme for illustration
  beta_J = 1.8,
  alpha_U = c(-1.5, 0.1),
  beta_U = -0.5
))

# Try fitting - Beta may fail
fit_beta_s3 <- tryCatch(
  {
    betareg(y ~ x1 * group | group, data = data_coef_s3)
  },
  error = function(e) NULL
)

fit_kw_s3 <- gkwreg(y ~ x1 * group | group, data = data_coef_s3, family = "kw")

if (!is.null(fit_beta_s3) && fit_beta_s3$converged) {
  coef_beta_s3 <- coef(fit_beta_s3)
  se_beta_s3 <- sqrt(diag(vcov(fit_beta_s3)))

  coef_kw_s3 <- coef(fit_kw_s3)
  se_kw_s3 <- fit_kw_s3$se

  coef_comparison_s3 <- data.frame(
    Parameter = names(coef_beta_s3),
    Beta_Est = round(coef_beta_s3, 4),
    Beta_SE = round(se_beta_s3, 4),
    Beta_Z = round(coef_beta_s3 / se_beta_s3, 2),
    Kw_Est = round(coef_kw_s3, 4),
    Kw_SE = round(se_kw_s3, 4),
    Kw_Z = round(coef_kw_s3 / se_kw_s3, 2),
    SE_Ratio = round(se_kw_s3 / se_beta_s3, 3)
  )

  knitr::kable(
    coef_comparison_s3,
    row.names = FALSE,
    caption = "Table 3: Parameter Estimates - Scenario 3 (Extreme Shapes, n=400, Beta Converged)",
    col.names = c("Parameter", "Est.", "SE", "z", "Est.", "SE", "z", "SE Ratio"),
    align = c("l", rep("r", 7)),
    digits = 4
  )
} else {
  cat("**Table 3: Parameter Estimates - Scenario 3**\n\n")
  cat("Beta regression failed to converge (as expected in 94.5% of cases).\n")
  cat("Only Kumaraswamy estimates are available:\n\n")

  coef_kw_s3 <- coef(fit_kw_s3)
  se_kw_s3 <- fit_kw_s3$se

  knitr::kable(
    data.frame(
      Parameter = names(coef_kw_s3),
      Estimate = round(coef_kw_s3, 4),
      SE = round(se_kw_s3, 4),
      z_stat = round(coef_kw_s3 / se_kw_s3, 2),
      p_value = round(2 * pnorm(-abs(coef_kw_s3 / se_kw_s3)), 4)
    ),
    row.names = FALSE,
    caption = "Table 3: Kumaraswamy Parameter Estimates - Scenario 3 (Beta Failed to Converge)",
    col.names = c("Parameter", "Estimate", "SE", "z", "p-value"),
    align = c("l", rep("r", 4)),
    digits = 4
  )
}

## ----results_table, echo=FALSE------------------------------------------------
all_scenarios <- rbind(
  data.frame(Scenario = "S1: Well-Specified Beta", comp_s1),
  data.frame(Scenario = "S2: Heavy Tails", comp_s2),
  data.frame(Scenario = "S3: Extreme Shapes", comp_s3)
)

knitr::kable(
  all_scenarios,
  row.names = FALSE,
  caption = "Table 4: Comprehensive Model Comparison Across Three Simulation Scenarios. Note: Scenario 3 Beta statistics are based on only `r conv_beta_s3`% successful fits; metrics may not be representative.",
  digits = c(0, 0, 0, 1, 2, 4, 4),
  align = c("l", "l", "r", "r", "r", "r", "r")
)

## ----comprehensive_comparison, echo=FALSE, fig.height=6, fig.cap="Comparative Performance Across Scenarios"----
par(mfrow = c(1, 3), mar = c(10, 4, 3, 1))

# AIC by scenario for Beta and Kumaraswamy only
scenarios <- c("S1", "S2", "S3")
beta_aic <- c(
  comp_s1$AIC[comp_s1$Model == "Beta (betareg)"],
  comp_s2$AIC[comp_s2$Model == "Beta (betareg)"],
  comp_s3$AIC[comp_s3$Model == "Beta (betareg)"]
)
kw_aic <- c(
  comp_s1$AIC[comp_s1$Model == "Kumaraswamy"],
  comp_s2$AIC[comp_s2$Model == "Kumaraswamy"],
  comp_s3$AIC[comp_s3$Model == "Kumaraswamy"]
)

# Handle NAs and extreme values for S3
# If Beta failed, use Kw AIC * 1.5 for visualization purposes
if (is.na(beta_aic[3]) || !is.finite(beta_aic[3])) {
  beta_aic[3] <- kw_aic[3] * 1.5
}

aic_matrix <- rbind(beta_aic, kw_aic)
colnames(aic_matrix) <- scenarios

barplot(aic_matrix,
  beside = TRUE, col = c("#D32F2F", "#388E3C"),
  names.arg = scenarios, las = 1,
  main = "Model Fit (AIC)", ylab = "AIC",
  cex.names = 1.2, cex.main = 1.2
)
legend("topleft", c("Beta", "Kumaraswamy"),
  fill = c("#D32F2F", "#388E3C"),
  bty = "n", cex = 1.1
)
text(2, max(beta_aic) * 0.9, "Lower\nis better", cex = 0.9, col = "gray30")

# Convergence rates
beta_conv <- c(
  comp_s1$Conv_Rate[comp_s1$Model == "Beta (betareg)"],
  comp_s2$Conv_Rate[comp_s2$Model == "Beta (betareg)"],
  comp_s3$Conv_Rate[comp_s3$Model == "Beta (betareg)"]
)
kw_conv <- c(
  comp_s1$Conv_Rate[comp_s1$Model == "Kumaraswamy"],
  comp_s2$Conv_Rate[comp_s2$Model == "Kumaraswamy"],
  comp_s3$Conv_Rate[comp_s3$Model == "Kumaraswamy"]
)

conv_matrix <- rbind(beta_conv, kw_conv)
barplot(conv_matrix,
  beside = TRUE, col = c("#D32F2F", "#388E3C"),
  names.arg = scenarios, las = 1, ylim = c(0, 110),
  main = "Convergence Reliability", ylab = "Success Rate (%)",
  cex.names = 1.2, cex.main = 1.2
)
abline(h = 95, lty = 2, lwd = 2, col = "darkred")
text(1.5, 98, "95% threshold", col = "darkred", cex = 0.9)

# Computational speedup
beta_time <- c(
  comp_s1$Time[comp_s1$Model == "Beta (betareg)"],
  comp_s2$Time[comp_s2$Model == "Beta (betareg)"],
  comp_s3$Time[comp_s3$Model == "Beta (betareg)"]
)
kw_time <- c(
  comp_s1$Time[comp_s1$Model == "Kumaraswamy"],
  comp_s2$Time[comp_s2$Model == "Kumaraswamy"],
  comp_s3$Time[comp_s3$Model == "Kumaraswamy"]
)

# Handle NAs in time measurements
speedup <- beta_time / kw_time
# If speedup is NA or infinite, use 1 (no speedup) for visualization
speedup[!is.finite(speedup)] <- 1

barplot(speedup,
  names.arg = scenarios, col = "#388E3C",
  las = 1, main = "Computational Speedup",
  ylab = "Speedup Factor (Kw vs Beta)",
  cex.names = 1.2, cex.main = 1.2
)
abline(h = 1, lty = 2, col = "gray40")
text(1:3, speedup + 0.5, sprintf("%.1fx", speedup), cex = 1.1, font = 2)

## ----efficiency_analysis, echo=FALSE------------------------------------------
timing <- aggregate(
  Time ~ Model,
  data = all_scenarios[!is.na(all_scenarios$Time), ],
  FUN = mean
)
timing <- timing[order(timing$Time), ]
timing$Speedup <- timing$Time[timing$Model == "Beta (betareg)"] / timing$Time

knitr::kable(
  timing,
  row.names = FALSE,
  caption = "Table 5: Average Computational Time and Speedup Relative to Beta Regression",
  col.names = c("Model", "Mean Time (sec)", "Speedup Factor"),
  align = c("l", "r", "r"),
  digits = 4
)

## ----session_info, echo=FALSE-------------------------------------------------
sessionInfo()