Reducing Predictor Collinearity with VIF in Python

Highly correlated environmental predictors destabilise species distribution models: coefficients swing wildly, response curves become uninterpretable, and the fitted model transfers poorly to new landscapes even when its training AUC looks fine. Pruning that redundancy is the last quality gate of Environmental Predictor Stacking, the stage of the Species Distribution Modeling with MaxEnt workflow that assembles the covariate cube. Bioclimatic layers are the worst offenders — annual mean temperature, mean temperature of the warmest quarter, and minimum temperature of the coldest month all encode the same broad thermal gradient. This page prunes a predictor set two complementary ways: a pairwise absolute-correlation filter and an iterative Variance Inflation Factor (VIF) filter that catches multi-way collinearity a pairwise check misses.

When to Use VIF vs Correlation Pruning

The two methods detect different structures. Pairwise |r| catches two-variable redundancy; VIF catches a predictor that is collinear with a combination of others even when no single pair is strongly correlated. Most robust pipelines apply both.

VIF for predictor $j$ is defined from the coefficient of determination $R_j^2$ of regressing that predictor on all the others:

$${\text{VIF}}_j=\frac{1}{1-R_j^2}$$

A VIF of 1 means the predictor is orthogonal to the rest; a VIF of 10 means 90% of its variance is explained by the other predictors ($R_j^2=0.9$), so it contributes almost no independent information. The retained set then feeds MaxEnt model training & tuning, where fewer, independent predictors make regularisation and feature-class selection far more stable.

Minimal Reproducible Example

Sample the predictor stack at a large set of random points (or the background points), assemble a tidy DataFrame, then run the pairwise filter followed by the iterative VIF filter. Sampling rather than using every pixel keeps the regressions fast and avoids memory blow-up on large stacks.

import numpy as np
import pandas as pd
import rasterio
from statsmodels.stats.outliers_influence import variance_inflation_factor

# 1. Sample the predictor stack at random valid cells.
with rasterio.open("predictor_stack.tif") as src:
    names = list(src.descriptions) or [f"b{i+1}" for i in range(src.count)]
    data = src.read(masked=True)                 # (bands, rows, cols)

flat = data.reshape(data.shape[0], -1)
valid = ~np.ma.getmaskarray(flat).any(axis=0)    # cells valid in every band
rng = np.random.default_rng(0)
sample_cols = rng.choice(np.flatnonzero(valid), size=5000, replace=False)
df = pd.DataFrame(flat[:, sample_cols].filled(np.nan).T, columns=names).dropna()


def prune_correlation(frame, threshold=0.7):
    """Drop one of each pair of predictors with |r| >= threshold."""
    corr = frame.corr().abs()
    upper = corr.where(np.triu(np.ones(corr.shape, dtype=bool), k=1))
    drop = [c for c in upper.columns if (upper[c] >= threshold).any()]
    return frame.drop(columns=drop), drop


def prune_vif(frame, threshold=10.0):
    """Iteratively drop the highest-VIF predictor until all are below threshold."""
    cols = list(frame.columns)
    while len(cols) > 1:
        X = frame[cols].values
        vifs = [variance_inflation_factor(X, i) for i in range(len(cols))]
        worst = int(np.argmax(vifs))
        if vifs[worst] < threshold:
            break
        cols.pop(worst)
    return frame[cols], cols


df_r, dropped_r = prune_correlation(df, threshold=0.7)
df_final, kept = prune_vif(df_r, threshold=10.0)

print(f"correlation dropped: {dropped_r}")
print(f"retained predictors: {kept}")

Run the correlation filter first: it is cheap and removes the obvious duplicates, which shrinks the VIF regressions and makes the iterative step converge faster. The order matters only for speed — VIF alone would eventually catch the same pairs.

Parameter Reference

Expected Output and Verification

The output is a reduced predictor list in which every retained variable is below both thresholds. Verify the final set directly — recompute VIFs and the correlation matrix and assert the invariants, because an iterative loop that breaks early is a common silent bug.

import numpy as np
from statsmodels.stats.outliers_influence import variance_inflation_factor

X = df_final.values
final_vifs = {c: variance_inflation_factor(X, i) for i, c in enumerate(df_final.columns)}

# Every retained predictor must be below the VIF threshold.
assert all(v < 10.0 for v in final_vifs.values()), f"VIF still high: {final_vifs}"

# No retained pair may exceed the correlation threshold.
corr = df_final.corr().abs().values
np.fill_diagonal(corr, 0.0)
assert corr.max() < 0.7, f"Correlated pair remains: max |r|={corr.max():.2f}"

print("Retained predictors and VIFs:")
for c, v in sorted(final_vifs.items(), key=lambda kv: kv[1]):
    print(f"  {c:<12} VIF={v:.2f}")

Expect the retained set to be roughly 4–8 predictors from an initial dozen bioclimatic layers, all with VIF comfortably below the threshold. If the loop leaves only one or two, your threshold is too strict for a genuinely multicollinear climate space — relax it and lean on the ecological-primacy rule instead of dropping mechanistically important variables. The pruned stack is what every later stage samples, including the honest evaluation in Model Validation & AUC Metrics.

Common Pitfalls

• Dropping by VIF alone without ecological judgment. VIF ranks by statistical redundancy, not biological relevance. When two collinear predictors both exceed the threshold, override the automatic “drop highest VIF” rule to keep the one with the clearer mechanistic link — for example, retain minimum temperature of the coldest month over annual mean temperature for a frost-limited species.
• Computing VIF on unsampled full rasters. Feeding millions of pixels into repeated OLS regressions is slow and memory-hungry, and it does not improve the estimate. Sample a few thousand valid cells instead.
• Leaving NaN or masked cells in the matrix. variance_inflation_factor returns inf or garbage if the design matrix contains NaNs. Mask to cells valid across every band and dropna() before pruning, as in the example.
• Correlation and VIF disagreeing. A pair can pass |r| < 0.7 yet still inflate VIF because three or more predictors form a near-linear combination. That is precisely why you run VIF after the pairwise filter rather than trusting correlation alone.

Frequently Asked Questions

Related

• Environmental Predictor Stacking — parent workflow that aligns and stacks the predictors pruned here.
• Stacking Climate Layers for SDM in Python — the sibling guide that builds the bioclimatic stack this filter thins.
• MaxEnt Model Training & Tuning — where the reduced predictor set is fitted and regularised.
• Species Distribution Modeling with MaxEnt — the full modelling workflow this guide belongs to.

Up one level: Environmental Predictor Stacking · Species Distribution Modeling with MaxEnt