Package Repository: https://github.com/mobuobi/multipathaic
Documentation: See Function Reference for complete API details.
Traditional forward stepwise selection follows a single greedy path, making locally optimal choices at each step without considering alternative competitive models. This can lead to:
The multipathaic package addresses these limitations through a three-stage approach:
This vignette explains the methodology and algorithms underlying the package.
Find a set of plausible models that:
Standard forward selection follows a single greedy path:
Step 0: {} (null model)
Step 1: {x₃} (best single variable)
Step 2: {x₃, x₇} (best addition to x₃)
Step 3: {x₃, x₇, x₁} (best addition to x₃, x₇)
...Problem: If \(\{x_5\}\) was almost as good as \(\{x_3\}\) at Step 1, we never explore the \(x_5\) branch.
Instead of keeping only the best model, we maintain a frontier of competitive models at each step by:
This creates a forest of paths rather than a single tree.
Algorithm: build_paths(X, y, K, eps, delta, L)
Input: Data (X, y), parameters (K, eps, delta, L)
Output: Forest of model paths
1. Initialize:
F₀ = {∅} (frontier contains only null model)
2. For k = 1 to K:
a. Generate candidates:
For each model M in Fₖ₋₁:
For each variable x not in M:
Create candidate M ∪ {x}
Compute AIC(M ∪ {x})
b. For each parent separately, find its best finite child AIC.
Retain a child only if both:
AIC(parent) - AIC(child) >= eps
AIC(child) <= best child AIC for this parent + delta
c. Merge duplicate variable sets, then retain at most L lowest-AIC children.
This is the next frontier Fₖ.
d. Stop if no children survive; otherwise continue to the next step.
3. Return the intercept-only model and all retained frontiers.
Rejected candidates are not included in all_models.Delta (δ): Controls exploration breadth - Small δ (e.g., 0.5): Conservative, stays close to best path - δ = 2: Standard choice (a relative support guideline, not an equivalence test) - Large δ (e.g., 5): Aggressive, explores many alternatives
L: Controls computational cost - Small L (e.g., 25): Fast, less exploration - L = 50: Balanced (recommended) - Large L (e.g., 100): Thorough but slower
# Simple example to illustrate path exploration
library(multipathaic)
# Generate toy data
set.seed(42)
n <- 100
p <- 8
X <- matrix(rnorm(n * p), n, p)
colnames(X) <- paste0("x", 1:p)
y <- X[,1] + 0.5*X[,2] + rnorm(n)
# Build paths
forest <- build_paths(X, y, family = "gaussian", K = 4, delta = 2, L = 10)# In practice, you'd see output like:
cat("Step 1: 8 candidates → Keep 5 models within AIC* + 2\n")
#> Step 1: 8 candidates → Keep 5 models within AIC* + 2
cat("Step 2: 35 candidates → Keep 10 models within AIC* + 2\n")
#> Step 2: 35 candidates → Keep 10 models within AIC* + 2
cat("Step 3: 70 candidates → Keep 10 models within AIC* + 2\n")
#> Step 3: 70 candidates → Keep 10 models within AIC* + 2
cat("Step 4: 70 candidates → Keep 10 models within AIC* + 2\n")
#> Step 4: 70 candidates → Keep 10 models within AIC* + 2
cat("\nTotal unique models explored: 35\n")
#>
#> Total unique models explored: 35Even if a model has low AIC, it may be unstable - small changes to the data could lead to completely different variable selections. We want to identify variables that are reliably selected across data perturbations.
For each variable \(x_j\), define \(z_j^{(b)}\) as the fraction of models in the deepest retained frontier of bootstrap draw \(b\) that contain \(x_j\). The reported score is the mean of these fractions over successful draws:
\[\hat{\pi}_j = \frac{1}{B_{\mathrm{ok}}}\sum_{b \in \mathrm{successful}} z_j^{(b)}.\]
Rows are sampled with replacement. A successful intercept-only draw contributes zero for every variable. Failed draws are reported and omitted; an all-failed run stops with an error. A score of 0.6 means mean within-frontier inclusion of 0.6, not necessarily selection in 60% of draws. Earlier terminal branches are omitted by this definition; changing that choice would define a different estimator.
build_paths() (K, eps, delta, L)Algorithm: stability(X, y, K, eps, delta, L, B, resample_fraction)
Input: Data (X, y), path parameters, bootstrap parameters (B, resample_fraction)
Output: Stability scores π for each variable
1. Initialize:
count[j] = 0 for all variables j = 1, ..., p
2. For b = 1 to B:
a. Draw bootstrap sample:
(X_b, y_b) = resample(X, y, fraction = resample_fraction)
b. Run multi-path selection:
forest_b = build_paths(X_b, y_b, K, eps, delta, L)
c. On success, increment B_ok and compute each variable's fraction
across models in the deepest retained frontier (zero for a null-only draw).
Add those fractions to count[j].
On failure, record the error and do not increment B_ok.
3. If B_ok is zero, stop with an error. Otherwise π[j] = count[j] / B_ok.
4. Return π, B_ok, requested B, and failure messages.Why bootstrap? - Non-parametric: No distributional assumptions - Mimics sampling variation: Tests robustness to data perturbations - Well-established in statistical learning (Breiman, 1996)
Why resample_fraction < 1? - Creates more variation between resamples - Better tests stability under realistic data limitations - Standard choice: 0.8 (80% subsampling)
# Example stability scores
cat("Variable Stability Scores:\n")
#> Variable Stability Scores:
cat("x2: π = 0.94 (highly stable)\n")
#> x2: π = 0.94 (highly stable)
cat("x5: π = 0.86 (highly stable)\n")
#> x5: π = 0.86 (highly stable)
cat("x1: π = 0.72 (moderately stable)\n")
#> x1: π = 0.72 (moderately stable)
cat("x7: π = 0.58 (marginally stable)\n")
#> x7: π = 0.58 (marginally stable)
cat("x3: π = 0.34 (unstable)\n")
#> x3: π = 0.34 (unstable)
cat("x4: π = 0.12 (unstable)\n")
#> x4: π = 0.12 (unstable)We now have: 1. A forest of models with varying AICs (from Algorithm 1) 2. Stability scores for each variable (from Algorithm 2)
We want to filter this forest to identify plausible models that are both: - Near-optimal in terms of AIC - Built from stable variables
A model \(\mathcal{M}\) is plausible if:
build_paths() (all explored models)stability() (stability scores)Algorithm: plausible_models(forest, pi, Delta, tau)
Input: Forest of models, stability scores π, thresholds (Delta, tau)
Output: Subset of plausible models
1. Find best AIC:
AIC* = min{AIC(M) : M in forest}
2. For each model M in forest:
a. Check AIC criterion:
If AIC(M) > AIC* + Delta:
REJECT (too far from best AIC)
Continue to next model
b. Compute average stability:
π̄(M) = mean{π[j] : variable j in M}
c. Check stability criterion:
If π̄(M) < tau:
REJECT (built from unstable variables)
Continue to next model
d. If both criteria pass:
ADD M to plausible set
3. Return:
- Plausible models (sorted by AIC)
- Best AIC value (AIC*)
- Variable summary (inclusion frequency, mean stability)For each plausible model, we report: - Variables: List of predictors in the model - Size: Number of variables (model complexity) - AIC: Akaike Information Criterion - Avg Stability: Mean stability score of variables
For each variable appearing in plausible models: - Inclusion %: Percentage of plausible models containing this variable - Avg Stability: Mean stability score from bootstrap - Interpretation: Variables with high inclusion % and high stability are the most reliable predictors
cat("PLAUSIBLE MODELS (Delta=2, tau=0.6):\n\n")
#> PLAUSIBLE MODELS (Delta=2, tau=0.6):
cat("Best AIC: 2847.32\n\n")
#> Best AIC: 2847.32
cat("Model 1: [x2, x5, x8, x1] | Size=4 | AIC=2847.32 | Avg_Stab=0.83\n")
#> Model 1: [x2, x5, x8, x1] | Size=4 | AIC=2847.32 | Avg_Stab=0.83
cat("Model 2: [x2, x5, x8] | Size=3 | AIC=2848.15 | Avg_Stab=0.87\n")
#> Model 2: [x2, x5, x8] | Size=3 | AIC=2848.15 | Avg_Stab=0.87
cat("Model 3: [x2, x5, x1] | Size=3 | AIC=2849.01 | Avg_Stab=0.84\n\n")
#> Model 3: [x2, x5, x1] | Size=3 | AIC=2849.01 | Avg_Stab=0.84
cat("VARIABLE SUMMARY:\n\n")
#> VARIABLE SUMMARY:
cat("Variable | Inclusion% | Avg_Stability\n")
#> Variable | Inclusion% | Avg_Stability
cat("---------|------------|---------------\n")
#> ---------|------------|---------------
cat("x2 | 100% | 0.94\n")
#> x2 | 100% | 0.94
cat("x5 | 100% | 0.86\n")
#> x5 | 100% | 0.86
cat("x8 | 67% | 0.82\n")
#> x8 | 67% | 0.82
cat("x1 | 67% | 0.72\n")
#> x1 | 67% | 0.72Interpretation: - x2 and x5 are core predictors (100% inclusion, high stability) - x8 and x1 are important but less consistent (67% inclusion) - Choose Model 1 for best fit or Model 2 for parsimony
The multipath_aic() function combines all three algorithms into a single convenient call:
result <- multipath_aic(
X = X_train,
y = y_train,
family = "gaussian",
# Path building parameters
K = 10,
eps = 1e-6,
delta = 2,
L = 50,
# Stability parameters
B = 100,
resample_fraction = 0.8,
# Filtering parameters
Delta = 2,
tau = 0.6,
verbose = TRUE
)cat("result$forest # All explored models from build_paths()\n")
#> result$forest # All explored models from build_paths()
cat("result$stability # Variable stability scores from stability()\n")
#> result$stability # Variable stability scores from stability()
cat("result$plausible_models # Filtered plausible models\n")
#> result$plausible_models # Filtered plausible models
cat(" ├─ $plausible_models # Data frame of models passing filters\n")
#> ├─ $plausible_models # Data frame of models passing filters
cat(" ├─ $best_aic # Lowest AIC value\n")
#> ├─ $best_aic # Lowest AIC value
cat(" └─ $summary # Variable inclusion summary\n")
#> └─ $summary # Variable inclusion summary| Parameter | Role | Typical Values | Recommendation |
|---|---|---|---|
| K | Max steps | 5-20 | min(p, 10) |
| eps | Stop threshold | 1e-8 to 1e-4 | 1e-6 |
| delta | Path branching | 1-4 | 2 |
| L | Frontier size | 25-100 | 50 |
| B | Bootstrap samples | 50-200 | 100 |
| resample_fraction | Bootstrap size | 0.6-0.9 | 0.8 |
| Delta | AIC filter | 1-4 | 2 |
| tau | Stability filter | 0.5-0.8 | 0.6 |
More exploration (finds more alternatives, slower): - Increase delta (e.g., 3-4) - Increase L (e.g., 75-100) - Increase B (e.g., 150-200)
Faster computation (more focused, quicker): - Decrease delta (e.g., 1) - Decrease L (e.g., 25) - Decrease B (e.g., 50)
High-dimensional problems (p > n): - Limit K to 5-8 (prevent overfitting) - Increase tau to 0.7-0.8 (stricter stability) - Use delta = 2 (moderate exploration)
Highly correlated predictors: - Increase delta to 3-4 (more path diversity) - Increase B to 100-200 (better stability estimates) - May need to lower tau to 0.5 (relaxed stability)
Small sample size (n < 100): - Use resample_fraction = 0.9 (larger bootstrap samples) - Reduce B to 50 (each resample is more similar) - Limit K to 5-7 (prevent overfitting)
The Akaike Information Criterion balances model fit with complexity:
\[\text{AIC} = -2 \log L(\hat{\theta}) + 2k\]
Where: - \(L(\hat{\theta})\) = maximized likelihood - \(k\) = number of parameters
Advantages: - Has asymptotic connections to predictive risk under regularity conditions; it does not replace held-out validation - Penalizes overfitting - Allows direct model comparison
Burnham & Anderson (2002) provide evidence-based guidelines:
Where \(\Delta_i = \text{AIC}_i - \text{AIC}_{\text{min}}\)
Using \(\Delta = 2\) retains only models with substantial empirical support.
Bootstrap resampling (Efron, 1979; Breiman, 1996): - Non-parametric: No distributional assumptions required - Mimics sampling variation: Tests robustness to data perturbations - Well-studied: Theoretical properties well-understood - Practical: Computationally feasible with modern hardware
High scores describe stability under this particular resampling and search procedure. They do not prove that a variable is causally important or control false discoveries. Strongly correlated predictors can substitute for each other, splitting their scores.
This is an exploratory default for the mean variable score in a model. A model with variable scores 1.0 and 0.2 passes at 0.6; individual variables need not exceed the threshold. Intercept-only models are exempt and have undefined average stability. There is no demonstrated universal optimal threshold or robustness guarantee. Assess sensitivity and validate on held-out observations.
| Aspect | Traditional Forward | Multi-Path AIC |
|---|---|---|
| Exploration | Single greedy path | Multiple competitive paths |
| Stability | Not assessed | Quantified via bootstrap |
| Robustness | High variance | Identifies stable predictors |
| Alternatives | None | Multiple plausible models |
| Aspect | All-Subsets | Multi-Path AIC |
|---|---|---|
| Feasibility | \(2^p\) models (infeasible for large p) | Explores \(\sim K \times L\) models |
| Stability | Not assessed | Built-in stability estimation |
| Computation | Exponential in p | Linear in K, manageable with L |
| Aspect | Regularization | Multi-Path AIC |
|---|---|---|
| Interpretability | Continuous shrinkage | Discrete variable selection |
| Stability | Depends on tuning | Explicitly measured |
| Multiple solutions | Single λ chosen | Multiple plausible models |
| Inference | Post-selection inference complex | Clear model sets for reporting |
result <- multipath_aic(X, y, family = "gaussian")tau or increase Deltatau or decrease DeltaL or BWhen publishing analyses using this method, report:
Currently supports: - Gaussian (linear regression) - Binomial (logistic regression)
Future extensions could include: - Poisson, negative binomial (count data) - Survival models (Cox regression) - Multinomial (multi-class classification)
Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716-723.
Breiman, L. (1996). Bagging predictors. Machine Learning, 24(2), 123-140.
Burnham, K. P., & Anderson, D. R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach (2nd ed.). Springer-Verlag.
Efron, B. (1979). Bootstrap methods: Another look at the jackknife. The Annals of Statistics, 7(1), 1-26.
Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.). Springer.
For Gaussian linear models:
\[\text{AIC} = n \left[\log(2\pi) + 1 + \log\left(\frac{\text{RSS}}{n}\right)\right] + 2k\]
Where: - RSS = residual sum of squares - n = sample size - k = number of estimated regression coefficients (including intercept) plus one for the Gaussian variance; this formula is for unweighted fits
For binomial logistic models:
\[\text{AIC} = -2 \sum_{i=1}^n \left[y_i \log(\hat{p}_i) + (1-y_i)\log(1-\hat{p}_i)\right] + 2k\]
Let \(F_b\) be the deepest retained frontier in successful draw \(b\).
\[z_j^{(b)} = \frac{1}{|F_b|}\sum_{M\in F_b}\mathbb{1}(j\in M),\qquad \hat{\pi}_j = \frac{1}{B_{\mathrm{ok}}}\sum_{b\in\mathrm{successful}}z_j^{(b)}.\]
For a null-only draw use \(z_j^{(b)}=0\). Increasing the number of independent draws reduces Monte Carlo error for this conditional resampling quantity; it does not recover a probability that a variable is truly important. Frequent failures can bias the successful-draw distribution and should be investigated.
For a model \(\mathcal{M} = \{j_1, j_2, \ldots, j_k\}\):
\[\bar{\pi}(\mathcal{M}) = \frac{1}{k} \sum_{j \in \mathcal{M}} \pi_j\]
This measures the average reliability of the variables in the model.
Package Version: 0.1.0
Last Updated: December 2024
Maintainer: Michael Obuobi
GitHub: https://github.com/mobuobi/multipathaic