Slide 1: Introduction to Dispersion Parameter
The dispersion parameter is a crucial concept in statistical modeling, particularly in generalized linear models (GLMs). It measures the variability of data points around the predicted values. Understanding this parameter helps in assessing model fit and making accurate predictions.
import numpy as np
import matplotlib.pyplot as plt
# Generate sample data
np.random.seed(42)
x = np.linspace(0, 10, 100)
y_low_dispersion = 2 * x + np.random.normal(0, 1, 100)
y_high_dispersion = 2 * x + np.random.normal(0, 5, 100)
# Plot data with different dispersion
plt.figure(figsize=(10, 6))
plt.scatter(x, y_low_dispersion, alpha=0.5, label='Low Dispersion')
plt.scatter(x, y_high_dispersion, alpha=0.5, label='High Dispersion')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.title('Data with Different Dispersion')
plt.show()Slide 2: Mathematical Definition
The dispersion parameter, often denoted as φ (phi), is defined as the variance of the response variable divided by its mean. In mathematical terms:
φ = Var(Y) / E(Y)
For some distributions, like the Poisson, the dispersion parameter is fixed at 1. For others, like the negative binomial, it's a parameter to be estimated.
import numpy as np
from scipy import stats
# Generate data from Poisson and Negative Binomial distributions
np.random.seed(42)
poisson_data = stats.poisson.rvs(mu=5, size=1000)
negbin_data = stats.nbinom.rvs(n=5, p=0.5, size=1000)
# Calculate dispersion for Poisson (should be close to 1)
poisson_dispersion = np.var(poisson_data) / np.mean(poisson_data)
# Calculate dispersion for Negative Binomial
negbin_dispersion = np.var(negbin_data) / np.mean(negbin_data)
print(f"Poisson dispersion: {poisson_dispersion:.2f}")
print(f"Negative Binomial dispersion: {negbin_dispersion:.2f}")Slide 3: Importance in GLMs
In Generalized Linear Models, the dispersion parameter plays a crucial role in determining the relationship between the mean and variance of the response variable. It affects the standard errors of coefficient estimates and confidence intervals.
import statsmodels.api as sm
import numpy as np
# Generate sample data
np.random.seed(42)
X = np.random.rand(100, 1)
y = np.random.poisson(np.exp(1 + 2 * X).ravel())
# Fit Poisson GLM
poisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson())
poisson_results = poisson_model.fit()
# Print summary
print(poisson_results.summary())Slide 4: Overdispersion
Overdispersion occurs when the observed variance is greater than the theoretical variance predicted by the model. It's common in count data and can lead to underestimated standard errors if not addressed.
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# Generate overdispersed data
np.random.seed(42)
lambda_param = 5
size = 1000
overdispersed_data = np.random.negative_binomial(n=3, p=3/(3+lambda_param), size=size)
# Plot histogram
plt.figure(figsize=(10, 6))
plt.hist(overdispersed_data, bins=range(max(overdispersed_data)+2), alpha=0.7)
plt.xlabel('Count')
plt.ylabel('Frequency')
plt.title('Histogram of Overdispersed Count Data')
plt.show()
# Calculate dispersion
dispersion = np.var(overdispersed_data) / np.mean(overdispersed_data)
print(f"Dispersion: {dispersion:.2f}")Slide 5: Underdispersion
Underdispersion is less common but can occur when the observed variance is less than the theoretical variance. It may indicate that the model is not capturing all the structure in the data.
import numpy as np
import matplotlib.pyplot as plt
# Generate underdispersed data
np.random.seed(42)
n = 1000
p = 0.3
underdispersed_data = np.random.binomial(n, p, size=1000)
# Plot histogram
plt.figure(figsize=(10, 6))
plt.hist(underdispersed_data, bins=30, alpha=0.7)
plt.xlabel('Count')
plt.ylabel('Frequency')
plt.title('Histogram of Underdispersed Count Data')
plt.show()
# Calculate dispersion
dispersion = np.var(underdispersed_data) / (n * p * (1 - p))
print(f"Dispersion: {dispersion:.2f}")Slide 6: Estimating Dispersion Parameter
There are several methods to estimate the dispersion parameter. One common approach is the Pearson chi-square statistic divided by the degrees of freedom.
import numpy as np
import statsmodels.api as sm
# Generate sample data
np.random.seed(42)
X = np.random.rand(100, 1)
y = np.random.negative_binomial(n=5, p=0.3, size=100)
# Fit Negative Binomial GLM
nb_model = sm.GLM(y, sm.add_constant(X), family=sm.families.NegativeBinomial())
nb_results = nb_model.fit()
# Calculate dispersion using Pearson chi-square
pearson_chi2 = nb_results.pearson_chi2
degrees_of_freedom = nb_results.df_resid
estimated_dispersion = pearson_chi2 / degrees_of_freedom
print(f"Estimated dispersion parameter: {estimated_dispersion:.2f}")Slide 7: Impact on Standard Errors
The dispersion parameter affects the standard errors of coefficient estimates. When overdispersion is present and not accounted for, standard errors may be underestimated, leading to overconfident inferences.
import numpy as np
import statsmodels.api as sm
# Generate overdispersed data
np.random.seed(42)
X = np.random.rand(1000, 1)
y = np.random.negative_binomial(n=5, p=0.3, size=1000)
# Fit Poisson GLM (ignoring overdispersion)
poisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson())
poisson_results = poisson_model.fit()
# Fit Negative Binomial GLM (accounting for overdispersion)
nb_model = sm.GLM(y, sm.add_constant(X), family=sm.families.NegativeBinomial())
nb_results = nb_model.fit()
print("Poisson Model (ignoring overdispersion):")
print(poisson_results.summary().tables[1])
print("\nNegative Binomial Model (accounting for overdispersion):")
print(nb_results.summary().tables[1])Slide 8: Dispersion in Different Distributions
Different probability distributions have different relationships between their mean and variance, which affects their dispersion characteristics.
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# Generate data from different distributions
np.random.seed(42)
size = 1000
poisson_data = stats.poisson.rvs(mu=5, size=size)
normal_data = stats.norm.rvs(loc=5, scale=2, size=size)
binomial_data = stats.binom.rvs(n=10, p=0.5, size=size)
# Calculate dispersion
poisson_dispersion = np.var(poisson_data) / np.mean(poisson_data)
normal_dispersion = np.var(normal_data) / np.mean(normal_data)
binomial_dispersion = np.var(binomial_data) / np.mean(binomial_data)
# Plot histograms
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))
ax1.hist(poisson_data, bins=20, alpha=0.7)
ax1.set_title(f'Poisson (φ={poisson_dispersion:.2f})')
ax2.hist(normal_data, bins=20, alpha=0.7)
ax2.set_title(f'Normal (φ={normal_dispersion:.2f})')
ax3.hist(binomial_data, bins=20, alpha=0.7)
ax3.set_title(f'Binomial (φ={binomial_dispersion:.2f})')
plt.tight_layout()
plt.show()Slide 9: Quasi-likelihood Approach
When the true underlying distribution is unknown, a quasi-likelihood approach can be used. This method allows for flexible modeling of the mean-variance relationship.
import numpy as np
import statsmodels.api as sm
# Generate overdispersed count data
np.random.seed(42)
X = np.random.rand(1000, 1)
y = np.random.negative_binomial(n=5, p=0.3, size=1000)
# Fit quasi-Poisson model
quasipoisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson())
quasipoisson_results = quasipoisson_model.fit(scale='X2')
print(quasipoisson_results.summary().tables[1])
print(f"\nEstimated dispersion parameter: {quasipoisson_results.scale:.2f}")Slide 10: Detecting Overdispersion
Several methods can be used to detect overdispersion, including visual inspection of residuals and formal statistical tests.
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy import stats
# Generate overdispersed data
np.random.seed(42)
X = np.random.rand(1000, 1)
y = np.random.negative_binomial(n=5, p=0.3, size=1000)
# Fit Poisson GLM
poisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson())
poisson_results = poisson_model.fit()
# Calculate residuals
residuals = poisson_results.resid_pearson
# Plot residuals
plt.figure(figsize=(10, 6))
plt.scatter(poisson_results.mu, residuals, alpha=0.5)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Fitted values')
plt.ylabel('Pearson residuals')
plt.title('Residual Plot for Detecting Overdispersion')
plt.show()
# Perform Cameron & Trivedi's test for overdispersion
y_pred = poisson_results.predict(sm.add_constant(X))
aux_ols = sm.OLS(((y - y_pred)**2 - y) / y_pred, y_pred).fit()
print(f"Cameron & Trivedi's test p-value: {aux_ols.pvalues[0]:.4f}")Slide 11: Handling Overdispersion
When overdispersion is detected, several approaches can be used to address it, including using a different distribution (e.g., negative binomial) or adjusting standard errors.
import numpy as np
import statsmodels.api as sm
# Generate overdispersed data
np.random.seed(42)
X = np.random.rand(1000, 1)
y = np.random.negative_binomial(n=5, p=0.3, size=1000)
# Fit Poisson GLM
poisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson())
poisson_results = poisson_model.fit()
# Fit Negative Binomial GLM
nb_model = sm.GLM(y, sm.add_constant(X), family=sm.families.NegativeBinomial())
nb_results = nb_model.fit()
# Fit Quasi-Poisson GLM
quasipoisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson())
quasipoisson_results = quasipoisson_model.fit(scale='X2')
print("Poisson Model:")
print(poisson_results.summary().tables[1])
print("\nNegative Binomial Model:")
print(nb_results.summary().tables[1])
print("\nQuasi-Poisson Model:")
print(quasipoisson_results.summary().tables[1])Slide 12: Real-life Example: Ecological Study
In an ecological study, researchers are investigating the number of bird species in different forest patches. The dispersion parameter helps account for variability in species counts beyond what's expected from area alone.
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
# Generate synthetic data
np.random.seed(42)
forest_area = np.random.uniform(1, 100, 100)
species_count = np.random.negative_binomial(n=5, p=5/(5+forest_area), size=100)
# Fit Poisson and Negative Binomial models
X = sm.add_constant(np.log(forest_area))
poisson_model = sm.GLM(species_count, X, family=sm.families.Poisson()).fit()
nb_model = sm.GLM(species_count, X, family=sm.families.NegativeBinomial()).fit()
# Plot data and fitted curves
plt.figure(figsize=(10, 6))
plt.scatter(forest_area, species_count, alpha=0.5)
plt.plot(forest_area, poisson_model.predict(X), 'r-', label='Poisson')
plt.plot(forest_area, nb_model.predict(X), 'g-', label='Negative Binomial')
plt.xlabel('Forest Area (ha)')
plt.ylabel('Number of Bird Species')
plt.legend()
plt.title('Bird Species Richness vs Forest Area')
plt.show()
print(f"Poisson AIC: {poisson_model.aic:.2f}")
print(f"Negative Binomial AIC: {nb_model.aic:.2f}")Slide 13: Real-life Example: Manufacturing Quality Control
In a manufacturing process, the number of defects per unit is being monitored. The dispersion parameter helps identify whether the defect occurrence is more variable than expected under a simple Poisson process.
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
# Generate synthetic data
np.random.seed(42)
production_run = np.arange(1, 101)
defect_count = np.random.negative_binomial(n=3, p=0.5, size=100)
# Fit Poisson and Negative Binomial models
X = sm.add_constant(production_run)
poisson_model = sm.GLM(defect_count, X, family=sm.families.Poisson()).fit()
nb_model = sm.GLM(defect_count, X, family=sm.families.NegativeBinomial()).fit()
# Plot data and fitted curves
plt.figure(figsize=(10, 6))
plt.scatter(production_run, defect_count, alpha=0.5)
plt.plot(production_run, poisson_model.predict(X), 'r-', label='Poisson')
plt.plot(production_run, nb_model.predict(X), 'g-', label='Negative Binomial')
plt.xlabel('Production Run')
plt.ylabel('Number of Defects')
plt.legend()
plt.title('Defect Count in Manufacturing Process')
plt.show()
# Compare model fits
print(f"Poisson AIC: {poisson_model.aic:.2f}")
print(f"Negative Binomial AIC: {nb_model.aic:.2f}")
print(f"Estimated dispersion parameter: {nb_model.scale:.2f}")Slide 14: Dispersion and Model Selection
The dispersion parameter plays a crucial role in model selection, especially when comparing models with different distributional assumptions. It affects measures like AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion).
import numpy as np
import statsmodels.api as sm
# Generate synthetic data
np.random.seed(42)
X = np.random.rand(1000, 2)
y = np.random.negative_binomial(n=5, p=0.3, size=1000)
# Fit different models
poisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson()).fit()
nb_model = sm.GLM(y, sm.add_constant(X), family=sm.families.NegativeBinomial()).fit()
quasipoisson_model = sm.GLM(y, sm.add_constant(X), family=sm.families.Poisson()).fit(scale='X2')
# Compare model fits
models = [poisson_model, nb_model, quasipoisson_model]
model_names = ['Poisson', 'Negative Binomial', 'Quasi-Poisson']
for name, model in zip(model_names, models):
print(f"{name}:")
print(f" AIC: {model.aic:.2f}")
print(f" BIC: {model.bic:.2f}")
print(f" Log-likelihood: {model.llf:.2f}")
print(f" Dispersion: {model.scale:.2f}\n")Slide 15: Additional Resources
For further exploration of the dispersion parameter and its applications in statistical modeling, consider the following resources:
- "Generalized Linear Models" by McCullagh and Nelder (1989) - A comprehensive text on GLMs, including detailed discussions on dispersion.
- "Extending the Linear Model with R" by Faraway (2016) - Provides practical examples of working with dispersion in R, with concepts applicable to Python.
- ArXiv paper: "On the Estimation of Dispersion Parameters in GLMs" by Kokonendji et al. (2020) ArXiv URL: https://arxiv.org/abs/2007.03521
- statsmodels documentation: https://www.statsmodels.org/stable/glm.html - Offers detailed information on implementing GLMs in Python, including handling dispersion.
These resources provide a mix of theoretical foundations and practical applications, allowing for a deeper understanding of the dispersion parameter and its role in statistical modeling.