statsmodels

# Statsmodels - Statistical Modeling & Inference

Statsmodels is the bridge between Python and the rigor of R-style statistical analysis. It allows users to estimate models using formulas (via patsy), perform extensive diagnostic tests, and produce detailed summary tables that are the standard in academic publishing.

## When to Use

- Estimating Linear Regression models with detailed diagnostics (OLS, WLS).
- Generalized Linear Models (GLM): Logistic, Poisson, Gamma regression.
- Time Series Analysis (ARIMA, SARIMAX, VAR, State Space models).
- Statistical hypothesis testing (t-tests, ANOVA, normality, heteroscedasticity).
- Survival analysis (Kaplan-Meier, Cox Proportional Hazards).
- Estimating treatment effects and causal inference.
- Non-parametric statistics (Kernel Density Estimation).

## Reference Documentation

**Official docs**: https://www.statsmodels.org/stable/  
**Formula API**: https://www.statsmodels.org/stable/example_formulas.html  
**Search patterns**: `sm.OLS`, `smf.ols`, `sm.tsa`, `results.summary()`, `statsmodels.api`

## Core Principles

### Statsmodels vs. scikit-learn

| Feature | scikit-learn | Statsmodels |
|---------|--------------|-------------|
| Goal | Prediction (Accuracy) | Inference (Explanation/p-values) |
| Interface | fit / predict | fit / summary |
| Pre-processing | Pipeline objects | Formulas (Patsy) or design matrices |
| Diagnostics | Cross-validation | Residue analysis, p-values, CI |

### The Two APIs

- **API (statsmodels.api)**: Requires explicit addition of a constant (intercept) and uses NumPy-like arrays.
- **Formula API (statsmodels.formula.api)**: Uses R-style formulas (`y ~ x1 + x2`) and works directly with Pandas DataFrames. (Recommended for most users).

## Quick Reference

### Installation

```bash
pip install statsmodels patsy
```

### Standard Imports

```python
import statsmodels.api as sm
import statsmodels.formula.api as smf
import pandas as pd
import numpy as np
```

### Basic Pattern - Ordinary Least Squares (OLS)

```python
import statsmodels.formula.api as smf

# 1. Define model with R-style formula
# 'y ~ x1 + x2' means: y = beta0 + beta1*x1 + beta2*x2
model = smf.ols('tip ~ total_bill + size', data=df_tips)

# 2. Fit the model
results = model.fit()

# 3. Print the comprehensive results table
print(results.summary())

# 4. Access specific values
p_values = results.pvalues
params = results.params # beta coefficients
```

## Critical Rules

### ✅ DO

- **Check Residuals** - Always plot and test residuals (`results.resid`) for normality and homoscedasticity.
- **Add a Constant** - If using the `sm.api` (not formula), remember `X = sm.add_constant(X)` or your model will pass through the origin (beta0 = 0).
- **Use Categorical Variables** - Use the `C()` operator in formulas (e.g., `y ~ C(region)`) to automatically create dummy variables.
- **Specify Covariance Type** - Use `cov_type='HC3'` or `'cluster'` if you suspect non-constant variance (heteroscedasticity).
- **Interpret R-squared carefully** - High R-squared doesn't imply a good model if the residuals are patterned.
- **Check for Multicollinearity** - Use VIF (Variance Inflation Factor) to ensure predictors aren't highly correlated.

### ❌ DON'T

- **Assume Prediction is Inference** - Just because a model has a high R-squared doesn't mean the coefficients represent real-world causal effects.
- **Ignore the Intercept** - Most physical and social processes require a constant term.
- **Overfit with too many predictors** - Use AIC/BIC metrics to penalize complex models.
- **Extrapolate beyond the range** - Statistical models are only valid within the domain of the training data.

## Anti-Patterns (NEVER)

```python
import statsmodels.api as sm

# ❌ BAD: Forgetting the intercept in Array API
model = sm.OLS(y, X) # This forces the line through (0,0)
results = model.fit()

# ✅ GOOD: Add the constant
X = sm.add_constant(X)
model = sm.OLS(y, X).fit()

# ❌ BAD: Ignoring categorical data strings
# smf.ols('price ~ color', data=df) # Might fail if 'color' isn't numeric

# ✅ GOOD: Explicitly tell Statsmodels it's categorical
smf.ols('price ~ C(color)', data=df).fit()

# ❌ BAD: Only looking at R-squared
# print(results.rsquared) # Tells only part of the story

# ✅ GOOD: Inspecting the whole summary
print(results.summary()) # Check p-values, F-stat, Jarque-Bera
```

## Regression Analysis

### Linear Models (OLS, WLS)

```python
# Multiple Linear Regression with interactions
# 'x1 * x2' includes x1, x2, and the interaction term x1:x2
model = smf.ols('y ~ x1 * x2 + np.log(x3)', data=df).fit()

# Weighted Least Squares (for heteroscedastic data)
wls_model = sm.WLS(y, X, weights=1.0/variance_estimates).fit()
```

### Generalized Linear Models (GLM)

```python
# Logistic Regression (Binary outcome)
logit_model = smf.logit('admit ~ gre + gpa + C(rank)', data=df).fit()

# Poisson Regression (Count data)
poisson_model = smf.poisson('num_awards ~ math + C(prog)', data=df).fit()

# Negative Binomial (For over-dispersed counts)
nb_model = smf.glm('y ~ x1', data=df, family=sm.families.NegativeBinomial()).fit()
```

## Time Series Analysis (tsa)

### Stationarity and Modeling

```python
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX

# 1. Test for Stationarity (Augmented Dickey-Fuller)
adf_result = adfuller(df['sales'])
print(f"ADF P-value: {adf_result[1]}") # p < 0.05 means stationary

# 2. SARIMAX (Seasonal ARIMA with eXogenous variables)
model = SARIMAX(df['sales'], 
                order=(1, 1, 1), 
                seasonal_order=(1, 1, 0, 12),
                exog=df['advertising'])
results = model.fit()

# 3. Forecasting
forecast = results.get_forecast(steps=12)
conf_int = forecast.conf_int()
```

## ANOVA and Hypothesis Testing

```python
from statsmodels.stats.anova import anova_lm

# Perform ANOVA on OLS model
model = smf.ols('yield ~ C(fertilizer) + C(soil)', data=df).fit()
anova_table = anova_lm(model, typ=2)

# Post-hoc tests (Tukey's HSD)
from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey = pairwise_tukeyhsd(endog=df['yield'], groups=df['fertilizer'], alpha=0.05)
print(tukey)
```

## Model Diagnostics

### Residual Analysis

```python
import statsmodels.stats.api as sms

# 1. Test for Heteroscedasticity (Breusch-Pagan)
name = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sms.het_breuschpagan(results.resid, results.model.exog)
print(dict(zip(name, test)))

# 2. Test for Normality (Omnibus)
# Included in results.summary() by default. 

# 3. Check for Outliers (Influence)
influence = results.get_influence()
cooks_d = influence.cook_distance[0]
```

## Practical Workflows

### 1. Robust Scientific Reporting Pipeline

```python
def analyze_experiment(df):
    """Rigorous analysis of an experimental dataset."""
    # 1. Fit model
    model = smf.ols('outcome ~ treatment + age + gender', data=df).fit()
    
    # 2. Diagnostic Plots (requires matplotlib)
    import matplotlib.pyplot as plt
    sm.graphics.plot_regress_exog(model, 'treatment')
    
    # 3. Check for Multicollinearity
    from statsmodels.stats.outliers_influence import variance_inflation_factor
    # (calculation for VIF...)
    
    return model.summary()
```

### 2. Market Mix Modeling (Attribution)

```python
def estimate_attribution(df):
    # Log-log model to calculate elasticities
    # Coefficients will be % change in sales for 1% change in spend
    model = smf.ols('np.log(sales) ~ np.log(tv_spend) + np.log(digital_spend)', data=df).fit()
    return model.params
```

### 3. Survival Analysis

```python
from statsmodels.duration.hazard_regression import PHReg

# Cox Proportional Hazards Model
model = PHReg.from_formula('time ~ age + C(treatment)', data=df, status=df['event'])
results = model.fit()
print(results.summary())
```

## Performance Optimization

### Using numba for Likelihoods

While Statsmodels is primarily written in Python and Cython, some of the newer time series modules utilize optimized numerical backends for faster fitting of state-space models.

### Formulas vs Design Matrices

For very large datasets (1M+ rows), creating the design matrix with patsy can be memory-intensive. In these cases, construct your X matrix manually and use `sm.OLS(y, X)`.

## Common Pitfalls and Solutions

### Singular Matrix Error

"LinAlgError: Singular matrix" means your predictors are perfectly correlated (e.g., including both temp_celsius and temp_fahrenheit).

```python
# ✅ Solution: Remove redundant columns
df = df.drop('redundant_col', axis=1)
```

### Categorical Leakage (The Dummy Variable Trap)

Including intercept and dummy variables for ALL categories creates perfect multicollinearity.

```python
# ❌ Problem: y ~ dummy_cat1 + dummy_cat2 + dummy_cat3 + intercept 
# ✅ Solution: Statsmodels/Patsy automatically drops one category (ref category)
# to avoid this. Don't try to force all dummies in!
```

### Non-Stationary Time Series

Predicting a non-stationary series leads to "spurious regression".

```python
# ✅ Solution: Difference your data first
df['diff_y'] = df['y'].diff()
# Or use ARIMA with integrated term (d=1)
```

Statsmodels is the gold standard for statistical validity in the Python ecosystem. It moves beyond black-box predictions to provide the transparency and mathematical rigor required for high-stakes scientific and economic decision-making.