Logistic Regression with Correlated Predictors

Yizeng Li

2025-10-28

Introduction

The glmsusie package implements the generalized sum of single effects (gSuSiE) model, which represents the overall effect as a sum of a small number of single-effect components. This approach extends the SuSiE method to general regression models including logistic regression.

In this vignette, we demonstrate gSuSiE’s performance with binary outcomes and highly correlated predictors:

• Simulate a block-wise correlation structure with 100 predictors
• Apply gSuSiE to perform variable selection via glmsusie() with binomial family
• Visualize coefficient estimates, posterior inclusion probabilities (PIPs), and credible sets (CSs)
• Evaluate predictive performance for binary classification

Simulate data

We generate n=2500 observations and p=100 predictors with block-wise correlation. Every 10 consecutive variables are highly correlated (\rho=0.95) within each block, while blocks are independent. Only 4 variables have nonzero effects, located in different correlation blocks.

set.seed(42)
n <- 2500     # sample size
p <- 100      # number of predictors
L <- 10       # number of single-effect components
block_size <- 10
n_blocks <- p / block_size
rho <- 0.95    # within-block correlation

# Create block-wise correlation matrix
Sigma <- matrix(0, p, p)
for (b in 1:n_blocks) {
  block_idx <- ((b-1)*block_size + 1):(b*block_size)
  # Within-block correlation matrix
  block_corr <- matrix(rho, block_size, block_size)
  diag(block_corr) <- 1
  Sigma[block_idx, block_idx] <- block_corr
}

# Generate correlated predictors
X <- MASS::mvrnorm(n, mu = rep(0, p), Sigma = Sigma)

# True sparse coefficients (one per block, spread across different blocks)
theta_true <- rep(0, p)
theta_true[c(3, 18, 47, 82)] <- c(1, -1, 1, -1)

# Generate binary response via logistic model
linear_pred <- drop(-3 + X %*% theta_true) # with intercept
prob <- plogis(linear_pred)  # inverse logit
y <- rbinom(n, 1, prob)

cat("True nonzero coefficients at positions:", which(theta_true != 0), "\n")
## True nonzero coefficients at positions: 3 18 47 82
cat("True coefficient values:", theta_true[theta_true != 0], "\n")
## True coefficient values: 1 -1 1 -1
cat("Response proportion (positives):", mean(y), "\n")
## Response proportion (positives): 0.1296

Fit gSuSiE model

We allow up to L=10 single effects and use the binomial family for logistic regression.

library(glmsusie)

# Model fitting
fit <- glmsusie(
  X      = X,
  y      = y,
  L      = L,
  family = binomial()
)

summary(fit)
## 
## Call:
## glmsusie(X = X, y = y, L = L, family = binomial())
## 
## Family: binomial 
## 
## Coefficients: (sorted by PIP)
##        Estimate    PIP
## X47  0.94628134 0.9990
## X18 -0.90801671 0.9802
## X82 -0.73486720 0.9109
## X7   0.24976181 0.4268
## X4   0.17222280 0.3477
## X2   0.12770384 0.2657
## X8   0.12946897 0.2530
## X5   0.05700890 0.1331
## X10  0.05091489 0.1236
## X3   0.04815316 0.1183
## ... (90 more coefficients not shown)
## 
## 95% Confidence Sets:
##                                 Set Coverage
## cs1                        {81, 82}   0.9743
## cs2 {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}   0.9737
## cs3                            {47}   0.9990
## cs4                            {18}   0.9801
## cs5          {2, 3, 4, 5, 7, 8, 10}   0.9597
## 
## Model converged after 3 iterations.
## Computation time: 18.97 seconds.

Results

Coefficient estimates

The method successfully identifies the true signal locations despite high correlation within blocks:

plot(fit, which = "coefficients")

Posterior inclusion probabilities

Shows the probability that each variable is included in any single-effect component:

plot(fit, which = "probabilities")

95% credible sets

Each credible set contains variables where at least one is likely active with 95% confidence. Note how the method handles correlated variables within blocks:

plot(fit, which = "sets")

Credible sets summary

# Summary of credible sets
cs_sets <- fit$cs$sets
if (length(cs_sets) > 0) {
  cat("Number of credible sets:", length(cs_sets), "\n")
  for (i in seq_along(cs_sets)) {
    cat("CS", i, "contains variables:", cs_sets[[i]], "\n")
  }
  
  # Check coverage of true variables
  true_vars <- which(theta_true != 0)
  covered_vars <- intersect(unlist(cs_sets), true_vars)
  cat("\nTrue variables covered by credible sets:", covered_vars, "\n")
} else {
  cat("No credible sets identified\n")
}
## Number of credible sets: 5 
## CS 1 contains variables: 81 82 
## CS 2 contains variables: 1 2 3 4 5 6 7 8 9 10 
## CS 3 contains variables: 47 
## CS 4 contains variables: 18 
## CS 5 contains variables: 2 3 4 5 7 8 10 
## 
## True variables covered by credible sets: 82 3 47 18

Coefficient comparison

Compare estimated coefficients with true values:

# Extract coefficients
coef_est <- coef(fit)
comparison <- data.frame(
  Variable = 1:p,
  True_Coef = theta_true,
  Estimated_Coef = coef_est,
  PIP = fit$pip
)

# Show results for true active variables
active_comparison <- comparison[theta_true != 0, ]
print(active_comparison)
##     Variable True_Coef Estimated_Coef       PIP
## X3         3         1     0.04815316 0.1183015
## X18       18        -1    -0.90801671 0.9801525
## X47       47         1     0.94628134 0.9990470
## X82       82        -1    -0.73486720 0.9108888

Conclusion

This example demonstrates that gSuSiE effectively handles logistic regression with highly correlated predictors. The method identifies relevant variables even when they are embedded within correlation blocks, providing both point estimates and uncertainty quantification through credible sets. The binomial family implementation successfully recovers the sparse signal structure in binary outcomes.