Tracer des intervalles de confiance pour les probabilités prédites à partir d'une régression logistique

20

Ok, j'ai une régression logistique et j'ai utilisé la predict()fonction pour développer une courbe de probabilité basée sur mes estimations.

## LOGIT MODEL:
library(car)
mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))

## PROBABILITY CURVE:
all.x <- expand.grid(won=unique(won), bid=unique(bid))
y.hat.new <- predict(mod1, newdata=all.x, type="response")
plot(bid<-000:1000,predict(mod1,newdata=data.frame(bid<-c(000:1000)),type="response"), lwd=5, col="blue", type="l")

C'est très bien, mais je suis curieux de tracer les intervalles de confiance pour les probabilités. J'ai essayé plot.ci()mais je n'ai pas eu de chance. Quelqu'un peut-il m'indiquer quelques façons d'y parvenir, de préférence avec le carpackage ou la base R.

ATMathew
la source
4
(+1) En réponse aux votes pour clore comme hors sujet: Apparemment, la base de ces votes est que la question semble poser une question purement logicielle ("comment tracer telle ou telle chose dans R"), un question qui devrait en effet apparaître sur SO. Notez cependant que dans la réponse actuelle sont enfouies des formules statistiques pour créer les points de traçage. Cela donne à penser que la question présente un intérêt statistique, donc j'hésite à voter pour la migration. Une bonne réponse ici soulignerait et expliquerait ce point statistique.
whuber

Réponses:

26

Le code que vous avez utilisé estime un modèle de régression logistique à l'aide de la glmfonction. Vous n'avez pas inclus de données, je vais donc en inventer.

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

Un modèle de régression logistique modélise la relation entre une variable de réponse binaire et, dans ce cas, un prédicteur continu. Le résultat est une probabilité transformée par logit en tant que relation linéaire avec le prédicteur. Dans votre cas, le résultat est une réponse binaire correspondant à gagner ou à ne pas gagner au jeu et il est prédit par la valeur du pari. Les coefficients de mod1sont donnés en cotes enregistrées (difficiles à interpréter), selon:

logit(p)=Journal(p(1-p))=β0+β1X1

Pour convertir les cotes enregistrées en probabilités, nous pouvons traduire ce qui précède en

p=exp(β0+β1X1)(1+exp(β0+β1X1))

Vous pouvez utiliser ces informations pour configurer le tracé. Tout d'abord, vous avez besoin d'une plage de la variable de prédiction:

plotdat <- data.frame(bid=(0:1000))

Ensuite, en utilisant predict, vous pouvez obtenir des prédictions basées sur votre modèle

preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)

Notez que les valeurs ajustées peuvent également être obtenues via

mod1$fitted

En spécifiant se.fit=TRUE, vous obtenez également l'erreur standard associée à chaque valeur ajustée. Le résultat data.frameest une matrice avec les composantes suivantes: les prédictions ajustées ( fit), les erreurs-types estimées ( se.fit) et un scalaire donnant la racine carrée de la dispersion utilisée pour calculer les erreurs-types ( residual.scale). Dans le cas d'un logit binomial, la valeur sera 1 (que vous pouvez voir en entrant preddat$residual.scaledans R). Si vous souhaitez voir un exemple de ce que vous avez calculé jusqu'à présent, vous pouvez taper head(data.frame(preddat)).

L'étape suivante consiste à configurer le tracé. J'aime d'abord créer une zone de traçage vierge avec les paramètres:

with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))

Vous pouvez maintenant voir où il est important de savoir comment calculer les probabilités ajustées. Vous pouvez tracer la ligne correspondant aux probabilités ajustées en suivant la deuxième formule ci-dessus. À l'aide de, preddat data.framevous pouvez convertir les valeurs ajustées en probabilités et les utiliser pour tracer une ligne par rapport aux valeurs de votre variable prédictive.

with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))

Enfin, répondez à votre question, les intervalles de confiance peuvent être ajoutés au graphique en calculant la probabilité des valeurs ajustées +/- 1.96multipliée par l'erreur standard:

with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

Le tracé résultant (à partir des données générées aléatoirement) devrait ressembler à ceci:

entrez la description de l'image ici

Par souci d'opportunité, voici tout le code en un seul morceau:

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))
plotdat <- data.frame(bid=(0:1000))
preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)
with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))
with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))
with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

(Remarque: il s'agit d'une réponse fortement modifiée dans le but de la rendre plus pertinente pour stats.stackexchange.)

smillig
la source
où est se.fitdéfinie la variable ?
Macro
Dans predict(..., se.fit=TRUE).
smillig
(-1) Ces IC sont pour chacun des cas individuels? Si c'est le cas, pour un résultat binaire, le seul CI sensible pour une probabilité prédite est [0,1]. Même si cela peut être une réponse techniquement compétente.
rolando2
Selon le commentaire de @ whuber, je pense qu'une bonne réponse devrait inclure une formule pour le calcul de la SE. Quelqu'un pourrait-il peut-être modifier et améliorer la réponse?
Heisenberg
1
Votre réponse semble ne donner que «l'intervalle de prédiction moyen». Comment ajouter «l'intervalle de prédiction ponctuelle»?
Bob Hopez
0

Voici une modification de la solution de @ smillig. J'utilise des outils tidyverse ici, et j'utilise également la linkinvfonction qui fait partie de l'objet modèle GLM mod1. De cette façon, vous n'avez pas à inverser manuellement la fonction logistique, et cette approche fonctionnera quel que soit le GLM spécifique auquel vous correspondez.

library(tidyverse)
library(magrittr)


set.seed(1234)

# create fake data on gambling. Does prob win depend on bid size? 
mydat <- data.frame(
  won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
  bid=runif(250, min=0, max=1000)
)

# logistic regression model: 
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

# new predictor values to use for prediction: 
plotdat <- data.frame(bid=(0:1000))

# df with predictions, lower and upper limits of CIs: 
preddat <- predict(mod1,
               type = "link",
               newdata=plotdat,
               se.fit=TRUE) %>% 
  as.data.frame() %>% 
  mutate(bid = (0:1000), 

         # model object mod1 has a component called linkinv that 
         # is a function that inverts the link function of the GLM:
         lower = mod1$family$linkinv(fit - 1.96*se.fit), 
         point.estimate = mod1$family$linkinv(fit), 
         upper = mod1$family$linkinv(fit + 1.96*se.fit)) 


# plotting with ggplot: 
preddat %>% ggplot(aes(x = bid, 
                   y = point.estimate)) + 
  geom_line(colour = "blue") + 
  geom_ribbon(aes(ymin = lower,
                  ymax = upper), 
              alpha = 0.5) + 
  scale_y_continuous(limits = c(0,1))
Nayef
la source
3
Bien que l'implémentation soit souvent mélangée à un contenu substantiel dans les questions, nous sommes censés être un site pour fournir des informations sur les statistiques, l'apprentissage automatique, etc., pas le code. Il peut également être utile de fournir du code, mais veuillez élaborer votre réponse substantielle dans le texte pour les personnes qui ne lisent pas assez bien cette langue pour reconnaître et extraire la réponse du code.
gung - Rétablir Monica