Viés de regressão logística de eventos raros: como simular os p's subestimados com um exemplo mínimo?

O CrossValidated tem várias perguntas sobre quando e como aplicar a correção de viés de evento raro por King e Zeng (2001) . Estou procurando algo diferente: uma demonstração mínima baseada em simulação de que o viés existe.

Em particular, o rei e o estado de Zeng

"... em dados de eventos raros, os desvios nas probabilidades podem ser substancialmente significativos com tamanhos de amostra na casa dos milhares e estão em uma direção previsível: as probabilidades estimadas de eventos são muito pequenas."

Aqui está minha tentativa de simular esse viés no 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)

Quando executo isso, tendem a obter escores z muito pequenos, e o histograma de estimativas está muito próximo da verdade p = 0,01.

o que estou perdendo? Será que minha simulação não é grande o suficiente para mostrar o verdadeiro viés (e evidentemente muito pequeno)? O viés exige que algum tipo de covariável (mais que a interceptação) seja incluído?

Atualização 1: King e Zeng incluem uma aproximação aproximada para o viés de na equação 12 de seu artigo. Observando o denominador, reduzi drasticamente a ser e refiz a simulação, mas ainda não é evidente nenhum viés nas probabilidades estimadas do evento. (Eu usei isso apenas como inspiração. Note-se que a minha pergunta acima é sobre probabilidades de eventos estimados, não β 0 ).$\beta_0$NN5$\hat \beta_0$

Atualização 2: Seguindo uma sugestão nos comentários, incluí uma variável independente na regressão, levando a resultados equivalentes:

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

Explicação: Eu me usei pcomo a variável independente, onde pé um vetor com repetições de um valor pequeno (0,01) e um valor maior (0,2). No final, simarmazena apenas as probabilidades estimadas correspondentes a e não há sinal de viés.$p = 0.01$

Atualização 3 (5 de maio de 2016): Isso não altera visivelmente os resultados, mas minha nova função de simulação interna é

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

$p=0$

Essa é uma pergunta interessante - fiz algumas simulações que publico abaixo, na esperança de que isso estimule uma discussão mais aprofundada.

Primeiro de tudo, alguns comentários gerais:

• O artigo que você cita é sobre viés de eventos raros. O que não estava claro para mim antes (também com relação aos comentários que foram feitos acima) é se há algo de especial nos casos em que você tem 10/10000 em oposição a 10/30 observações. No entanto, após algumas simulações, eu concordo que existe.

• Um problema que eu tinha em mente (eu encontrei isso com frequência, e recentemente houve um artigo em Métodos em Ecologia e Evolução sobre isso, mas não consegui encontrar a referência) é que você pode obter casos degenerados com GLMs em dados pequenos situações em que o MLE está FAAAR longe da verdade, ou mesmo no - / + infinito (devido ao link não-linear, suponho). Não está claro para mim como se deve tratar esses casos na estimativa de viés, mas pelas minhas simulações eu diria que eles parecem essenciais para o viés de evento raro. Minha intuição seria removê-los, mas não está claro até que ponto eles precisam ser removidos. Talvez algo a ter em mente para a correção de preconceitos.

• Além disso, esses casos degenerados parecem propensos a causar problemas numéricos (portanto, aumentei o maxit na função glm, mas também se poderia pensar em aumentar o epsilon para garantir que o relatório realmente seja verdadeiro).

De qualquer forma, aqui está um código que calcula a diferença entre estimativas e verdade para interceptação, inclinação e previsões em uma regressão logística, primeiro para uma amostra pequena / situação de incidência moderada:

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)

O viés resultante e os erros padrão para interceptação, inclinação e previsão são

-0.120429315  0.296453122 -0.001619793
 0.016105833  0.032835468  0.002040664

Concluo que há boas evidências para um leve viés negativo na interceptação e um leve viés positivo na inclinação, embora uma olhada nos resultados plotados mostre que o viés é pequeno comparado à variação dos valores estimados.

Se estou definindo os parâmetros para uma situação de evento raro

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

Estou recebendo um viés maior pela interceptação, mas ainda NENHUM na previsão

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

No histograma dos valores estimados, vemos o fenômeno de estimativas degeneradas de parâmetros (se deveríamos chamá-las assim)

Vamos remover todas as linhas para as quais as estimativas de interceptação são <20

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

O viés diminui e as coisas ficam um pouco mais claras nas figuras - as estimativas de parâmetros claramente não são normalmente distribuídas. Eu me pergunto que isso significa para a validade dos ICs que são relatados.

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

Concluo que o viés de evento raro na interceptação é impulsionado pelos próprios eventos raros, ou seja, aquelas estimativas raras e extremamente pequenas. Não tenho certeza se queremos removê-los ou não, não tenho certeza qual seria o ponto de corte.

Uma coisa importante a ser observada, porém, é que, de qualquer forma, parece não haver viés nas previsões na escala de resposta - a função de link simplesmente absorve esses valores extremamente pequenos.

O viés de eventos raros ocorre apenas quando há regressores. Não ocorrerá em um modelo somente de interceptação como o simulado aqui. Consulte esta publicação para obter detalhes: http://statisticalhorizons.com/linear-vs-logistic#comment-276108

A Figura 7 no artigo parece abordar mais diretamente a questão do viés nas previsões. Não entendo completamente a figura (especificamente, a interpretação "as probabilidades estimadas de eventos são muito pequenas" parece uma simplificação excessiva), mas consegui reproduzir algo semelhante a ele com base na descrição concisa da simulação na Seção 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')

O primeiro gráfico é minha replicação da figura 7, com a adição de curvas tracejadas que representam toda a gama de resultados em 10 tentativas.

De acordo com o artigo, xaqui está uma variável preditora na regressão extraída de uma normal padrão. Assim, como ilustrado no segundo gráfico, a frequência relativa de observações para x > 3(onde o viés mais visível ocorre no primeiro gráfico) é cada vez menor.

O terceiro gráfico mostra as probabilidades "verdadeiras" de simulação no processo de geração em função de x. Parece que o maior viés ocorre onde xé raro ou inexistente.

Tomados em conjunto, eles sugerem que o OP mal interpretou completamente a alegação central do artigo, confundindo "evento raro" (ie x > 3) com "evento para o qual P(Y = 1)é muito pequeno". Presumivelmente, o artigo diz respeito ao primeiro, e não ao segundo.