---
title: "Performance Comparison Guide"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Performance Comparison Guide}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
library(heteroTests)
library(ggplot2)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
set.seed(2024)
```

Heteroscedasticity diagnostics differ in their sensitivity to alternative
variance structures and assumptions about the error distribution. This vignette
compares the Type I error rates and empirical power of several core tests using
a reproducible simulation design.

## Simulation design

We simulate data from a two-regressor linear model
\[
y_i = 1 + 2 x_{i1} + 3 x_{i2} + \varepsilon_i,
\]
with covariates drawn from $x_{i1} \sim \mathcal{N}(0,1)$ and
$x_{i2} = 0.5 x_{i1} + u_i$, $u_i \sim \mathcal{N}(0,1)$. The error variance is
controlled by the scenario:

* **Homoskedastic:** $\varepsilon_i \sim \mathcal{N}(0, 1)$
* **Linear heteroskedasticity:** $\varepsilon_i \sim \mathcal{N}\left(0,
  (1 + 1.5 |x_{i1}|)^2\right)$
* **Clustered heteroskedasticity:** with probability 0.5 the observation enters a
  high-variance regime $\mathcal{N}(0, 2.5^2)$, otherwise it remains at
  $\mathcal{N}(0, 1^2)$

```{r}
generate_sample <- function(n = 180, scenario = c("homoskedastic", "linear", "clustered")) {
  scenario <- match.arg(scenario)
  x1 <- rnorm(n)
  x2 <- 0.5 * x1 + rnorm(n)
  if (scenario == "homoskedastic") {
    errors <- rnorm(n, sd = 1)
  } else if (scenario == "linear") {
    errors <- rnorm(n, sd = 1 + 1.5 * abs(x1))
  } else {
    regime <- rbinom(n, 1, 0.5)
    sds <- ifelse(regime == 1, 2.5, 1)
    errors <- rnorm(n, sd = sds)
  }
  y <- 1 + 2 * x1 + 3 * x2 + errors
  data.frame(y = y, x1 = x1, x2 = x2)
}
```

To avoid simulation failures we wrap each diagnostic in a safe helper that
returns `NA` if input requirements are not met.

```{r}
safe_pvalue <- function(fn) {
  res <- tryCatch(fn(), error = function(e) NULL)
  if (is.null(res) || is.null(res$p.value)) {
    NA_real_
  } else {
    res$p.value
  }
}

run_diagnostics_once <- function(data) {
  model <- lm(y ~ x1 + x2, data = data)
  c(
    white = safe_pvalue(function() performWhiteTest(model, data)),
    breusch_pagan = safe_pvalue(function() performBPTest(model, data)),
    koenker = safe_pvalue(function() performKoenkerTest(model, data)),
    harvey = safe_pvalue(function() performHarveyTest(model))
  )
}
```

We repeat each scenario 200 times to estimate rejection frequencies at the 5%
level.

```{r}
run_simulation <- function(scenario, reps = 200) {
  replicate(reps, run_diagnostics_once(generate_sample(scenario = scenario)), simplify = "matrix")
}

calc_rejection_rate <- function(result_matrix, alpha = 0.05) {
  rowMeans(result_matrix < alpha, na.rm = TRUE)
}
```

## Type I error control

```{r}
type1_matrix <- run_simulation("homoskedastic")
type1_rates <- calc_rejection_rate(type1_matrix)
type1_rates
```

```{r fig.width=6, fig.height=4}
type1_df <- data.frame(
  test = names(type1_rates),
  rate = as.numeric(type1_rates),
  scenario = "Type I error"
)

ggplot(type1_df, aes(x = reorder(test, rate), y = rate)) +
  geom_col(fill = "#009E73", alpha = 0.85) +
  geom_hline(yintercept = 0.05, linetype = "dashed", colour = "#D55E00") +
  coord_cartesian(ylim = c(0, 0.15)) +
  coord_flip() +
  labs(
    x = "Test",
    y = "Empirical rejection rate",
    title = "Observed Type I error at 5% nominal level"
  ) +
  theme_minimal()
```

*Interpretation.* Breusch–Pagan and White tests slightly exceed the nominal rate
in moderately sized samples, while Koenker's robust statistic remains closest to
5%. Harvey's log-linear formulation holds size closely on its default variance
model; its `studentize = TRUE` variant is the safer choice under heavy-tailed errors.

## Power against structured alternatives

We evaluate power under the linear and clustered variance patterns.

```{r}
linear_matrix <- run_simulation("linear")
cluster_matrix <- run_simulation("clustered")
linear_rates <- calc_rejection_rate(linear_matrix)
cluster_rates <- calc_rejection_rate(cluster_matrix)
```

```{r fig.width=7, fig.height=4.5}
power_df <- rbind(
  data.frame(test = names(linear_rates), rate = as.numeric(linear_rates), scenario = "Linear variance"),
  data.frame(test = names(cluster_rates), rate = as.numeric(cluster_rates), scenario = "Clustered variance")
)

ggplot(power_df, aes(x = test, y = rate, fill = scenario)) +
  geom_col(position = "dodge", alpha = 0.85) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    x = "Test",
    y = "Detection probability",
    fill = "Scenario",
    title = "Empirical power across variance alternatives"
  ) +
  theme_minimal()
```

*Interpretation.* White's omnibus test excels when heteroskedasticity follows a
nonlinear pattern with multiple interaction terms (clustered scenario), whereas
Breusch–Pagan is most powerful for smooth linear variance inflation. Koenker's
studentised variant sacrifices some power for robustness. Harvey's
log-linear regression on the model regressors performs strongly when the
variance tracks those regressors, which occurs here because `x1` dominates the signal.

## Practical guidance

* Prefer Breusch–Pagan or Park when theory suggests variance scales with known
  regressors; Harvey targets the same regressors under a multiplicative rather
  than an additive variance model.
* Use Koenker's statistic or Glejser-type tests when heavy-tailed errors are
  plausible.
* Combine diagnostics with residual plots and domain knowledge; simulation
  studies complement but do not replace case-specific validation.
* For mission-critical decisions, calibrate $p$-values through bootstrap
  procedures (`performWhiteTestBootstrap()`) or resampling wrappers in the
  package.

The code above is fully reproducible and can be adapted to benchmark additional
tests available in heteroTests.