Poisson 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 variable selection to various regression models including Poisson regression for count data analysis.

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

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

Simulate data

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

set.seed(42)
n <- 1000     # sample size
p <- 100      # number of predictors
L <- 10       # number of single-effect components
block_size <- 10
n_blocks <- p / block_size
rho <- 0.98   # 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, -0.5, 0.5, -0.5)

# Generate count response via Poisson model
linear_pred <- drop(-1 + X %*% theta_true)  # with intercept for reasonable counts
lambda <- exp(linear_pred)  # log link
y <- rpois(n, lambda)

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 -0.5 0.5 -0.5
cat("Response summary:\n")
## Response summary:
print(summary(y))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.953   1.000  39.000
cat("Mean count:", mean(y), "\n")
## Mean count: 0.953
cat("Variance:", var(y), "\n")
## Variance: 6.287078

Fit gSuSiE model

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

# Load glmsusie library
library(glmsusie)

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

summary(fit)
## 
## Call:
## glmsusie(X = X, y = y, L = L, family = poisson())
## 
## Family: poisson 
## 
## Coefficients: (sorted by PIP)
##         Estimate    PIP
## X3   1.071676622 1.0000
## X47  0.528881765 0.9940
## X18 -0.503746091 0.9274
## X82 -0.283454779 0.7249
## X88 -0.025207794 0.0641
## X83 -0.021628894 0.0536
## X85 -0.016117410 0.0429
## X89 -0.014803722 0.0397
## X12 -0.014921813 0.0276
## X81 -0.009984642 0.0266
## ... (90 more coefficients not shown)
## 
## 95% Confidence Sets:
##                          Set Coverage
## cs1 {81, 82, 83, 85, 88, 89}   0.9519
## cs2                      {3}   1.0000
## cs3                     {47}   0.9940
## cs4                 {12, 18}   0.9550
## 
## Model converged after 3 iterations.
## Computation time: 15.8 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")
  cat("Coverage rate:", length(covered_vars), "/", length(true_vars), "\n")
} else {
  cat("No credible sets identified\n")
}
## Number of credible sets: 4 
## CS 1 contains variables: 81 82 83 85 88 89 
## CS 2 contains variables: 3 
## CS 3 contains variables: 47 
## CS 4 contains variables: 12 18 
## 
## True variables covered by credible sets: 82 3 47 18 
## Coverage rate: 4 / 4

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      1.0716766 0.9999998
## X18       18      -0.5     -0.5037461 0.9273998
## X47       47       0.5      0.5288818 0.9939535
## X82       82      -0.5     -0.2834548 0.7248531

Conclusion

This example demonstrates that gSuSiE effectively handles Poisson 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 Poisson family implementation successfully recovers the sparse signal structure in count data while properly handling the log-linear relationship between predictors and the rate parameter.