Logistische Regressionsverzerrung bei seltenen Ereignissen: Wie kann man die unterschätzten p mit einem minimalen Beispiel simulieren?

CrossValidated hat verschiedene Fragen, wann und wie die Selten-Ereignis-Bias-Korrektur von King und Zeng (2001) angewendet werden soll. . Ich suche etwas anderes: eine minimale simulationsbasierte Demonstration, dass der Bias existiert.

Insbesondere König und Zeng Zustand

"... in Daten zu seltenen Ereignissen können die Wahrscheinlichkeiten bei Stichprobengrößen zu Tausenden erheblich sein und in eine vorhersagbare Richtung weisen: Geschätzte Ereigniswahrscheinlichkeiten sind zu klein."

Hier ist mein Versuch, eine solche Verzerrung in R zu simulieren:

# FUNCTIONS
do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)  
    # Extract the fitted probability.
    #    If p is constant, glm does y ~ 1, the intercept-only model.
    #    If p is not constant, assume its smallest value is p[1]:
    glm(y ~ p, family = 'binomial')$fitted[1]
}
mean.of.K.estimates = function(p, K){
    mean(replicate(K, do.one.sim(p) ))
}

# MONTE CARLO
N = 100
p = rep(0.01, N)
reps = 100
# The following line may take about 30 seconds
sim = replicate(reps, mean.of.K.estimates(p, K=100))
# Z-score:
abs(p[1]-mean(sim))/(sd(sim)/sqrt(reps))
# Distribution of average probability estimates:
hist(sim)

Wenn ich dies durchführe, neige ich dazu, sehr kleine Z-Scores zu erhalten, und das Histogramm der Schätzungen ist sehr nahe daran, über der Wahrheit p = 0,01 zu zentrieren.

Was vermisse ich? Ist es so, dass meine Simulation nicht groß genug ist, um die wahre (und offensichtlich sehr kleine) Tendenz aufzuzeigen? Benötigt die Verzerrung eine Art Kovariate (mehr als der Achsenabschnitt), um einbezogen zu werden?

Update 1: King und Zeng geben in Gleichung 12 ihrer Arbeit eine grobe Näherung für die Verzerrung von an. Unter Hinweis auf die im Nenner, ich drastisch reduziert werden und wieder lief die Simulation, aber noch keine Verzerrung der geschätzten Ereigniswahrscheinlichkeiten ist evident. (Ich habe dies nur als Inspiration verwendet. Beachten Sie, dass es sich bei meiner obigen Frage um geschätzte Ereigniswahrscheinlichkeiten handelt, nicht um .)$\beta_0$NN5$\hat \beta_0$

Update 2: Aufgrund eines Vorschlags in den Kommentaren habe ich eine unabhängige Variable in die Regression aufgenommen, was zu äquivalenten Ergebnissen führte:

p.small = 0.01
p.large = 0.2
p = c(rep(p.small, round(N/2) ), rep(p.large, N- round(N/2) ) )
sim = replicate(reps, mean.of.K.estimates(p, K=100))

Erklärung: Ich habe mich pselbst als unabhängige Variable verwendet, wobei pes sich um einen Vektor mit Wiederholungen eines kleinen Werts (0,01) und eines größeren Werts (0,2) handelt. simSpeichert am Ende nur die geschätzten Wahrscheinlichkeiten, die$p = 0.01$ und es gibt kein Anzeichen einer Verzerrung.

Update 3 (5. Mai 2016): Das ändert nichts an den Ergebnissen, aber meine neue innere Simulationsfunktion ist

do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(0)
    }else{
        # Extract the fitted probability.
        #    If p is constant, glm does y ~ 1, the intercept only model.
        #    If p is not constant, assume its smallest value is p[1]:
        return(glm(y ~ p, family = 'binomial')$fitted[1])
    }
}

Erläuterung: Die MLE existiert nicht, wenn y identisch Null ist ( dank Kommentaren hier für die Erinnerung ). R gibt keine Warnung aus, weil seine " positive Konvergenztoleranz " tatsächlich erfüllt wird. Im Allgemeinen existiert die MLE und ist minus unendlich, was . daher mein Funktionsupdate. Die einzige andere zusammenhängende Sache, an die ich denken kann, ist das Verwerfen der Läufe der Simulation, bei denen y identisch Null ist. Dies würde jedoch eindeutig zu Ergebnissen führen, die der ursprünglichen Behauptung, dass "geschätzte Ereigniswahrscheinlichkeiten zu klein sind", noch mehr entgegenwirken.$p=0$

Dies ist eine interessante Frage - ich habe ein paar Simulationen durchgeführt, die ich unten poste, in der Hoffnung, dass dies die weitere Diskussion anregt.

Zunächst einige allgemeine Bemerkungen:

• In dem von Ihnen zitierten Artikel geht es um die Selten-Ereignis-Voreingenommenheit. Was mir vorher nicht klar war (auch in Bezug auf die Kommentare, die oben gemacht wurden), ist, ob es etwas Besonderes an Fällen gibt, in denen Sie 10/10000 gegenüber 10/30 Beobachtungen haben. Nach einigen Simulationen würde ich jedoch zustimmen, dass dies der Fall ist.

• Ein Problem, an das ich gedacht habe (ich bin dies oft angetroffen, und es gab kürzlich einen Artikel in Methoden in Ökologie und Evolution dazu, ich konnte jedoch keine Referenz finden), ist, dass man mit GLMs in kleinen Daten degenerierte Fälle bekommen kann Situationen, in denen die MLE FAAAR von der Wahrheit entfernt ist oder sogar bei - / + unendlich (aufgrund der nichtlinearen Verknüpfung, nehme ich an). Mir ist nicht klar, wie man diese Fälle in der Bias-Schätzung behandeln soll, aber aus meinen Simulationen würde ich sagen, dass sie für die Selten-Ereignis-Bias entscheidend sind. Meine Intuition wäre, sie zu entfernen, aber dann ist nicht ganz klar, wie weit sie entfernt werden müssen. Vielleicht etwas zu beachten für die Bias-Korrektur.

• Auch scheinen diese entarteten Fälle anfällig für numerische Probleme zu sein (ich habe daher das Maximum in der glm-Funktion erhöht, aber man könnte auch darüber nachdenken, das Epsilon zu erhöhen, um sicherzustellen, dass man das wahre MLE tatsächlich meldet).

Wie auch immer, hier ein Code, der die Differenz zwischen Schätzungen und Wahrheiten für Achsenabschnitt, Steigung und Vorhersagen in einer logistischen Regression berechnet, zunächst für eine Situation mit geringer Stichprobengröße / mittlerer Inzidenz:

set.seed(123)
replicates = 1000
N= 40
slope = 2 # slope (linear scale)
intercept = - 1 # intercept (linear scale)

bias <- matrix(NA, nrow = replicates, ncol = 3)
incidencePredBias <- rep(NA, replicates)

for (i in 1:replicates){
  pred = runif(N,min=-1,max=1) 
  linearResponse = intercept + slope*pred
  data = rbinom(N, 1, plogis(linearResponse))  
  fit <- glm(data ~ pred, family = 'binomial', control = list(maxit = 300))
  bias[i,1:2] = fit$coefficients - c(intercept, slope)
  bias[i,3] = mean(predict(fit,type = "response")) - mean(plogis(linearResponse))
}

par(mfrow = c(1,3))
text = c("Bias intercept", "Bias slope", "Bias prediction")

for (i in 1:3){
  hist(bias[,i], breaks = 100, main = text[i])
  abline(v=mean(bias[,i]), col = "red", lwd = 3)  
}

apply(bias, 2, mean)
apply(bias, 2, sd) / sqrt(replicates)

Die resultierenden Bias- und Standardfehler für Achsenabschnitt, Steigung und Vorhersage sind

-0.120429315  0.296453122 -0.001619793
 0.016105833  0.032835468  0.002040664

Ich würde daraus schließen, dass es ziemlich gute Hinweise auf eine leichte negative Abweichung im Achsenabschnitt und eine leichte positive Abweichung in der Steigung gibt, obwohl ein Blick auf die grafischen Ergebnisse zeigt, dass die Abweichung im Vergleich zur Varianz der geschätzten Werte gering ist.

Wenn ich die Parameter auf eine Situation mit seltenen Ereignissen einstelle

N= 4000
slope = 2 # slope (linear scale)
intercept = - 10 # intercept (linear scale)

Ich habe eine größere Neigung für den Abschnitt, aber immer noch KEINE für die Vorhersage

   -1.716144e+01  4.271145e-01 -3.793141e-06
    5.039331e-01  4.806615e-01  4.356062e-06

Im Histogramm der geschätzten Werte sehen wir das Phänomen der entarteten Parameterschätzungen (wenn wir sie so nennen sollten)

Entfernen wir alle Zeilen, für die die Intercept-Schätzungen <20 sind

apply(bias[bias[,1] > -20,], 2, mean)
apply(bias[bias[,1] > -20,], 2, sd) / sqrt(length(bias[,1] > -10))

Die Verzerrung nimmt ab und die Zahlen verdeutlichen die Situation - Parameterschätzungen sind eindeutig nicht normalverteilt. Ich frage mich, dass dies für die Gültigkeit der gemeldeten CIs bedeutet.

-0.6694874106  1.9740437782  0.0002079945
1.329322e-01 1.619451e-01 3.242677e-06

Ich würde den Schluss ziehen, dass die Tendenz zu seltenen Ereignissen auf dem Achsenabschnitt von seltenen Ereignissen selbst getrieben wird, nämlich diesen seltenen, extrem kleinen Schätzungen. Nicht sicher, ob wir sie entfernen möchten oder nicht, nicht sicher, wie hoch der Cutoff sein würde.

Es ist jedoch wichtig zu beachten, dass in beiden Fällen keine Verzerrung der Vorhersagen auf der Antwortskala zu bestehen scheint - die Link-Funktion absorbiert einfach diese extrem kleinen Werte.

Seltene Vorurteile treten nur auf, wenn es Regressoren gibt. Es wird nicht in einem Intercept-Only-Modell wie dem hier simulierten auftreten. Weitere Informationen finden Sie in diesem Beitrag: http://statisticalhorizons.com/linear-vs-logistic#comment-276108

Abbildung 7 in der Veröffentlichung scheint die Frage der Verzerrung in den Vorhersagen direkt zu behandeln. Ich verstehe die Abbildung nicht vollständig (insbesondere scheint die Interpretation "Geschätzte Ereigniswahrscheinlichkeiten sind zu klein" eine zu starke Vereinfachung zu sein), aber ich habe es geschafft, basierend auf der knappen Beschreibung ihrer Simulation in Abschnitt 6.1 etwas Ähnliches zu reproduzieren:

n_grid = 40
x_grid = seq(0, 7, length.out = n_grid)
beta0 = -6
beta1 = 1

inverse_logit = function(x) 1/(1 + exp(-x))

do.one.sim = function(){
    N = 5000
    x = rnorm(N)
    p = inverse_logit(beta0 + beta1*x)
    # Draw fake data based on probabilities p
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(rep(0, n_grid))
    }else{
        # Extract the error
        mod = glm(y ~ x, family = 'binomial')
        truth = inverse_logit(beta0 + beta1*x_grid)
        pred = predict(mod, newdata = data.frame(x = x_grid),
            type = 'response')
        return(pred - truth)
    }
}
mean.of.K.estimates = function(K){
    rowMeans(replicate(K, do.one.sim()))
}

set.seed(1)
bias = replicate(10, mean.of.K.estimates(100))
maxes = as.numeric(apply(bias, 1, max))
mins = as.numeric(apply(bias, 1, min))

par(mfrow = c(3, 1), mar = c(4,4,2,2))
plot(x_grid, rowMeans(bias), type = 'l',
    ylim = c(min(bias), max(bias)),
    xlab = 'x', ylab = 'bias')
lines(x_grid, maxes, lty = 2)
lines(x_grid, mins, lty = 2)
plot(x_grid, dnorm(x_grid), type = 'l',
    xlab = 'x', ylab = 'standard normal density')
plot(x_grid, inverse_logit(beta0 + beta1*x_grid),
    xlab = 'x', ylab = 'true simulation P(Y = 1)',
    type = 'l')

Der erste Plot ist meine Nachbildung von Abbildung 7, wobei gestrichelte Kurven hinzugefügt wurden, die den gesamten Bereich der Ergebnisse über 10 Versuche darstellen.

Nach dem Papier ist xhier eine Prädiktorvariable in der Regression aus einer Standardnormalen. Wie in der zweiten Auftragung dargestellt, ist die relative Häufigkeit der Beobachtungen x > 3(bei der die sichtbarste Verzerrung in der ersten Auftragung auftritt) daher sehr gering.

Das dritte Diagramm zeigt die "wahren" Simulationswahrscheinlichkeiten im Erzeugungsprozess als Funktion von x. Es scheint, dass die größte Verzerrung dort auftritt, wo sie xselten oder nicht vorhanden ist.

Zusammengenommen deuten diese darauf hin, dass das OP die zentrale Behauptung des Papiers völlig falsch interpretierte, indem es "seltenes Ereignis" (dh x > 3) mit "Ereignis, für das P(Y = 1)es sehr klein ist" verwechselte . Vermutlich betrifft das Papier das erstere anstelle des letzteren.