Por que o intervalo credível bayesiano nessa regressão polinomial é enviesado, enquanto o intervalo de confiança está correto?

9

Considere o gráfico abaixo no qual simulei os dados da seguinte maneira. Observamos um resultado binário para o qual a verdadeira probabilidade de ser 1 é indicada pela linha preta. A relação funcional entre uma covariável e é 3 polinomial ordem com ligação logístico (por isso é não linear num modo duplo-). x p ( y o b s = 1 | x )yobsxp(yobs=1|x)

A linha verde é o ajuste de regressão logística GLM, onde é introduzido como polinômio de terceira ordem. As linhas verdes tracejadas são os intervalos de confiança de 95% em torno da previsão , onde os coeficientes de regressão ajustados. Eu usei e para isso.p ( y o b s = 1 | x , β ) βxp(yobs=1|x,β^)β^R glmpredict.glm

Da mesma forma, a linha do pruple é a média do posterior com intervalo credível de 95% para de um modelo de regressão logística bayesiana usando um uniforme anterior. Eu usei o pacote com a função para isso (a configuração fornece o uniforme não informativo anterior).p(yobs=1|x,β)MCMCpackMCMClogitB0=0

Os pontos vermelhos indicam observações no conjunto de dados para o qual , os pontos pretos são observações com . Observe que, como comum na classificação / análise discreta, mas não é observado.y o b s = 0 y p ( y o b s = 1 | x )yobs=1yobs=0yp(yobs=1|x)

insira a descrição da imagem aqui

Várias coisas podem ser vistas:

  1. Simulei de propósito que é escasso na mão esquerda. Quero que o intervalo de confiança e credibilidade se amplie aqui devido à falta de informações (observações).x
  2. Ambas as previsões são tendenciosas para cima à esquerda. Esse viés é causado pelos quatro pontos vermelhos que observações, o que sugere erroneamente que a verdadeira forma funcional aqui. O algoritmo possui informações insuficientes para concluir que a verdadeira forma funcional está dobrada para baixo.yobs=1
  3. O intervalo de confiança aumenta conforme o esperado, enquanto o intervalo credível não . De fato, o intervalo de confiança inclui o espaço completo dos parâmetros, como deveria devido à falta de informações.

Parece que o intervalo credível está errado / muito otimista aqui para uma parte de . É realmente um comportamento indesejável que o intervalo credível seja reduzido quando as informações ficam escassas ou totalmente ausentes. Geralmente não é assim que um intervalo confiável reage. Alguém pode explicar:x

  1. Quais são as razões para isso?
  2. Que medidas posso tomar para obter um intervalo melhor credível? (ou seja, um que inclua pelo menos a verdadeira forma funcional, ou melhor, alcance o intervalo de confiança)

O código para obter intervalos de previsão no gráfico é impresso aqui:

fit <- glm(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
x_pred <- seq(0, 1, by=0.01)
pred <- predict(fit, newdata = data.frame(x=x_pred), se.fit = T)
plot(plogis(pred$fit), type='l')
matlines(plogis(pred$fit + pred$se.fit %o% c(-1.96,1.96)), type='l', col='black', lty=2)


library(MCMCpack)
mcmcfit <- MCMClogit(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
gibbs_samps <- as.mcmc(mcmcfit)
x_pred_dm <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=x_pred))
gibbs_preds <- apply(gibbs_samps, 1, `%*%`, t(x_pred_dm))
gibbs_pis <- plogis(apply(gibbs_preds, 1, quantile, c(0.025, 0.975)))
matlines(t(gibbs_pis), col='red', lty=2)

Acesso a dados : https://pastebin.com/1H2iXiew, obrigado @DeltaIV e @AdamO

tomka
fonte
Se alguém puder me explicar como compartilhar uma tabela com os dados, eu posso fazê-lo.
Tomka
Você pode usar dputo quadro de dados que contém os dados e incluir a dputsaída como código na sua postagem.
DeltaIV
11
@ tomka oh eu vejo. Eu não sou daltônico, mas é muito difícil para mim ver a diferença verde / azul!
AdamO 29/09
11
@AdamO espero que este é melhor
Tomka

Respostas:

6

Para um modelo frequencista, a variância dos amplia previsão em proporção ao quadrado da distância do centro de gravidade . Seu método de calcular intervalos de previsão para um GLM Bayesiano usa quantis empíricos com base na curva de probabilidade ajustada, mas não leva em consideração a alavancagem deXXX

Um GLM binomial frequentista não é diferente de um GLM com link de identidade, exceto que a variação é proporcional à média.

Observe que qualquer representação polinomial das probabilidades do logit leva a previsões de risco que convergem para 0 como e 1 como ou vice-versa, dependendo do sinal do termo de ordem polinomial mais alta .X XX

Para previsões freqüentes, o aumento proporcional do desvio ao quadrado (alavancagem) na variação das previsões domina essa tendência. É por isso que a taxa de convergência para intervalos de previsão aproximadamente igual a [0, 1] é mais rápida que a convergência polinomial de logit de terceira ordem para probabilidades de 0 ou 1 de maneira singular.

O mesmo não ocorre com os quantis ajustados posteriormente Bayesianos. Não há uso explícito de desvio ao quadrado; portanto, contamos simplesmente com a proporção de tendências dominantes de 0 ou 1 para construir intervalos de previsão de longo prazo.

Isto é feito aparente extrapolando muito longe para os extremos de .X

Usando o código que forneci acima, obtemos:

> x_pred_dom <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=c(1000)))
> gibbs_preds <- plogis(apply(gibbs_samps[1000:10000, ], 1, `%*%`, t(x_pred_dom))) # a bunch of 0/1s basically past machine precision
> prop.table(table(gibbs_preds))
gibbs_preds
         0          1 
0.97733585 0.02266415 
> 

Então, 97,75% das vezes, o terceiro termo polinomial era negativo. Isso pode ser verificado nas amostras de Gibbs:

> prop.table(table(gibbs_samps[, 4]< 0))

 FALSE   TRUE 
0.0225 0.9775 

Portanto, a probabilidade prevista converge para 0 quando vai para o infinito. Se inspecionarmos as SEs do modelo bayesiano, descobrimos que a estimativa do terceiro termo polinomial é -185,25 com se 108,81, ou seja, 1,70 DPs de 0; portanto, usando leis de probabilidade normais, ela deve ficar abaixo de 0 95,5% das vezes ( uma previsão muito diferente com base em 10.000 iterações). Apenas outra maneira de entender esse fenômeno.X

Por outro lado, o ajuste freqüente atinge 0,1, conforme o esperado:

freq <- predict(fit, newdata = data.frame(x=1000), se.fit=T)
plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)

dá:

> plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)
     [,1]
[1,]    0
[2,]    1
AdamO
fonte
Não obstante: o modelo bayesiano não está confiante demais em áreas dos dados quais não vê exemplos? Eu sei que posteriores bayesianos ou distribuições preditivas geralmente têm um comportamento muito diferente (isto é, mais parecido com o intervalo de conf.). Suspeito que exista algum impacto do anterior. Se você manipula , especifica a precisão de um anterior normal e pode observar um grande impacto no intervalo credível. xB0MCMClogit
tomka 29/09/17
@ tomka Eu não sei como responder exatamente, pois parece tangencial à pergunta em questão. O mais importante é ressaltar que esses métodos de cálculo de IPs não são realmente comparáveis, especialmente no que se refere à extrapolação. Obviamente, com a inferência bayesiana, se você usa um prior informativo, ganha eficiência quando o prior está certo e perde quando o prior está errado.
AdamO 29/09/17
Só para que você saiba que ainda estou pensando em sua resposta. Eu ainda sinto que é estranho que o posterior não reaja à esparsidade ao aumentar. Acredito que para outros priores um melhor comportamento na região esparsa possa ser alcançado. Não posso definir exatamente isso no momento; Talvez eu aprimore a pergunta com um exemplo em que o intervalo confiável funciona da maneira que eu esperaria, mesmo no caso de extrapolação (estou pensando na regressão bayesiana linear normal, em particular). Quando eu fizer, eu vou deixar você saber.
Tomka