Why am I not getting a uniform p-value distribution from logistic regression with random predictors?

The code below generates a set of test data consisting of a series of "signal" probabilities with binomial noise around it. The code then uses 5000 sets of random numbers as "explanatory" series and calculates the logistic regression p-value for each.

I find that the random explanatory series are statistically significant at the 5% level in 57% of the cases. If you read through the longer part of the post below I attribute this to the presence of the strong signal in the data.

So here's the primary question: what test should I use when evaluating the statistical significance of an explanatory variable when the data contains a strong signal? The simple p-value seems to be quite misleading.

Here's a more detailed explanation of the issue.

I'm puzzled by results I get for logistic regression p-values when the predictor is actually just a set of random numbers. My initial thought was that the distribution of p-values should be flat in this case; the R code below actually shows a huge spike at low p-values. Here's the code:

set.seed(541713)

lseries <-   50
nbinom  <-  100
ntrial  <- 5000
pavg    <- .1  # median probability
sd      <- 0 # data is pure noise
sd      <- 1 # data has a strong signal
orthogonalPredictor <- TRUE   # random predictor that is orthogonal to the true signal
orthogonalPredictor <- FALSE  # completely random predictor

qprobs  <- c(.05,.01) # find the true quantiles for these p-values

yactual <- sd * rnorm(lseries)  # random signal
pactual <- 1 / (1 + exp(-(yactual + log(pavg / (1-pavg)))))
heads   <- rbinom(lseries, nbinom, pactual)
## test data, binomial noise around pactual, the probability "signal"
flips   <- cbind(heads, nbinom - heads)
# summary(glm(flips ~ yactual, family = "binomial"))

pval <- numeric(ntrial)

for (i in 1:ntrial){
yrandom <- rnorm(lseries)
if (orthogonalPredictor){ yrandom <- residuals(lm(yrandom ~ yactual)) }
s       <- summary(glm(flips ~ yrandom, family="binomial"))
pval[i] <- s$coefficients[2,4] } hist(pval, breaks=100) print(quantile(pval, probs=c(.01,.05))) actualCL <- sapply(qprobs, function(c){ sum(pval <= c) / length(pval) }) print(data.frame(nominalCL=qprobs, actualCL))  The code generates test data consisting of binomial noise around a strong signal, then in a loop fits a logistic regression of the data against a set of random numbers, and accumulates the p-values of the random predictors; results are displayed as a histogram of p-values, the actual p-value quantiles for confidence levels of 1% and 5%, and the actual false positive rate corresponding to those confidence levels. I think one reason for the unexpected results is that a random predictor will generally have some correlation with the true signal and this mostly accounts for the results. However if you set orthogonalPredictor to TRUE there will be zero correlation between the random predictors and actual signal, but the problem is still there at a reduced level. My best explanation for that is that since the true signal is not anywhere in the model being fitted then the model is misspecified and p-values are suspect anyway. But this is a catch-22 - who ever has a set of predictors available that exactly fits the data? So here are some additional questions: 1. What is the precise null hypothesis for logistic regression? Is it that the data is purely binomial noise around a constant level (i.e., there isn't a true signal)? If you set sd to 0 in the code then there is no signal and the histogram does look flat. 2. The implicit null hypothesis in the code is that the predictor has no more explanatory power than a set of random numbers; it's tested by using say the empirical 5% quantile for the p-value as displayed by the code. Is there a better way of testing this hypothesis, or at least one that's not so numerically intensive? ------ Additional information This code mimics the following problem: Historical default rates show significant variation over time (signal) driven by economic cycles; the actual default counts at a given point in time are binomial around these default probabilities. I was trying to find explanatory variables for the signal when I became suspicious of the p-values. In this test the signal is randomly ordered over time rather than showing the economic cycles, but that shouldn't matter to logistic regression. So there is no overdispersion, the signal is really meant to be a signal. • This is interesting, but I don't think it will attract a lot of interest because it's 50% code. Could you intersperse the code with written explanations, formulas, etc? This is not Stack Overflow and nobody will debug your code for you. Also, definitely eliminate at least 2 of your 4 questions. Policy here is one question for post: two may be tolerable, 4 is like making sure you'll never get an answer. – DeltaIV Dec 4 '17 at 10:51 • Thanks, I've never used Cross Validated before. I'll make some edits. – Ron Hylton Dec 6 '17 at 1:36 2 Answers There are several issues here. In particular, there seem to be some confusions about how to simulate a standard logistic regression. Briefly, you don't add noise around... the probability "signal". As a result of the way you did this, there is a huge amount of variability in the resulting 'binomial'(-esque) data, way more than there should be. Here are the probabilities in your dataset: plot(flips[,1]/rowSums(flips))  If those .4+ observations end up on one side or the other, they will act as 'outliers' (they aren't really) and drive a type 1 error in a model that doesn't take into account the fact that these data aren't really binomial. Here is a version that uses a simple hack to allow the model to detect and account for overdispersion: set.seed(5082) pval <- numeric(ntrial) for (i in 1:ntrial){ yrandom <- rnorm(lseries) s <- summary(glm(flips ~ yrandom, family="quasibinomial")) # changed family pval[i] <- s$coefficients[2,4]
}

hist(pval, breaks=100)
print(quantile(pval, probs=c(.01,.05)))
#          1%          5%
# 0.006924617 0.046977246
actualCL <- sapply(qprobs, function(c){ sum(pval <= c) / length(pval) })
print(data.frame(nominalCL=qprobs, actualCL))
#   nominalCL actualCL
# 1      0.05   0.0536
# 2      0.01   0.0128


This is the model summary from the last iteration. Note that the dispersion is estimated to be $\approx 12\times$ what it should be for a true binomial:

s
# Call:
# glm(formula = flips ~ yrandom, family = "quasibinomial")
#
# Deviance Residuals:
#    Min      1Q  Median      3Q     Max
# -5.167  -2.925  -1.111   1.101   8.110
#
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept) -1.96910    0.14942 -13.178   <2e-16 ***
# yrandom     -0.02736    0.14587  -0.188    0.852
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# (Dispersion parameter for quasibinomial family taken to be 11.97867)
#
#     Null deviance: 532.38  on 49  degrees of freedom
# Residual deviance: 531.96  on 48  degrees of freedom
# AIC: NA
#
# Number of Fisher Scoring iterations: 5


Here is another version, where I fit the same model that you do, but just generate the data without the added noise around the signal. (Note that code that is otherwise the same is omitted for brevity.)

set.seed(541713)
...
pactual <- 1 / (1 + exp(-(log(pavg / (1-pavg)))))  # deleted yactual
...
for (i in 1:ntrial){
yrandom <- rnorm(lseries)
if (orthogonalPredictor){ yrandom <- residuals(lm(yrandom ~ yactual)) }
s       <- summary(glm(flips ~ yrandom, family="binomial"))
pval[i] <- scoefficients[2,4] } hist(pval, breaks=100) print(quantile(pval, probs=c(.01,.05))) # 1% 5% # 0.01993318 0.07027473 actualCL <- sapply(qprobs, function(c){ sum(pval <= c) / length(pval) }) print(data.frame(nominalCL=qprobs, actualCL)) # nominalCL actualCL # 1 0.05 0.0372 # 2 0.01 0.0036  • I've added an explanation at the bottom of the original post of the actual problem that the code is mimicking. The signal really is a signal, not overdispersion, and the actual problem really is binomial noise around the signal. – Ron Hylton Dec 7 '17 at 1:57 • @RonHylton, this is a misunderstanding. Your data really are overdispersed relative to the binomial distribution. You need to account for that somehow (there are several ways--using the quasibinomial is just the easiest). If you account for the overdispersion that is clearly there (this is true of how you are generating your data from the standpoint of statistical theory, & it is empirically demonstrated as well), then you no longer have type 1 error inflation. – gung - Reinstate Monica Dec 7 '17 at 2:11 • @RonHylton, you have no treatment effect (& you have no groups in this simulation--only a continuous predictor, yrandom), that's why these are type 1 errors, which you are rightly worried about. There is considerably more variability in the conditional response distributions than there is under a binomial. That is the definition of overdispersion. To be more specific, your conditional distributions are beta-binomial, not binomial. – gung - Reinstate Monica Dec 8 '17 at 2:01 • FWIW, the data generating mechanism in your code probably does not match the situation you describe & would be better addressed with a heteroscedasticity and autocorrelation consistent sandwich error than by using family=quasibinomial. Most likely you don't only have overdispersion, but probably errors that are correlated over time. – gung - Reinstate Monica Dec 8 '17 at 2:02 • @RonHylton, I'm not sure what to say here. Your code simulates a null--there is no signal. You have type 1 error inflation b/c you have overdispersion & the model doesn't account for that; if you do, the type 1 error inflation goes away. Your code does not match your description here, or in your edit; in truth, your code is not the normal way simulations would be set up. I don't follow your thinking about how your description or code relate to a GLMM. All I can say is that you have overdispersion & a "standard logistic regression" used in the "perfectly normal way" doesn't account for that. – gung - Reinstate Monica Dec 10 '17 at 2:43 The code's relation to the experimental goal is confusing. Do you expect a significant predictor for any selection of orthogonalPredictor or sd? I do not. Based on my interpretation, it looks like the experiment does not align with what we are trying to test. Where the noise is generated, it's implicitly being repeatedly (i.e. non-randomly) attached to individual observations, which provides a signal to the regression. Here's what I think was intended: lseries <- 50 nbinom <- 100 ntrial <- 5000 pavg <- .1 # median probability run_experiment <- function(sd = 0, orthogonalPredictor = FALSE, predictor_noise_sd = NA) { qprobs <- c(.05,.01) # find the true quantiles for these p-values yactual <- sd * rnorm(lseries) # random signal pactual <- 1 / (1 + exp(-(yactual + log(pavg / (1-pavg))))) heads <- rbinom(lseries, nbinom, pactual) ## test data, binomial noise around pactual, the probability "signal" flips_expanded <- rbind(data.frame(flip_result = rep(rep(1, length(heads)), heads), y_actual = rep(yactual, heads)), data.frame(flip_result = rep(rep(0, length(heads)), nbinom-heads), y_actual = rep(yactual, nbinom-heads)) ) summary(glm(flip_result ~ y_actual, flips_expanded, family = "binomial")) pval <- numeric(ntrial) for (i in 1:ntrial){ flips_expandedy <- rnorm(nrow(flips_expanded))
if (orthogonalPredictor){ flips_expanded$y <- residuals(lm(y ~ y_actual, flips_expanded)) } if (!is.na(predictor_noise_sd)) {flips_expanded$y <- rnorm(nrow(flips_expanded), flips_expanded$y_actual, predictor_noise_sd)} s <- summary(glm(flip_result ~ y, flips_expanded, family="binomial")) pval[i] <- s$coefficients[2,4]
}

hist(pval, breaks=100)
print(quantile(pval, probs=c(.01,.05)))
actualCL <- sapply(qprobs, function(c){ sum(pval <= c) / length(pval) })
print(data.frame(nominalCL=qprobs, actualCL))
}


The critical changes are:

• Expanding the data frame to per-observation format instead of a condensed format (flips_expanded)
• Also experimenting with a correlated predictor

For no correlation between y_actual & our predictor y:

> run_experiment()
1%         5%
0.01077116 0.05045712
nominalCL actualCL
1      0.05   0.0496
2      0.01   0.0096


And creating a fairly strong correlation:

> run_experiment(1,FALSE,10)
1%           5%
0.0001252817 0.0019125482
nominalCL actualCL
1      0.05   0.3002
2      0.01   0.1286


• This isn't correct. The problem isn't that the rows are grouped, that can be just fine. Moreover, the OP already knows that it works under some combinations of parameter settings. With sd =1 & orthogonalPredictor = FALSE, it (naively) shouldn't work because the x-variable (yrandom) is randomly generated & not related to the y-variable. – gung - Reinstate Monica Dec 6 '17 at 14:23
• If there isn't a problem with the way noise is added to grouped rows, why do I find the uniform p-value distribution when ungrouping and adding noise per observation? It doesn't seem to me that that was what the OP was trying to experiment on. – khol Dec 6 '17 at 15:36
• Overdispersion is undetectable with ungrouped data. If you ungrouped the data, though, you would need to indicate somehow that various observations come from the same unit (patient--assuming these were real data). – gung - Reinstate Monica Dec 6 '17 at 15:49
• I've added an explanation at the bottom of the original post of the actual problem that the code is mimicking. I think this should make the experimental goal clearer. – Ron Hylton Dec 7 '17 at 1:59
• Understood now that I see what you're trying to simulate, this is not the correct answer. – khol Dec 7 '17 at 3:10