Performance Comparison Guide

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)\)

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.

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.

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

type1_matrix <- run_simulation("homoskedastic")
type1_rates <- calc_rejection_rate(type1_matrix)
type1_rates
#>         white breusch_pagan       koenker        harvey 
#>         0.045         0.060         0.065         0.045

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()
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.

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.

linear_matrix <- run_simulation("linear")
#> Warning: Residual outliers detected at rows 87 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 87 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 87 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 87 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 62 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 62 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 62 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 62 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 39 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 39 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 39 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 39 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 90 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 90 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 90 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 90 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 122 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 122 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 122 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 122 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 11 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 11 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 11 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 11 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 130 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 130 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 130 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 130 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 51 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 51 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 51 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 51 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 137 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 137 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 137 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 137 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 14 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 14 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 14 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 14 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 107 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 107 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 107 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 107 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 145 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 145 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 145 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 145 (|z| > 5). Inspect leverage
#> before running Harvey.
cluster_matrix <- run_simulation("clustered")
#> Warning: Residual outliers detected at rows 91 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 91 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 91 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 91 (|z| > 5). Inspect leverage
#> before running Harvey.
#> Warning: Residual outliers detected at rows 114 (|z| > 5). Inspect leverage
#> before running White.
#> Warning: Residual outliers detected at rows 114 (|z| > 5). Inspect leverage
#> before running Breusch-Pagan.
#> Warning: Residual outliers detected at rows 114 (|z| > 5). Inspect leverage
#> before running Koenker.
#> Warning: Residual outliers detected at rows 114 (|z| > 5). Inspect leverage
#> before running Harvey.
linear_rates <- calc_rejection_rate(linear_matrix)
cluster_rates <- calc_rejection_rate(cluster_matrix)

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.

Performance Comparison Guide

Simulation design

Type I error control

Power against structured alternatives

Practical guidance