Pourquoi l'intervalle crédible bayésien dans cette régression polynomiale est-il biaisé alors que l'intervalle de confiance est correct?

9

Considérez le graphique ci-dessous dans lequel j'ai simulé des données comme suit. Nous regardons un résultat binaire pour lequel la vraie probabilité d'être 1 est indiquée par la ligne noire. La relation fonctionnelle entre une covariable x et est un polynôme de troisième ordre avec lien logistique (il est donc non linéaire dans un double sens).yobsxp(yobs=1|x)

La ligne verte est l'ajustement de régression logistique GLM où est introduit comme polynôme de troisième ordre. Les lignes vertes en pointillés sont les intervalles de confiance à 95% autour de la prédiction , où les coefficients de régression ajustés. J'ai utilisé et pour ça.p ( y o b s = 1 | x , β ) βxp(yobs=1|x,β^)β^R glmpredict.glm

De même, la ligne du pruple est la moyenne du postérieur avec un intervalle crédible à 95% pour d'un modèle de régression logistique bayésienne utilisant un a priori uniforme. J'ai utilisé le package avec fonction pour cela (le réglage donne la priorité uniforme non informative).p(yobs=1|x,β)MCMCpackMCMClogitB0=0

Les points rouges désignent des observations dans l'ensemble de données pour lesquelles , les points noirs sont des observations avec . Notez que comme courant dans la classification / analyse discrète mais pas est observé.y o b s = 0 y p ( y o b s = 1 | x )yobs=1yobs=0yp(yobs=1|x)

entrez la description de l'image ici

Plusieurs choses peuvent être vues:

  1. J'ai simulé exprès que est rare sur la main gauche. Je veux que la confiance et l'intervalle crédible s'élargissent ici en raison du manque d'informations (observations).x
  2. Les deux prédictions sont biaisées vers le haut à gauche. Ce biais est causé par les quatre points rouges dénotant observations, ce qui suggère à tort que la véritable forme fonctionnelle remonterait ici. L'algorithme ne dispose pas d'informations suffisantes pour conclure que la véritable forme fonctionnelle est courbée vers le bas.yobs=1
  3. L'intervalle de confiance devient plus large comme prévu, alors que l'intervalle crédible ne pas . En fait, l'intervalle de confiance renferme tout l'espace des paramètres, comme il se doit en raison du manque d'informations.

Il semble que l'intervalle crédible soit faux / trop optimiste ici pour une partie de . Il est vraiment indésirable que l'intervalle crédible se rétrécisse lorsque les informations deviennent rares ou totalement absentes. Habituellement, ce n'est pas ainsi que réagit un intervalle crédible. Quelqu'un peut-il expliquer:x

  1. Quelles en sont les raisons?
  2. Quelles mesures puis-je prendre pour arriver à un meilleur intervalle crédible? (c'est-à-dire qui englobe au moins la vraie forme fonctionnelle, ou mieux devient aussi large que l'intervalle de confiance)

Le code pour obtenir les intervalles de prédiction dans le graphique est imprimé ici:

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)

Accès aux données : https://pastebin.com/1H2iXiew merci @DeltaIV et @AdamO

tomka
la source
Si quelqu'un pouvait m'expliquer comment partager un tableau avec les données, je peux le faire.
tomka
Vous pouvez utiliser dputsur la trame de données contenant les données, puis inclure la dputsortie sous forme de code dans votre publication.
DeltaIV
1
@tomka oh je vois. Je ne suis pas daltonien mais c'est très difficile pour moi de voir la différence vert / bleu!
AdamO
1
@AdamO espère que cela va mieux
tomka

Réponses:

6

XX

Un GLM fréquentiste binomial n'est pas différent d'un GLM avec lien d'identité, sauf que la variance est proportionnelle à la moyenne.

XX

Pour la prédiction fréquentiste, l'augmentation proportionnelle au carré de l'écart (effet de levier) de la variance des prédictions domine cette tendance. C'est pourquoi le taux de convergence vers des intervalles de prédiction approximativement égaux à [0, 1] est plus rapide que la convergence logit polynomiale du troisième ordre avec des probabilités de 0 ou 1 singulièrement.

Ce n'est pas le cas pour les quantiles ajustés postérieurs bayésiens. Il n'y a pas d'utilisation explicite de l'écart au carré, nous nous appuyons donc simplement sur la proportion de tendances dominantes 0 ou 1 pour construire des intervalles de prédiction à long terme.

X

En utilisant le code que j'ai fourni ci-dessus, nous obtenons:

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

Donc, 97,75% du temps, le troisième terme polynomial était négatif. Ceci est vérifiable à partir des échantillons Gibbs:

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

 FALSE   TRUE 
0.0225 0.9775 

X

En revanche, l'ajustement fréquentiste souffle jusqu'à 0,1 comme prévu:

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

donne:

> plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)
     [,1]
[1,]    0
[2,]    1
AdamO
la source
xB0MCMClogit
@tomka Je ne sais pas exactement comment répondre à cette question, car cela semble tangentiel à la question posée. Le plus important est de souligner que ces méthodes de calcul des IP ne sont pas vraiment comparables, surtout en ce qui concerne l'extrapolation. Bien sûr, avec l'inférence bayésienne, si vous utilisez un prieur informatif, vous gagnez en efficacité lorsque le prieur a raison et perdez lorsque le prieur est faux.
AdamO
Juste pour vous faire savoir que je pense toujours à votre réponse. J'ai toujours l'impression qu'il est étrange que le postérieur ne réagisse pas à la rareté en s'élargissant. Je crois que pour d'autres prieurs, un meilleur comportement dans la région clairsemée peut être atteint. Je ne peux pas le préciser exactement pour le moment; Je vais peut-être enrichir la question avec un exemple où l'intervalle crédible fonctionne comme je m'y attendais, même en cas d'extrapolation (je pense en particulier à la régression bayésienne linéaire normale). Quand je le ferai, je vous le ferai savoir.
tomka