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 .)β 0N
N
5
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 p
siebie jako zmiennej niezależnej, gdzie p
jest wektorem z powtórzeniami małej wartości (0,01) i większej wartości (0,2). Ostatecznie sim
przechowuje tylko oszacowane prawdopodobieństwa odpowiadające i nie ma żadnych oznak błędu.
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”.