Complete Workflows: Real-World Examples"

In corrselect: Correlation-Based and Model-Based Predictor Pruning

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  dev = "svglite",
  fig.ext = "svg",
  cache = TRUE  # Cache expensive computations to speed up builds
)
library(corrselect)

Overview

This vignette demonstrates integration of corrselect into complete modeling workflows, from raw predictors to final models under correlation constraints. Four applied settings illustrate interface selection and workflow composition. See vignette("quickstart") for interface descriptions and vignette("theory") for mathematical foundations.

Functions used: corrPrune(), modelPrune(), see vignette("quickstart") for full signatures and parameter details.

Each workflow showcases different aspects of the package:

1. Ecological Modeling: Two-stage pruning (correlation + VIF) for environmental predictors

2. Survey Data: Protecting key variables while reducing redundancy in questionnaires

3. High-Dimensional Data: Greedy algorithms for gene expression (p >> n scenarios)

4. Mixed Models: VIF-based pruning of fixed effects in hierarchical data

---

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:

1. corrPrune() removes pairwise correlations > 0.7

2. modelPrune() 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)
head(names(bioclim_example))

# Response variable
summary(bioclim_example$species_richness)

The 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)))

# Which variables were kept?
head(attr(bio_clean, "selected_vars"), 10)

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:

1. Full model: All 19 bioclimatic variables (baseline)

2. After corrPrune: Variables with correlations < 0.7

3. 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)

Key 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"
)

---

Workflow 2: Survey Data Analysis

Goal: Reduce questionnaire length while preserving construct coverage and protecting key demographic variables.

Challenge: Survey instruments often contain redundant items within constructs (e.g., multiple satisfaction questions that are highly correlated). Reducing the number of items improves response rates and reduces respondent burden without losing measurement quality.

Strategy: Use corrPrune() with force_in to:

1. Ensure critical variables (like age) appear in the final model

2. Remove redundant Likert items within constructs

3. Maintain balanced representation across satisfaction, engagement, and loyalty domains

Data

data(survey_example)

# Data structure
dim(survey_example)
str(survey_example[, 1:10])  # First 10 columns

The dataset contains 200 respondents, 30 Likert-scale items (10 each for satisfaction, engagement, loyalty), plus demographics (age, gender, education) and an overall satisfaction outcome.

Prune with protected variables

We use force_in = "age" to ensure age remains in the analysis regardless of its correlation with other variables. This is useful when domain knowledge identifies theoretically important covariates that must not be removed.

# Exclude respondent_id, overall_satisfaction, and factor variables
survey_predictors <- survey_example[, !(names(survey_example) %in%
                                         c("respondent_id", "overall_satisfaction",
                                           "gender", "education"))]

# Convert ordered factors (Likert items 1-7) to numeric for correlation analysis
survey_numeric <- as.data.frame(lapply(survey_predictors, function(x) {
  if (is.ordered(x)) as.numeric(as.character(x)) else as.numeric(x)
}))

# Prune with protected variables
survey_clean <- corrPrune(
  data = survey_numeric,
  threshold = 0.6,
  force_in = "age"
)

# How many items remain?
cat(sprintf("Reduced from %d → %d variables\n",
            ncol(survey_numeric),
            ncol(survey_clean)))

# Which items were kept?
selected <- attr(survey_clean, "selected_vars")
print(selected)

The pruning reduced the questionnaire substantially while ensuring age was retained. The remaining items span all three constructs, avoiding the loss of entire domains.

Construct coverage

It's important to verify that pruning didn't eliminate entire constructs. We check how many items from each domain (satisfaction, engagement, loyalty) survived the correlation threshold.

# Count items per construct
satisfaction_kept <- sum(grepl("satisfaction_", selected))
engagement_kept <- sum(grepl("engagement_", selected))
loyalty_kept <- sum(grepl("loyalty_", selected))

cat(sprintf("Satisfaction: %d/10 items kept\n", satisfaction_kept))
cat(sprintf("Engagement: %d/10 items kept\n", engagement_kept))
cat(sprintf("Loyalty: %d/10 items kept\n", loyalty_kept))

Good balance: all three constructs retained representation, ensuring the reduced questionnaire still measures all intended dimensions.

Visualization

The barplots show (1) how many items survived pruning within each construct and (2) the overall variable reduction.

par(mfrow = c(1, 2))

# Items kept per construct
construct_data <- rbind(
  c(10, 10, 10),
  c(satisfaction_kept, engagement_kept, loyalty_kept)
)

barplot(construct_data,
        beside = TRUE,
        names.arg = c("Satisfaction", "Engagement", "Loyalty"),
        col = c("lightgray", "lightblue"),
        legend.text = c("Original (10)", "After pruning"),
        args.legend = list(x = "topright", bty = "n"),
        main = "Items per Construct",
        ylab = "Number of Items",
        ylim = c(0, 12))

# Percentage reduction
barplot(c(ncol(survey_numeric), ncol(survey_clean)),
        names.arg = c("Before", "After"),
        col = c("salmon", "lightgreen"),
        main = "Total Variables",
        ylab = "Count",
        ylim = c(0, max(ncol(survey_numeric)) * 1.2))
text(0.7, ncol(survey_numeric) + 1, ncol(survey_numeric), pos = 3)
text(1.9, ncol(survey_clean) + 1, ncol(survey_clean), pos = 3)

Model satisfaction

Now we fit a regression model predicting overall satisfaction from the pruned item set. Despite using fewer predictors, the model should maintain good explanatory power because we removed only redundant items.

# Add response back
survey_model_data <- data.frame(
  overall_satisfaction = survey_example$overall_satisfaction,
  survey_clean
)

# Fit regression model
model_survey <- lm(overall_satisfaction ~ ., data = survey_model_data)

# Summary
summary(model_survey)

Model comparison

Comparing the pruned model against the full 33-variable model shows that we retain most of the explanatory power (R²) while dramatically reducing model complexity.

# Full model (all 30 items + demographics)
full_survey_data <- data.frame(
  overall_satisfaction = survey_example$overall_satisfaction,
  survey_predictors
)

model_full_survey <- lm(overall_satisfaction ~ ., data = full_survey_data)

# Compare
data.frame(
  Model = c("Full (33 vars)", "Pruned (10 vars)"),
  R2 = c(summary(model_full_survey)$r.squared,
         summary(model_survey)$r.squared),
  Adj_R2 = c(summary(model_full_survey)$adj.r.squared,
             summary(model_survey)$adj.r.squared),
  Num_Predictors = c(33, 10)
)

---

Workflow 3: High-Dimensional Data

Goal: Reduce dimensionality in a gene expression dataset where the number of predictors far exceeds the number of samples (p >> n).

Challenge: With 200 genes and only 100 samples, exact enumeration of all maximal subsets becomes computationally expensive. Standard regression is also impossible due to rank deficiency.

Strategy: Use mode = "greedy" for fast, approximate pruning:

• The greedy algorithm scales linearly with the number of variables

• While not guaranteed to find the largest subset, it provides a high-quality solution orders of magnitude faster

• Ideal for exploratory analysis in high-dimensional settings

Data

data(genes_example)

# Data structure
dim(genes_example)

# Disease prevalence
table(genes_example$disease_status)

The dataset contains gene expression measurements for 200 genes across 100 samples, with a binary disease outcome. This is a classic p >> n scenario where regularization or dimensionality reduction is essential.

Greedy pruning

The greedy algorithm is a backward-elimination heuristic: it starts with all variables and iteratively removes the "worst" one -- prioritizing the variable with the most threshold violations, then the highest max correlation, then the highest average correlation -- until every remaining pair satisfies the threshold. This heuristic runs in O(p²) time compared to the exponential complexity of exact methods.

# Extract gene expression data (exclude ID and outcome)
gene_expr <- genes_example[, -(1:2)]

# Greedy pruning with timing (median of 3 runs)
greedy_times <- sapply(1:3, function(i) {
  system.time(corrPrune(
    data = gene_expr,
    threshold = 0.8,
    mode = "greedy"  # Fast for large p
  ))["elapsed"]
})
greedy_ms <- median(greedy_times) * 1000

genes_pruned <- corrPrune(data = gene_expr, threshold = 0.8, mode = "greedy")

# Reduction
cat(sprintf("Reduced from %d → %d genes (%.1f ms)\n",
            ncol(gene_expr),
            ncol(genes_pruned),
            greedy_ms))

The greedy algorithm completed in milliseconds while ensuring all pairwise correlations remain below 0.8.

Dimensionality reduction

# Barplot showing reduction
reduction_data <- c(ncol(gene_expr), ncol(genes_pruned))
barplot(reduction_data,
        names.arg = c("Original", "After Pruning"),
        main = "Gene Dimensionality Reduction",
        ylab = "Number of Genes",
        col = c("salmon", "lightblue"),
        ylim = c(0, max(reduction_data) * 1.2))
text(0.7, reduction_data[1] + 10, paste(reduction_data[1], "genes"), pos = 3)
text(1.9, reduction_data[2] + 10, paste(reduction_data[2], "genes\n(",
     round(100 * reduction_data[2] / reduction_data[1], 1), "% retained)"), pos = 3)

Exact vs greedy comparison

To demonstrate the speed advantage of the greedy algorithm, we benchmark both approaches on a smaller subset. The performance gap widens dramatically as the number of variables increases.

# Subset for comparison (use smaller subset for vignette build speed)
gene_subset <- gene_expr[, 1:20]  # Reduced from 50 to 20 for faster builds

# Benchmark exact mode (median of 3 runs)
exact_times <- sapply(1:3, function(i) {
  system.time(corrPrune(gene_subset, threshold = 0.8, mode = "exact"))["elapsed"]
})
exact_time <- median(exact_times) * 1000  # seconds -> milliseconds

# Benchmark greedy mode (median of 3 runs)
greedy_times <- sapply(1:3, function(i) {
  system.time(corrPrune(gene_subset, threshold = 0.8, mode = "greedy"))["elapsed"]
})
greedy_time <- median(greedy_times) * 1000  # seconds -> milliseconds

# Run once more to get actual results for comparison
exact_result <- corrPrune(gene_subset, threshold = 0.8, mode = "exact")
greedy_result <- corrPrune(gene_subset, threshold = 0.8, mode = "greedy")

# Compare
cat(sprintf("Exact mode: %d genes kept (%.1f ms)\n", ncol(exact_result), exact_time))
cat(sprintf("Greedy mode: %d genes kept (%.1f ms)\n", ncol(greedy_result), greedy_time))
cat(sprintf("Speedup: %.1fx faster\n", exact_time / greedy_time))

The greedy mode is substantially faster. For the full 200-gene dataset, exact enumeration would be prohibitively slow, while greedy mode completes in milliseconds.

Classification

Finally, we demonstrate that the pruned gene set is suitable for downstream classification, despite the dramatic dimensionality reduction.

# Prepare classification data
classification_data <- data.frame(
  disease_status = genes_example$disease_status,
  genes_pruned
)

# Logistic regression
model_genes <- glm(disease_status ~ .,
                   data = classification_data,
                   family = binomial())

# Prediction accuracy
predictions <- ifelse(predict(model_genes, type = "response") > 0.5,
                     "Disease", "Healthy")
accuracy <- mean(predictions == genes_example$disease_status)

cat(sprintf("Classification accuracy: %.1f%%\n", accuracy * 100))

---

Workflow 4: Mixed Models

Goal: Apply correlation-based pruning to fixed effects in a mixed-effects model with longitudinal data.

Challenge: Longitudinal data has hierarchical structure (observations nested within subjects, subjects nested within sites). Standard VIF calculations don't account for random effects, but we still need to control multicollinearity among fixed-effect predictors.

Strategy: Use modelPrune() with engine = "lme4":

• Only fixed effects are pruned based on VIF

• Random effects (1|subject) and (1|site) are preserved in the model formula

• This maintains the hierarchical structure while reducing collinearity

Note: This workflow requires the lme4 package and is shown with eval=FALSE for portability.

Data

data(longitudinal_example)

# Data structure
dim(longitudinal_example)
head(longitudinal_example)

# Study design
cat(sprintf("Subjects: %d\n", length(unique(longitudinal_example$subject))))
cat(sprintf("Sites: %d\n", length(unique(longitudinal_example$site))))
cat(sprintf("Observations per subject: %d\n",
            nrow(longitudinal_example) / length(unique(longitudinal_example$subject))))

The dataset has 500 observations from 50 subjects across 5 sites, with 10 measurements per subject. We have 5 correlated fixed-effect predictors (x1-x5).

Prune fixed effects

The modelPrune() function with engine = "lme4" respects the random-effects structure. Only the fixed effects (x1-x5) are candidates for removal; the random intercepts for subject and site remain untouched.

# Note: This example requires lme4 package
library(lme4)

# Define formula with random effects
# Note: Only fixed effects (x1-x5) will be pruned
#       Random effects (1|subject), (1|site) are preserved

pruned_mixed <- modelPrune(
  formula = outcome ~ x1 + x2 + x3 + x4 + x5 + (1|subject) + (1|site),
  data = longitudinal_example,
  engine = "lme4",
  limit = 5
)

# Which fixed effects were kept?
selected_fixed <- attr(pruned_mixed, "selected_vars")
cat("Fixed effects kept:\n")
print(selected_fixed)

# Which were removed?
removed_fixed <- attr(pruned_mixed, "removed_vars")
cat("\nFixed effects removed:\n")
print(removed_fixed)

The algorithm sequentially removes fixed effects with VIF > 5 until all remaining predictors satisfy the limit. The random effects structure is never modified.

Final model

The final model contains only the fixed effects that passed the VIF threshold, along with the original random effects.

final_mixed <- attr(pruned_mixed, "final_model")
summary(final_mixed)

VIF verification

We can manually verify that the pruning successfully reduced multicollinearity by comparing VIF values before and after pruning.

# Note: This example requires lme4 package
library(lme4)

# Fit full model
full_formula <- as.formula(paste("outcome ~",
                                 paste(paste0("x", 1:5), collapse = " + "),
                                 "+ (1|subject) + (1|site)"))

model_full_mixed <- lmer(full_formula, data = longitudinal_example)

# Extract fixed effects design matrices
X_full <- getME(model_full_mixed, "X")
X_pruned <- getME(final_mixed, "X")

# Compute VIF
compute_vif <- function(X) {
  X_scaled <- scale(X[, -1])  # Remove intercept
  sapply(seq_len(ncol(X_scaled)), function(i) {
    r2 <- summary(lm(X_scaled[, i] ~ X_scaled[, -i]))$r.squared
    1 / (1 - r2)
  })
}

vif_full <- compute_vif(X_full)
vif_pruned <- compute_vif(X_pruned)

# Compare
comparison_vif <- data.frame(
  Predictor = colnames(X_pruned)[-1],
  VIF_Before = vif_full,
  VIF_After = vif_pruned
)
print(comparison_vif)

All remaining predictors now have VIF < 5, indicating acceptable multicollinearity levels for mixed-effects modeling.

VIF reduction visualization

The plot below shows how pruning reduced VIF values across the retained fixed effects:

# Extract VIF values
predictors <- comparison_vif$Predictor
vif_before <- comparison_vif$VIF_Before
vif_after <- comparison_vif$VIF_After

# Set up bar positions
x <- seq_along(predictors)
width <- 0.35

# Create plot
par(mar = c(5, 4, 4, 2))
plot(
  x, vif_before,
  type = "n",
  xaxt = "n",
  xlab = "Fixed Effects",
  ylab = "VIF",
  main = "VIF Reduction After Pruning (Mixed Model)",
  ylim = c(0, max(vif_before) * 1.15)
)

# Add horizontal line at VIF = 5 threshold
abline(h = 5, col = "red", lty = 2, lwd = 2)

# Before bars (darker)
rect(
  x - width, 0,
  x, vif_before,
  col = rgb(0.8, 0.2, 0.2, 0.6),
  border = "white"
)

# After bars (lighter)
rect(
  x, 0,
  x + width, vif_after,
  col = rgb(0.2, 0.5, 0.8, 0.6),
  border = "white"
)

# Add x-axis labels
axis(1, at = x, labels = predictors, las = 2)

# Add legend
legend(
  "topright",
  legend = c("Before Pruning", "After Pruning", "VIF = 5 Threshold"),
  fill = c(rgb(0.8, 0.2, 0.2, 0.6), rgb(0.2, 0.5, 0.8, 0.6), NA),
  border = c("white", "white", NA),
  lty = c(NA, NA, 2),
  lwd = c(NA, NA, 2),
  col = c(NA, NA, "red"),
  bty = "o",
  bg = "white"
)

---

Summary

These four workflows demonstrate how corrselect integrates into diverse analytical pipelines:

Key takeaways:

1. Ecological modeling: Two-stage pruning (correlation → VIF) balances fit and stability

2. Survey analysis: force_in protects theoretically important variables while reducing redundancy

3. High-dimensional data: Greedy mode enables fast pruning when p >> n

4. Mixed models: VIF-based pruning respects hierarchical structure by targeting only fixed effects

Workflow selection guide:

| Scenario | Function | Key Parameters | |----------|----------|----------------| | Small to moderate p (< 50) | corrPrune() | mode = "auto" or "exact" | | High-dimensional (p >> n) | corrPrune() | mode = "greedy" | | Model-based refinement | modelPrune() | limit (VIF threshold) | | Protected variables | corrPrune() | force_in | | Mixed-effects models | modelPrune() | engine = "lme4" |

When to use each approach:

• corrPrune(): First-stage dimensionality reduction based purely on pairwise correlations

• modelPrune(): Second-stage refinement that accounts for model-specific multicollinearity (VIF)

• Combining both: Often yields best results (correlation pruning → VIF refinement → final model)

See Also

• vignette("quickstart") - Interface overview

• vignette("advanced") - Custom engines and algorithmic control

• vignette("comparison") - Comparison with alternatives

• vignette("theory") - Mathematical foundations

References

Thresholds:

• O'Brien, R. M. (2007). A caution regarding rules of thumb for variance inflation factors. Quality & Quantity, 41(5), 673-690.

• Dormann, C. F., et al. (2013). Collinearity: a review of methods to deal with it. Ecography, 36(1), 27-46.

Methods: - See package documentation and JOSS paper for algorithm details

Session Info

sessionInfo()

corrselect documentation built on July 19, 2026, 1:06 a.m.