Por que não estou obtendo uma distribuição uniforme do valor p da regressão logística com preditores aleatórios?

O código abaixo gera um conjunto de dados de teste que consiste em uma série de probabilidades de "sinal" com ruído binomial ao seu redor. O código então usa 5000 conjuntos de números aleatórios como séries "explicativas" e calcula o valor p da regressão logística para cada um.

Acho que as séries explicativas aleatórias são estatisticamente significativas no nível de 5% em 57% dos casos. Se você ler a parte mais longa do post abaixo, atribuo isso à presença do sinal forte nos dados.

Então, aqui está a pergunta principal: qual teste devo usar ao avaliar a significância estatística de uma variável explicativa quando os dados contiverem um sinal forte? O simples valor p parece ser bastante enganador.

Aqui está uma explicação mais detalhada do problema.

Estou intrigado com os resultados que recebo para os valores de p de regressão logística quando o preditor é realmente apenas um conjunto de números aleatórios. Meu pensamento inicial era que a distribuição dos valores-p deveria ser plana neste caso; o código R abaixo mostra, na verdade, um enorme aumento nos baixos valores de p. Aqui está o código:

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

O código gera dados de teste consistindo em ruído binomial em torno de um sinal forte e, em um loop, ajusta uma regressão logística dos dados contra um conjunto de números aleatórios e acumula os valores p dos preditores aleatórios; os resultados são exibidos como um histograma de valores p, os quantis reais do valor p para níveis de confiança de 1% e 5% e a taxa real de falso positivo correspondente a esses níveis de confiança.

Eu acho que uma razão para os resultados inesperados é que um preditor aleatório geralmente terá alguma correlação com o sinal verdadeiro e isso explica principalmente os resultados. No entanto, se você definir orthogonalPredictorpara TRUE, haverá zero correlação entre os preditores aleatórios e o sinal real, mas o problema ainda estará lá em um nível reduzido. Minha melhor explicação para isso é que, como o sinal verdadeiro não está em nenhum lugar do modelo sendo ajustado, o modelo é especificado incorretamente e os valores de p são suspeitos de qualquer maneira. Mas esse é um problema - quem já tem um conjunto de preditores disponíveis que se encaixa exatamente nos dados? Então, aqui estão algumas perguntas adicionais:

1. Qual é a hipótese nula precisa para regressão logística? Os dados são um ruído puramente binomial em torno de um nível constante (ou seja, não há um sinal verdadeiro)? Se você definir sd como 0 no código, não haverá sinal e o histograma parecerá plano.

2. A hipótese nula implícita no código é que o preditor não tem mais poder explicativo do que um conjunto de números aleatórios; é testado usando, digamos, o quantil empírico de 5% para o valor-p, conforme exibido pelo código. Existe uma maneira melhor de testar essa hipótese, ou pelo menos uma que não seja tão numericamente intensa?

------ Informação adicional

Esse código imita o seguinte problema: Taxas de inadimplência históricas mostram variação significativa ao longo do tempo (sinal) impulsionada por ciclos econômicos; as contagens padrão reais em um determinado momento são binomiais em torno dessas probabilidades padrão. Eu estava tentando encontrar variáveis ​​explicativas para o sinal quando suspeitei dos valores de p. Neste teste, o sinal é ordenado aleatoriamente ao longo do tempo, em vez de mostrar os ciclos econômicos, mas isso não deve importar para a regressão logística. Portanto, não há sobredispersão, o sinal realmente deve ser um sinal.

Existem várias questões aqui. Em particular, parece haver algumas confusões sobre como simular uma regressão logística padrão. Resumidamente, você não adiciona noise around... the probability "signal". Como resultado da maneira como você fez isso, há uma enorme quantidade de variabilidade nos dados 'binomiais' (esque) resultantes, muito mais do que deveria. Aqui estão as probabilidades no seu conjunto de dados:

plot(flips[,1]/rowSums(flips))

Se essas observações .4+ terminarem de um lado ou de outro, elas agirão como 'outliers' (na verdade não são) e gerarão um erro do tipo 1 em um modelo que não leva em consideração o fato de que esses dados não são realmente binomiais. Aqui está uma versão que usa um hack simples para permitir que o modelo detecte e leve em conta a sobredispersão:

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

Este é o resumo do modelo da última iteração. Note-se que a dispersão é estimada em$\approx 12\times$ o que deveria ser para um binômio verdadeiro:

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

Aqui está outra versão, onde eu me encaixo no mesmo modelo que você, mas apenas gere os dados sem o ruído adicional ao redor do sinal. (Observe que o código que é o mesmo é omitido por questões de brevidade.)

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] <- s$coefficients[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

A relação do código com o objetivo experimental é confusa. Você espera um preditor significativo para qualquer seleção de ortogonalPredictor ou sd? Eu não.

Com base na minha interpretação, parece que o experimento não está alinhado com o que estamos tentando testar. Onde o ruído é gerado, está implicitamente sendo repetidamente (ou seja, não aleatoriamente) anexado a observações individuais, o que fornece um sinal para a regressão.

Aqui está o que eu acho que foi planejado:

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_expanded$y <- 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))
}

As mudanças críticas são:

• Expandir o quadro de dados para o formato por observação em vez de um formato condensado (flips_expanded)
• Também experimentando um preditor correlacionado

Para nenhuma correlação entre y_actual e nosso preditor y:

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

E criando uma correlação bastante forte:

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