Warum erhalte ich keine gleichmäßige p-Wert-Verteilung durch logistische Regression mit zufälligen Prädiktoren?

Der folgende Code generiert einen Satz von Testdaten, die aus einer Reihe von "Signal" -Wahrscheinlichkeiten mit Binomialrauschen bestehen. Der Code verwendet dann 5000 Sätze von Zufallszahlen als "erklärende" Reihen und berechnet den logistischen Regressions-p-Wert für jede.

Ich finde, dass die zufälligen Erklärungsreihen in 57% der Fälle bei 5% statistisch signifikant sind. Wenn Sie den längeren Teil des Beitrags unten durchlesen, schreibe ich dies auf das Vorhandensein des starken Signals in den Daten zurück.

Hier ist also die Hauptfrage: Welchen Test sollte ich verwenden, um die statistische Signifikanz einer erklärenden Variablen zu bewerten, wenn die Daten ein starkes Signal enthalten? Der einfache p-Wert scheint ziemlich irreführend zu sein.

Hier finden Sie eine detailliertere Erläuterung des Problems.

Ich bin verwirrt über die Ergebnisse, die ich für logistische Regressions-p-Werte erhalte, wenn der Prädiktor eigentlich nur eine Reihe von Zufallszahlen ist. Mein erster Gedanke war, dass die Verteilung der p-Werte in diesem Fall flach sein sollte; Der R-Code unten zeigt tatsächlich eine große Spitze bei niedrigen p-Werten. Hier ist der 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))

Der Code erzeugt Testdaten, die aus Binomialrauschen um ein starkes Signal bestehen, passt dann in einer Schleife eine logistische Regression der Daten gegen einen Satz von Zufallszahlen an und akkumuliert die p-Werte der Zufallsprädiktoren; Die Ergebnisse werden als Histogramm der p-Werte, der tatsächlichen p-Wert-Quantile für Konfidenzniveaus von 1% und 5% und der tatsächlichen Falsch-Positiv-Rate entsprechend diesen Konfidenzniveaus angezeigt.

Ich denke, ein Grund für die unerwarteten Ergebnisse ist, dass ein zufälliger Prädiktor im Allgemeinen eine gewisse Korrelation mit dem wahren Signal aufweist und dies hauptsächlich für die Ergebnisse verantwortlich ist. Wenn Sie jedoch einstellen orthogonalPredictor, besteht TRUEkeine Korrelation zwischen den zufälligen Prädiktoren und dem tatsächlichen Signal, aber das Problem liegt immer noch auf einem reduzierten Niveau. Meine beste Erklärung dafür ist, dass das Modell falsch spezifiziert ist und p-Werte ohnehin verdächtig sind, da sich das wahre Signal nirgendwo im angepassten Modell befindet. Aber dies ist ein Catch-22 - wer hat jemals eine Reihe von Prädiktoren zur Verfügung, die genau zu den Daten passen? Hier sind einige zusätzliche Fragen:

1. Was ist die genaue Nullhypothese für die logistische Regression? Ist es so, dass die Daten rein binomiales Rauschen um einen konstanten Pegel sind (dh es gibt kein echtes Signal)? Wenn Sie sd im Code auf 0 setzen, gibt es kein Signal und das Histogramm sieht flach aus.

2. Die implizite Nullhypothese im Code lautet, dass der Prädiktor nicht mehr Erklärungskraft hat als eine Reihe von Zufallszahlen. Es wird getestet, indem beispielsweise das empirische 5% -Quantil für den vom Code angezeigten p-Wert verwendet wird. Gibt es eine bessere Möglichkeit, diese Hypothese zu testen, oder zumindest eine, die nicht so numerisch intensiv ist?

------ Zusätzliche Information

Dieser Code ahmt das folgende Problem nach: Historische Ausfallraten zeigen signifikante zeitliche Schwankungen (Signale), die durch Konjunkturzyklen bedingt sind. Die tatsächlichen Standardzählungen zu einem bestimmten Zeitpunkt sind binomial um diese Standardwahrscheinlichkeiten. Ich habe versucht, erklärende Variablen für das Signal zu finden, als ich den p-Werten misstrauisch wurde. In diesem Test wird das Signal über die Zeit zufällig geordnet, anstatt die Konjunkturzyklen anzuzeigen, aber das sollte für die logistische Regression keine Rolle spielen. Es gibt also keine Überdispersion, das Signal soll eigentlich ein Signal sein.

Hier gibt es mehrere Probleme. Insbesondere scheint es einige Verwirrungen darüber zu geben, wie eine logistische Standardregression simuliert werden kann. Kurz, Sie fügen nicht hinzu noise around... the probability "signal". Aufgrund der Art und Weise, wie Sie dies getan haben, gibt es eine enorme Variabilität in den resultierenden 'Binomial'-Daten (- esque), weit mehr als es sein sollte. Hier sind die Wahrscheinlichkeiten in Ihrem Datensatz:

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

Wenn diese .4+ Beobachtungen auf der einen oder anderen Seite landen, fungieren sie als "Ausreißer" (sie sind nicht wirklich) und verursachen einen Typ-1-Fehler in einem Modell, das die Tatsache, dass diese Daten nicht berücksichtigt, nicht berücksichtigt sind nicht wirklich binomisch. Hier ist eine Version, die einen einfachen Hack verwendet, damit das Modell eine Überdispersion erkennen und berücksichtigen kann:

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

Dies ist die Modellzusammenfassung der letzten Iteration. Es ist zu beachten, dass die Dispersion auf schätzungsweise geschätzt wird$\approx 12\times$ was es für ein echtes Binomial sein sollte:

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

Hier ist eine andere Version, in der ich das gleiche Modell wie Sie anpasse, aber nur die Daten ohne das zusätzliche Rauschen um das Signal generiere. (Beachten Sie, dass Code, der ansonsten derselbe ist, der Kürze halber weggelassen wird.)

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

Die Beziehung des Codes zum experimentellen Ziel ist verwirrend. Erwarten Sie einen signifikanten Prädiktor für eine Auswahl von orthogonalPredictor oder SD? Ich nicht.

Nach meiner Interpretation scheint das Experiment nicht mit dem übereinzustimmen, was wir zu testen versuchen. Wenn das Rauschen erzeugt wird, wird es implizit wiederholt (dh nicht zufällig) an einzelne Beobachtungen angehängt, was ein Signal für die Regression liefert.

Folgendes war meiner Meinung nach beabsichtigt:

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

Die kritischen Änderungen sind:

• Erweitern des Datenrahmens auf das Beobachtungsformat anstelle eines komprimierten Formats (flips_expanded)
• Experimentieren Sie auch mit einem korrelierten Prädiktor

Für keine Korrelation zwischen y_actual und unserem Prädiktor y:

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

Und eine ziemlich starke Korrelation schaffen:

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