# Problem:

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.

# Question:

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

#set random seed
Random.seed!(2019)

#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)
shuffle!(true_b)
correct_position = findall(!iszero, true_b)

#simulate y
μ = X * true_b
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
end

#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

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