Calcul des intervalles de prédiction pour la régression logistique

20

J'aimerais comprendre comment générer des intervalles de prédiction pour les estimations de régression logistique.

On m'a conseillé de suivre les procédures décrites dans Collett's Modeling Binary Data , 2nd Ed p.98-99. Après avoir implémenté cette procédure et l'avoir comparée aux R predict.glm, je pense en fait que ce livre montre la procédure de calcul des intervalles de confiance , pas des intervalles de prédiction.

La mise en œuvre de la procédure de Collett, avec une comparaison avec predict.glm, est illustrée ci-dessous.

Je voudrais savoir: comment aller d'ici pour produire un intervalle de prédiction au lieu d'un intervalle de confiance?

#Derived from Collett 'Modelling Binary Data' 2nd Edition p.98-99
#Need reproducible "random" numbers.
seed <- 67

num.students <- 1000
which.student <- 1

#Generate data frame with made-up data from students:
set.seed(seed) #reset seed
v1 <- rbinom(num.students,1,0.7)
v2 <- rnorm(length(v1),0.7,0.3)
v3 <- rpois(length(v1),1)

#Create df representing students
students <- data.frame(
    intercept = rep(1,length(v1)),
    outcome = v1,
    score1 = v2,
    score2 = v3
)
print(head(students))

predict.and.append <- function(input){
    #Create a vanilla logistic model as a function of score1 and score2
    data.model <- glm(outcome ~ score1 + score2, data=input, family=binomial)

    #Calculate predictions and SE.fit with the R package's internal method
    # These are in logits.
    predictions <- as.data.frame(predict(data.model, se.fit=TRUE, type='link'))

    predictions$actual <- input$outcome
    predictions$lower <- plogis(predictions$fit - 1.96 * predictions$se.fit)
    predictions$prediction <- plogis(predictions$fit)
    predictions$upper <- plogis(predictions$fit + 1.96 * predictions$se.fit)


    return (list(data.model, predictions))
}

output <- predict.and.append(students)

data.model <- output[[1]]

#summary(data.model)

#Export vcov matrix 
model.vcov <- vcov(data.model)

# Now our goal is to reproduce 'predictions' and the se.fit manually using the vcov matrix
this.student.predictors <- as.matrix(students[which.student,c(1,3,4)])

#Prediction:
this.student.prediction <- sum(this.student.predictors * coef(data.model))
square.student <- t(this.student.predictors) %*% this.student.predictors
se.student <- sqrt(sum(model.vcov * square.student))

manual.prediction <- data.frame(lower = plogis(this.student.prediction - 1.96*se.student), 
    prediction = plogis(this.student.prediction), 
    upper = plogis(this.student.prediction + 1.96*se.student))

print("Data preview:")
print(head(students))
print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by Collett's procedure:"))
manual.prediction
print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by R's predict.glm:"))    
print(output[[2]][which.student,c('lower','prediction','upper')])
carbocation
la source
Une question de base, pourquoi sqrt (sum (model.vcov * square.student)) est-il supposé comme erreur standard? N'est-ce pas l'écart-type et doit être divisé par sqrt (n)? Dans l'affirmative, quel n doit être utilisé, n utilisé pour ajuster le modèle ou n du nouveau bloc de données utilisé pour prédire?
Rafael

Réponses:

6

0<=y<=1

Greg Snow
la source
6
Je recherche un intervalle de prédiction à 95% d'une prédiction qui est dans un espace de log-odds. Plus tard, je transforme cela en espace de probabilité. Un intervalle de prédiction de 100% ne serait jamais intéressant pour aucune procédure, non? Par exemple, un intervalle de prédiction de 100% pour la régression linéaire comprendrait -Inf à Inf ... Quoi qu'il en soit, comme vous pouvez le voir dans mon code, l'intervalle de prédiction est calculé dans l'espace des cotes logarithmiques, qui est ensuite transformé en espace de probabilité plus tard . Je ne pense donc pas que ma question soit inutile.
carbocation
2
Les log-odds peuvent être convertis en une probabilité et vous pouvez calculer un intervalle de confiance sur la probabilité (ou les log-odds). Mais un intervalle de prédiction est sur la variable de réponse qui est 0 ou 1. Si votre résultat est la survie avec 0 = mort et 1 = vivant, alors vous pouvez prédire la probabilité d'être vivant pour un ensemble donné de covariables et calculer un intervalle de confiance sur cette probabilité. Mais le résultat est 0/1, vous ne pouvez pas avoir un patient vivant à 62%, il doit être 0 ou 1, donc les seuls intervalles de prédiction possibles sont 0-0, 0-1 et 1-1 (ce qui est pourquoi la plupart des gens s'en tiennent aux intervalles de confiance).
Greg Snow
8
Si vous avez une situation où la réponse est binomiale (qui pourrait être un agrégat de 0-1 dans les mêmes conditions), un intervalle de prédiction peut avoir du sens.
Glen_b -Reinstate Monica
7
La régression logistique est la régression d'une probabilité, essayant de modéliser la probabilité d'un événement en fonction des variables du régresseur. Les intervalles de prédiction dans ce paramètre sont considérés comme des intervalles sur l'échelle de probabilité ou l'échelle log-odds, ce qui rend les sénes parfaits.
kjetil b halvorsen
2
@Cesar, la formule d'intervalle de prédiction est dérivée en supposant que Y est normalement distribué autour de la ligne, mais dans la régression logistique, nous n'avons pas de distribution normale, nous avons un Bernoulli ou Binomial. L'application des formules sur cette page entraînerait soit un intervalle de confiance (peut déjà le faire) soit un intervalle de confiance artificiellement élargi qui ne correspond pas à la définition d'un intervalle de prédiction (prédire les résultats réels sur l'échelle de résultats d'origine). Comme Glen_b l'a mentionné, un intervalle de prédiction peut avoir un sens si le résultat est vraiment binomial.
Greg Snow