Workflow 1: Ecological Modeling
Goal: Build an interpretable species distribution model from highly correlated bioclimatic variables.
Challenge: WorldClim’s 19 bioclimatic variables contain many temperature and precipitation metrics that are mathematically related (e.g., mean temperature, minimum temperature, maximum temperature). Using all variables leads to multicollinearity, unstable coefficients, and poor model interpretation.
Strategy: Two-stage pruning:
corrPrune()removes pairwise correlations > 0.7modelPrune()refines further using variance inflation factors (VIF)
This approach balances model fit with interpretability and numerical stability.
Data
data(bioclim_example)
# Data structure
dim(bioclim_example)
#> [1] 100 20
head(names(bioclim_example))
#> [1] "species_richness" "BIO1" "BIO2" "BIO3"
#> [5] "BIO4" "BIO5"
# Response variable
summary(bioclim_example$species_richness)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 5.0 96.0 118.5 120.2 144.2 206.0The dataset contains 100 sampling locations with 19 bioclimatic predictors and a species richness response. The bioclimatic variables are standard WorldClim metrics (BIO1-BIO19) measuring temperature and precipitation patterns.
Correlation-based pruning
We start by removing variables with pairwise correlations exceeding
0.7. The mode = "auto" setting uses the exact algorithm
(Bron-Kerbosch) to enumerate all maximal subsets, then selects the
largest one.
# Remove highly correlated predictors
bio_clean <- corrPrune(
data = bioclim_example[, -1], # Exclude response
threshold = 0.7,
mode = "auto"
)
# How much did we reduce?
cat(sprintf("Reduced from %d → %d variables\n",
ncol(bioclim_example) - 1,
ncol(bio_clean)))
#> Reduced from 19 → 12 variables
# Which variables were kept?
head(attr(bio_clean, "selected_vars"), 10)
#> [1] "BIO1" "BIO3" "BIO6" "BIO9" "BIO11" "BIO12" "BIO13" "BIO14" "BIO16"
#> [10] "BIO17"The pruning successfully reduced dimensionality while ensuring no remaining pair exceeds the 0.7 threshold. The selected variables span both temperature and precipitation domains, maintaining ecological interpretability.
Correlation distribution
The histogram below shows how corrPrune() reshapes the
correlation structure. Before pruning (red), many variable pairs exceed
the 0.7 threshold. After pruning (blue), all pairwise correlations fall
below the threshold.
cor_before <- cor(bioclim_example[, -1])
cor_after <- cor(bio_clean)
vals_before <- abs(cor_before[upper.tri(cor_before)])
vals_after <- abs(cor_after[upper.tri(cor_after)])
# Common x limits
xlim <- c(0, 1)
# Shared breaks for fair comparison
breaks <- seq(0, 1, length.out = 30)
# Before histogram
hist(vals_before,
breaks = breaks,
freq = FALSE,
main = "Distribution of Absolute Correlations",
xlab = "Absolute Correlation",
col = rgb(0.8, 0.2, 0.2, 0.4),
border = "white",
xlim = xlim)
# After histogram
hist(vals_after,
breaks = breaks,
freq = FALSE,
col = rgb(0.2, 0.5, 0.8, 0.4),
border = "white",
add = TRUE)
# Threshold line
abline(v = 0.7, col = "black", lty = 2, lwd = 2)
legend("topright",
legend = c("Before", "After", "Threshold"),
fill = c(rgb(0.8, 0.2, 0.2, 0.4),
rgb(0.2, 0.5, 0.8, 0.4),
NA),
border = c("white", "white", NA),
lty = c(NA, NA, 2),
lwd = c(NA, NA, 2),
col = c(NA, NA, "black"),
bty = "o",
bg = "white")
Fit models
We fit three models to compare the effect of pruning:
- Full model: All 19 bioclimatic variables (baseline)
- After corrPrune: Variables with correlations < 0.7
- After modelPrune: Further refined using VIF < 5
The VIF criterion complements correlation-based pruning by detecting multicollinearity involving more than two variables simultaneously.
# Model 1: Full model (19 variables)
model_full <- lm(species_richness ~ ., data = bioclim_example)
# Model 2: After corrPrune (correlation-based pruning)
bio_clean_full <- data.frame(
species_richness = bioclim_example$species_richness,
bio_clean
)
model_corrprune <- lm(species_richness ~ ., data = bio_clean_full)
# Model 3: Sequential VIF-based refinement
bio_final <- modelPrune(
formula = species_richness ~ .,
data = bio_clean_full,
limit = 5
)
model_final <- attr(bio_final, "final_model")Model comparison
The table below shows that pruning dramatically improves numerical stability (condition number κ) while maintaining model fit (adjusted R²). A lower κ indicates better-conditioned matrices and more stable coefficient estimates.
# Variable counts
n_full <- 19
n_corrprune <- length(attr(bio_clean, "selected_vars"))
n_final <- length(attr(bio_final, "selected_vars"))
# Compute condition numbers (measure of collinearity)
X_full <- model.matrix(model_full)[, -1]
X_corrprune <- model.matrix(model_corrprune)[, -1]
X_final <- model.matrix(model_final)[, -1]
kappa_full <- kappa(X_full, exact = TRUE)
kappa_corrprune <- kappa(X_corrprune, exact = TRUE)
kappa_final <- kappa(X_final, exact = TRUE)
# Summary table
comparison <- data.frame(
Step = c("Full", "corrPrune", "+ modelPrune"),
Predictors = c(n_full, n_corrprune, n_final),
Adj_R2 = c(
summary(model_full)$adj.r.squared,
summary(model_corrprune)$adj.r.squared,
summary(model_final)$adj.r.squared
),
Kappa = c(kappa_full, kappa_corrprune, kappa_final)
)
print(comparison)
#> Step Predictors Adj_R2 Kappa
#> 1 Full 19 0.9821862 914.8146
#> 2 corrPrune 12 0.8962848 538.1157
#> 3 + modelPrune 12 0.8962848 538.1157Key insights:
- The full model has substantial multicollinearity (high κ)
- Correlation-based pruning reduces κ while losing minimal fit
- VIF-based refinement can further improve stability with only marginal R² decrease
Visualization
The plot below illustrates the pruning tradeoff: as we remove variables, condition number (κ) drops dramatically while adjusted R² remains high. This demonstrates that many of the original 19 variables were redundant for prediction.
# Extract data
n_vars <- comparison$Predictors
adj_r2 <- comparison$Adj_R2
kappa <- comparison$Kappa
# Left y-axis: Adjusted R²
par(mar = c(5, 4, 4, 4)) # extra space on the right for second axis
plot(
n_vars, adj_r2,
type = "b",
pch = 19,
cex = 1.5,
col = rgb(0.2, 0.5, 0.8, 1),
lwd = 2,
xlab = "Number of Predictors",
ylab = "Adjusted R²",
ylim = c(0, 1),
main = "Pruning Reduces Collinearity While Preserving Fit"
)
# Right y-axis: log10(κ)
log_kappa <- log10(kappa)
# Nice ylim for log10(κ)
ylim_right <- range(log_kappa, finite = TRUE)
ylim_right <- ylim_right * c(0.9, 1.1)
par(new = TRUE)
plot(
n_vars, log_kappa,
type = "b",
pch = 17,
cex = 1.5,
col = rgb(0.8, 0.2, 0.2, 1),
lwd = 2,
xaxt = "n",
yaxt = "n",
xlab = "",
ylab = "",
ylim = ylim_right
)
# Right-hand axis ticks using pretty() on log scale
log_ticks <- pretty(log_kappa)
kappa_labels <- round(10^log_ticks)
axis(4, at = log_ticks, labels = kappa_labels)
mtext("Condition Number (κ)", side = 4, line = 3)
# Legend centered at top
legend(
"top",
inset = 0.02,
legend = c("Adjusted R² (higher better)", "κ (lower better)"),
col = c(
rgb(0.2, 0.5, 0.8, 1),
rgb(0.8, 0.2, 0.2, 1)
),
pch = c(19, 17),
lwd = 2,
horiz = TRUE,
bty = "o",
bg = "white",
x.intersp = 0.8
)
Coefficient stability
Multicollinearity inflates coefficient variance, making estimates unstable. The plot below compares coefficients between the full model (19 variables) and the final pruned model. Notice how:
- Variables dropped in the pruned model (shown as red-only bars) had unstable estimates
- Variables retained in both models (overlapping bars) show more consistent magnitudes
- The pruned model yields clearer, more interpretable effect sizes
# Extract coefficients (excluding intercept)
coef_full <- coef(model_full)[-1]
coef_final <- coef(model_final)[-1]
# Use *all* variables (not just common ones)
all_vars <- union(names(coef_full), names(coef_final))
# Align coefficients to the same full variable list
vals_full <- coef_full[all_vars]
vals_pruned <- coef_final[all_vars]
# Replace missing values (variables dropped in a model) with 0
vals_full[is.na(vals_full)] <- 0
vals_pruned[is.na(vals_pruned)] <- 0
# Colours
col_full <- rgb(0.8, 0.2, 0.2, 0.5)
col_pruned <- rgb(0.2, 0.5, 0.8, 0.5)
x <- seq_along(all_vars)
# Symmetric Y range
ylim <- range(c(vals_full, vals_pruned)) * 1.15
# Empty plot first
plot(
x, vals_full,
type = "n",
xaxt = "n",
xlab = "",
ylab = "Coefficient",
main = "Coefficient Comparison (Full vs modelPrune)",
ylim = ylim
)
axis(1, at = x, labels = all_vars, las = 2, cex.axis = 0.7)
# Full model bars
rect(
x - 0.4, 0,
x + 0.4, vals_full,
col = col_full, border = NA
)
# Pruned model bars
rect(
x - 0.4, 0,
x + 0.4, vals_pruned,
col = col_pruned, border = NA
)
legend(
"topright",
legend = c(
sprintf("Full (%d vars)", n_full),
sprintf("Final (%d vars)", n_final)
),
fill = c(col_full, col_pruned),
border = "white",
bty = "o",
bg = "white"
)
