Poisson, Overdispersion, Hampel

code
python
Confusions
Author

Zhenglei Gao

Published

September 25, 2025

The Problem

  • The Hampel-style rule here is a simple threshold on sample variance minus mean; it’s scale-dependent and will often flag more in high-mean settings where sampling variability inflates var − mean slightly above 0.
  • The ratio rule uses \(\text{VMR}=\hat{\sigma}^2/\bar{y}\).
  • The Pearson test uses \(S=\sum (y_i-\bar y)^2/\bar y\) compared against \(\chi^2_{n-1}\) at 0.05.

Simulation

  • Scenario mu=10, var=16 (VMR=1.6), Negative Binomial: Strong overdispersion signal. VMR>1 flags most samples and increases with n (~0.78 at n=10 to ~1.00 at n=100). Pearson power rises with n (~0.28 to ~0.95). Hampel (var−mean>5) is moderately sensitive and increases with n (~0.47 to ~0.64).
  • Same scenario, Poisson(true): Pearson stays near 0.05 across n (as expected under the null). VMR>1 hovers around 0.44–0.48 regardless of n (descriptive, not a test). Hampel flags less as n grows (~0.14 at n=10 to ~0.00 at n=100).
  • Scenario mu=1000, var=1006 (VMR≈1.006): Poisson and NB look nearly identical. Pearson stays near 0.04–0.05. VMR>1 and Hampel both sit around 0.44–0.50 across n, reflecting variability and the fact that the Hampel threshold (5) is very small relative to the scale.
  • Sample means and variances in all cells align closely with their targets. Overall: Pearson controls type I error and gains power with n; VMR>1 is a descriptive flag that often triggers about half the time near-Poisson; the Hampel rule is scale-dependent and can overflag at high means.

Hampel Identifier Formula

The Hampel identifier is a heuristic rule based on the difference between the sample variance and the sample mean. It flags overdispersion if this difference exceeds a fixed threshold.

Let \(\bar{y}\) be the sample mean and \(\hat{\sigma}^2\) be the sample variance. The rule is:

Flag overdispersion if \(\hat{\sigma}^2 - \bar{y} > c\)

where \(c\) is a predetermined constant. In the CPCAT scenarios, \(c=5\).

Simulation to Compare Methods

Here is an R script that implements the simulation. It compares the Hampel Identifier, the Variance-to-Mean Ratio (VMR), and the Pearson Chi-squared Test across the specified scenarios.

# Load necessary libraries 
# library(dplyr)

# --- Simulation Setup ---

# Number of simulation repetitions for each scenario
N_sim <- 10000

# Sample sizes to test
sample_sizes <- c(10, 20, 50, 100)

# Significance level for the Pearson test
alpha <- 0.05

# Hampel identifier threshold
hampel_threshold <- 5

# --- Simulation Function ---
# This function runs one simulation for a given scenario.
# It generates data and tests for overdispersion using the three methods.
run_single_simulation <- function(n, dist_type, mu, variance = NULL) {
  
  # Generate data based on the distribution type
  if (dist_type == "poisson") {
    y <- rpois(n, lambda = mu)
  } else if (dist_type == "neg_binom") {
    if (is.null(variance)) {
      stop("Variance must be provided for Negative Binomial.")
    }
    # rnbinom in R is parameterized by 'mu' and 'size' (theta).
    # The relationship is: Var = mu + mu^2 / size
    # We solve for size: size = mu^2 / (Var - mu)
    if (variance <= mu) {
      # This would happen for underdispersion, rnbinom requires size > 0
      # We can approximate with Poisson if Var is very close to mu
      y <- rpois(n, lambda = mu)
    } else {
      size_param <- mu^2 / (variance - mu)
      y <- rnbinom(n, mu = mu, size = size_param)
    }
  } else {
    stop("Unknown distribution type.")
  }
  
  # Calculate sample statistics
  y_bar <- mean(y)
  y_var <- var(y)
  
  # Avoid division by zero if mean is 0
  if (y_bar == 0) {
    return(c(hampel_flag = FALSE, vmr_flag = FALSE, pearson_flag = FALSE))
  }
  
  # --- Apply the three rules ---
  
  # 1. Hampel Identifier: var - mean > threshold
  hampel_flag <- (y_var - y_bar) > hampel_threshold
  
  # 2. Variance-to-Mean Ratio (VMR): var / mean > 1
  vmr_flag <- (y_var / y_bar) > 1
  
  # 3. Pearson Chi-squared Test
  # Statistic S = sum((y_i - y_bar)^2 / y_bar) = (n-1) * VMR
  pearson_stat <- sum((y - y_bar)^2) / y_bar
  # Critical value from Chi-squared distribution with n-1 degrees of freedom
  crit_value <- qchisq(1 - alpha, df = n - 1)
  pearson_flag <- pearson_stat > crit_value
  
  # Return logical flags
  return(c(hampel_flag = hampel_flag, vmr_flag = vmr_flag, pearson_flag = pearson_flag))
}


# --- Main Simulation Loop ---

# Define the scenarios to run
scenarios <- list(
  list(name = "Strong OD (NB, mu=10, var=16)", dist = "neg_binom", mu = 10, var = 16),
  list(name = "Equidispersed (Poisson, mu=10)", dist = "poisson", mu = 10, var = 10),
  list(name = "Slight OD (NB, mu=1000, var=1006)", dist = "neg_binom", mu = 1000, var = 1006)
)

# Create a data frame to store results
results <- data.frame()

# Loop through each scenario and sample size
for (scenario in scenarios) {
  for (n in sample_sizes) {
    # Replicate the simulation N_sim times
    sim_outputs <- replicate(N_sim, run_single_simulation(n, scenario$dist, scenario$mu, scenario$var))
    
    # Calculate the proportion of times each method flagged overdispersion
    flag_proportions <- rowMeans(sim_outputs)
    
    # Store the results
    temp_result <- data.frame(
      Scenario = scenario$name,
      SampleSize_n = n,
      Hampel_Rate = flag_proportions["hampel_flag"],
      VMR_Rate = flag_proportions["vmr_flag"],
      Pearson_Rate = flag_proportions["pearson_flag"]
    )
    results <- rbind(results, temp_result)
  }
}

# --- Display Results ---
print(results, row.names = FALSE)

Simulation Results and Analysis

Running the code above produces the following results.

                         Scenario SampleSize_n Hampel_Rate VMR_Rate Pearson_Rate
 Strong OD (NB, mu=10, var=16)           10      0.4674   0.7818       0.2801
 Strong OD (NB, mu=10, var=16)           20      0.5408   0.8988       0.5186
 Strong OD (NB, mu=10, var=16)           50      0.6033   0.9859       0.8497
 Strong OD (NB, mu=10, var=16)          100      0.6382   0.9996       0.9877
  Equidispersed (Poisson, mu=10)           10      0.1388   0.4770       0.0465
  Equidispersed (Poisson, mu=10)           20      0.0385   0.4789       0.0505
  Equidispersed (Poisson, mu=10)           50      0.0016   0.4851       0.0494
  Equidispersed (Poisson, mu=10)          100      0.0000   0.4905       0.0520
Slight OD (NB, mu=1000, var=1006)           10      0.4452   0.4795       0.0491
Slight OD (NB, mu=1000, var=1006)           20      0.4571   0.4820       0.0483
Slight OD (NB, mu=1000, var=1006)           50      0.4727   0.4897       0.0543
Slight OD (NB, mu=1000, var=1006)          100      0.4850   0.4912       0.0528

Interpretation of Results

The simulation confirms the expectation perfectly:

  1. Scenario 1: Strong Overdispersion (mu=10, var=16)
    • Pearson Test: Power correctly increases with sample size n (from ~0.28 to ~0.99), effectively detecting the strong overdispersion in larger samples.
    • VMR > 1: This simple rule is highly sensitive, flagging overdispersion most of the time, as the true VMR is 1.6. Its detection rate naturally approaches 100% as n increases.
    • Hampel Identifier: With a threshold of 5, it is moderately sensitive. Its detection power increases with n but is significantly lower than the VMR rule.
  2. Scenario 2: Equidispersion (Poisson, mu=10)
    • Pearson Test: This is the key test. It correctly maintains a Type I error rate close to the target of 0.05 across all sample sizes.
    • VMR > 1: This is a descriptive flag, not a formal test. Since the sample VMR will randomly fall above or below the true value of 1, it flags overdispersion about 48-49% of the time, regardless of n.
    • Hampel Identifier: Its flagging rate decreases significantly as n grows. This is because for a true Poisson distribution, the sample variance converges to the mean. With larger n, it becomes increasingly rare for the sample variance to randomly exceed the mean by more than 5.
  3. Scenario 3: Slight Overdispersion at High Mean (mu=1000, var=1006)
    • Pearson Test: The true overdispersion is very slight (VMR ≈ 1.006). The Pearson test has very low power to detect this and correctly reports a significance rate near the null level of 0.05.
    • VMR > 1: Similar to the Poisson case, it flags about 48-49% of the time, as the true VMR is only slightly above 1.
    • Hampel Identifier: This highlights the rule’s scale-dependency. Because the variance (1006) and mean (1000) are large, the difference var - mean is centered around 6. A fixed threshold of 5 is easily exceeded by sampling variability. Consequently, the Hampel rule flags overdispersion about 45-48% of the time, essentially performing no better than a coin flip and dramatically overstating the severity of the dispersion.

Conclusion: The Pearson test is a formal, reliable statistical test that controls Type I error and gains power with sample size. The VMR rule is a simple descriptive measure. The Hampel identifier is a non-robust, scale-dependent heuristic that can be misleading, especially with high-mean data.

Post Checks

  • Parameterization check:
    • The Negative Binomial parameterization is correct: \(r = \mu^2 / (v - \mu)\) with \(p = r/(r+\mu)\) gives mean \(\mu\) and variance \(v = \mu + \mu^2/r\).
    • Scenario A: \(\mu=1000, v=1006 \Rightarrow r \approx 1{,}000{,}000/6 \approx 166{,}667\) (very large), so NB ≈ Poisson.
    • Scenario B: \(\mu=10, v=16 \Rightarrow r \approx 100/6 \approx 16.667\), so clear overdispersion vs Poisson.
  • Scenario A (VMR ≈ 1.006):
    • Poisson(true) and NB(target var) are essentially indistinguishable: mean VMR ≈ 1.006–1.009, Pearson p<0.05 around 5%. That’s exactly what we’d expect given the overdispersion is tiny.
    • Prop(VMR > 1) ≈ 0.47 for both, close to 0.5. Using “VMR > 1” as a flag is not a test; it will hover around 0.5 when the true VMR is 1.
  • Scenario B (VMR = 1.6):
    • Poisson(true) behaves well: mean VMR ≈ 1, false-positive rate ≈ 4–5%.
    • NB(target var) shows clear overdispersion: mean VMR ≈ 1.6, and the Pearson test flags about 46% at n=20. That’s reasonable power for this effect size at that n.
    • Prop(VMR > 1) ≈ 0.89 reflects that most samples show overdispersion, though not all, due to sampling variability.

Why the Pearson test results make sense

  • The Pearson statistic is \(S = \sum_i (y_i - \bar y)^2 / \bar y\), with test using \(S \sim \chi^2_{n-1}\) under the Poisson null and dispersion estimate \(\hat\phi = S/(n-1)\).
  • Under an alternative with variance \(\phi \mu\) (i.e., VMR = \(\phi\)), a good approximation is \(S \approx \phi \cdot \chi^2_{n-1}\). Hence power at level 0.05 is roughly:
    • \(\text{Power} \approx 1 - F_{\chi^2_{n-1}}\big(\chi^2_{0.95, n-1} / \phi\big)\).
  • Plugging \(\phi=1.6\):
    • n=20 (df=19): power ≈ 0.46–0.48, matching ~0.459.
    • n=30 (df=29): power should rise to roughly 0.60–0.70.
    • n=100 (df=99): power should be high, roughly 0.85–0.95.