I have an $n \times p$ design matrix $\mathbf{X}$ where $p > n$ and a response vector $\mathbf{y}$ that are count data, which suggests using the Poisson GLM model. I suspect only a few predictors actually affect the response, and I want to determine what they are.

Partial Solution:

One solution is to use something like LASSO. Another is to fit a bunch of univariate regressions, and select predictors based on a p-value cutoff (as people usually do in my field). The latter approach is typically conservative if we use $0.05 / p$ as the cutoff.


I was testing the latter approach (i.e. fit a bunch of univariate regressions) with simulated data, and realized that it is very non-conservative in that more than 10% of all predictors are "significant". Doesn't this violate the conservative nature of the Bonferroni correction? Below I performed a simulation with $p=10000$ predictors where only $10$ predictor affects the response. After all the univariate analysis, I got over $1000$ "causative" predictors. If I use do normal or logistic regression, the problem vanishes. So it seems to be a problem particular to Poisson data. What is happening here?

Minimal code example (in Julia) to reproduce my claim:

#sample size, number of covariates, number of true causative predictors, Distribution, Link
n = 1000
p = 10000
k = 10
d = Poisson
l = LogLink()

#set random seed

#random normal design matrix
X = randn(n, p)

#define true model with non-zero entries 0.1, 0.2, ..., 1.0
true_b = zeros(p)
true_b[1:10] .= collect(0.1:0.1:1.0)
correct_position = findall(!iszero, true_b)

#simulate y
μ = X * true_b
prob = linkinv.(l, μ)
clamp!(prob, -20, 20)
y = [rand(d(i)) for i in prob]
y = Float64.(y)

1000-element Array{Float64,1}: 13.0 0.0 0.0 2.0 1.0 0.0 8.0 1.0 0.0 3.0 0.0 2.0 0.0 ⋮
7.0 1.0 15.0 2.0 14.0 0.0 8.0 2.0 11.0 15.0 16.0 14.0

#now perform univariate tests to test for association
data = DataFrame(X=zeros(n), Y=y)
placeholder = zeros(n)
pvalues = zeros(p)
for i in 1:p
    copyto!(placeholder, @view(X[:, i]))
    data.X .= placeholder
    result = glm(@formula(Y ~ X), data, Poisson(), LogLink())
    pvalues[i] = coeftable(result).cols[4][2].v

#significance = bonferonni correction
significance = 0.05 / p
passing_position = findall(pvalues .<= significance)

#determine whether the 10 causative predictors were found
marginal_found = [correct_position[i] in passing_position for i in 1:k]

10-element Array{Bool,1}: true true true false true false true true true true

#compute some summary statistic
true_positive = count(!iszero, true_b[passing_snps_position2])
false_positive = length(passing_snps_position2) - juliaglm_tp
false_negative = k - juliaglm_tp

@show true_positive
@show false_positive
@show false_negative

true_positive = 8

false_positive = 1386

false_negative = 2


1 Answer 1


I found the solution a long time ago. The main reason for this is because we have too many $0$'s in the response vector. This is expected because a lot of the mean parameter $\mu_i = e^{x^t_i\beta}$ is small. Excess $0$'s creates overinflated variance, and a standard Poisson model suffers because it assumes its mean is equal to its variance.

One solution would be to use a zero-inflated Poisson regression, which is what I ended up using. The R package pscl is a good choice for fitting, which works nicely and only gave me 5~6 significant predictors (much better than >1000 predictors as provided by a Poisson model).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.