Rzadkie uprzedzenie regresji logistycznej zdarzenia: jak symulować niedoszacowane wartości p na minimalnym przykładzie?

CrossValidated ma kilka pytań na temat tego, kiedy i jak zastosować korektę błędu rzadkich zdarzeń autorstwa Kinga i Zenga (2001) . Szukam czegoś innego: minimalnej demonstracji opartej na symulacji, że istnieje uprzedzenie.

W szczególności państwo King i Zeng

„... w danych dotyczących rzadkich zdarzeń tendencje w prawdopodobieństwach mogą mieć istotne znaczenie przy wielkościach próbek w tysiącach i są w przewidywalnym kierunku: oszacowane prawdopodobieństwa zdarzeń są zbyt małe”.

Oto moja próba symulacji takiego obciążenia w R:

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

Kiedy to uruchamiam, mam tendencję do uzyskiwania bardzo małych wyników Z, a histogram oszacowań jest bardzo bliski wyśrodkowaniu ponad prawdą p = 0,01.

czego mi brakuje? Czy to, że moja symulacja nie jest wystarczająco duża, pokazuje prawdziwe (i ewidentnie bardzo małe) odchylenie? Czy odchylenie wymaga uwzględnienia jakiegoś współzmiennego (więcej niż przecięcia)?

Aktualizacja 1: King i Zeng zawierają przybliżone przybliżenie błędu w równaniu 12 dokumentu. Zwracając uwagę na mianownik, drastycznie zmniejszyłem się i ponownie przeprowadziłem symulację, ale nadal nie widać błędu w szacowanych prawdopodobieństwach zdarzeń. (Użyłem tego tylko jako inspiracji. Pamiętaj, że moje pytanie powyżej dotyczy szacunkowych prawdopodobieństw zdarzeń, a nie .)β$\beta_0$NN5$\hat \beta_0$

Aktualizacja 2: Zgodnie z sugestią zawartą w komentarzach uwzględniłem zmienną niezależną w regresji, co prowadzi do równoważnych wyników:

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

Objaśnienie: Użyłem psiebie jako zmiennej niezależnej, gdzie pjest wektorem z powtórzeniami małej wartości (0,01) i większej wartości (0,2). Ostatecznie simprzechowuje tylko oszacowane prawdopodobieństwa odpowiadające i nie ma żadnych oznak błędu.$p = 0.01$

Aktualizacja 3 (5 maja 2016 r.): To nie zmienia zauważalnie wyników, ale moja nowa funkcja wewnętrznej symulacji to

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

Objaśnienie: MLE, gdy y jest identycznie zero, nie istnieje ( dzięki komentarzom tutaj dla przypomnienia ). R nie rzuca ostrzeżenia, ponieważ jego „ pozytywna tolerancja zbieżności ” faktycznie się spełnia. Mówiąc bardziej swobodnie, MLE istnieje i ma minus nieskończoności, co odpowiada ; stąd moja aktualizacja funkcji. Jedyną inną spójną rzeczą, jaką mogę wymyślić, jest odrzucenie tych przebiegów symulacji, w których y wynosi identycznie zero, ale to wyraźnie doprowadziłoby do wyników jeszcze bardziej sprzecznych z początkowym twierdzeniem, że „szacowane prawdopodobieństwo zdarzenia jest zbyt małe”.$p=0$

To interesujące pytanie - wykonałem kilka symulacji, które zamieszczam poniżej w nadziei, że zachęci to do dalszej dyskusji.

Przede wszystkim kilka ogólnych uwag:

• Artykuł, który cytujesz, dotyczy stronniczości rzadkich zdarzeń. To, co wcześniej nie było dla mnie jasne (również w odniesieniu do komentarzy, które zostały przedstawione powyżej), to czy jest coś specjalnego w przypadkach, w których masz 10/10000 w przeciwieństwie do obserwacji 10/30. Jednak po kilku symulacjach zgodziłbym się z tym.

• Problemem, o którym myślałem (często się z tym spotykałem, a ostatnio był w nim artykuł w Methods in Ecology and Evolution, ale nie mogłem znaleźć odniesienia), że można uzyskać zdegenerowane przypadki z GLM w małych danych sytuacje, w których MLE jest FAAAR z dala od prawdy, a nawet w nieskończoności (- przypuszczam, że z powodu połączenia nieliniowego). Nie jest dla mnie jasne, jak należy traktować te przypadki w oszacowaniu błędu, ale z moich symulacji powiedziałbym, że wydają się kluczowe dla błędu rzadkich zdarzeń. Moją intuicją byłoby ich usunięcie, ale wtedy nie jest całkiem jasne, jak daleko muszą być usunięte. Może coś, o czym należy pamiętać przy korekcji błędów.

• Te zdegenerowane przypadki również wydają się powodować problemy numeryczne (dlatego zwiększyłem maksimum w funkcji glm, ale można również pomyśleć o zwiększeniu epsilon, aby upewnić się, że rzeczywiście zgłasza się prawdziwy MLE).

W każdym razie, oto kod, który oblicza różnicę między szacunkami a prawdą dla przechwytywania, nachylenia i prognoz w regresji logistycznej, najpierw dla sytuacji o małej wielkości próby / umiarkowanej częstości:

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)

Wynikowe odchylenie i standardowe błędy dla przechwytywania, nachylenia i prognozowania wynoszą

-0.120429315  0.296453122 -0.001619793
 0.016105833  0.032835468  0.002040664

Doszedłbym do wniosku, że istnieją dość dobre dowody na niewielkie ujemne odchylenie w punkcie przecięcia i lekkie dodatnie odchylenie na zboczu, chociaż spojrzenie na wykreślone wyniki pokazuje, że odchylenie jest niewielkie w porównaniu z wariancją wartości szacunkowych.

Jeśli ustawiam parametry na rzadkie zdarzenia

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

Dostaję większe odchylenie dla przechwytywania, ale nadal NIE ma prognozy

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

Na histogramie wartości szacunkowych widzimy zjawisko oszacowania parametrów zdegenerowanych (jeśli powinniśmy je tak nazwać)

Usuńmy wszystkie wiersze, dla których oszacowania przechwytywania wynoszą <20

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

Błąd systematyczny maleje, a na rysunkach wszystko staje się bardziej wyraźne - oszacowania parametrów wyraźnie nie są normalnie rozłożone. Zastanawiam się, że to oznacza ważność zgłaszanych elementów CI.

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

Doszedłbym do wniosku, że stronniczość rzadkich zdarzeń na przechwytywaniu jest spowodowana przez same rzadkie zdarzenia, a mianowicie te rzadkie, bardzo małe oszacowania. Nie jestem pewien, czy chcemy je usunąć, czy nie, nie wiem, co to będzie granica.

Należy jednak zauważyć, że tak czy inaczej, wydaje się, że nie ma stronniczości w prognozach w skali odpowiedzi - funkcja link po prostu absorbuje te niezwykle małe wartości.

Uprzedzenie rzadkich zdarzeń występuje tylko wtedy, gdy występują regresory. Nie wystąpi w modelu tylko przechwytującym, takim jak ten symulowany tutaj. Zobacz ten post, aby uzyskać szczegółowe informacje: http://statistichorizons.com/linear-vs-logistic#comment-276108

Wydaje się, że ryc. 7 w artykule najbardziej bezpośrednio odnosi się do błędu w prognozach. Nie do końca rozumiem ten rysunek (a konkretnie interpretacja „szacunkowe prawdopodobieństwa zdarzeń są zbyt małe” wydaje się nadmiernym uproszczeniem), ale udało mi się odtworzyć coś podobnego na podstawie ich zwięzłego opisu ich symulacji w Rozdziale 6.1:

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

Pierwszy wykres to moja replikacja ich ryciny 7, z dodatkiem przerywanych krzywych reprezentujących pełny zakres wyników w 10 próbach.

Zgodnie z artykułem, xtutaj jest zmienna predykcyjna w regresji zaczerpnięta ze standardowej normy. Tak więc, jak pokazano na drugim wykresie, względna częstotliwość obserwacji dla x > 3(gdzie najbardziej widoczne odchylenie występuje na pierwszym wykresie) jest coraz mniejsza.

Trzeci wykres pokazuje „prawdziwe” prawdopodobieństwa symulacji w procesie generowania w funkcji x. Wydaje się, że największe uprzedzenie występuje tam, gdzie xjest rzadkie lub nie istnieje.

Podsumowując, sugerują one, że OP całkowicie błędnie zinterpretował główne twierdzenie artykułu, myląc „rzadkie zdarzenie” (tj. x > 3) Z „zdarzeniem, dla którego P(Y = 1)jest bardzo małe”. Prawdopodobnie artykuł dotyczy tego pierwszego zamiast drugiego.