Obliczanie przedziałów prognoz dla regresji logistycznej

20

Chciałbym zrozumieć, jak generować przedziały prognoz dla oszacowań regresji logistycznej.

Poradzono mi, aby postępować zgodnie z procedurami zawartymi w Collett's Modeling Binary Data , 2nd Ed str. 98-99. Po wdrożeniu tej procedury i porównaniu jej z R predict.glm, tak naprawdę uważam, że ta książka pokazuje procedurę obliczania przedziałów ufności , a nie przedziałów prognozowania.

Implementacja procedury Collett, w porównaniu do predict.glm, pokazano poniżej.

Chciałbym wiedzieć: jak przejść odtąd do tworzenia przedziału prognozy zamiast przedziału ufności?

#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')])
karbokokacja
źródło
Podstawowe pytanie, dlaczego sqrt (suma (model.vcov * square.student)) przyjęto jako błąd standardowy? Czy nie jest to odchylenie standardowe i należy je podzielić przez sqrt (n)? Jeśli tak, to które n należy zastosować, n zastosować do modelu lub n nowej ramki danych użytej do przewidywania?
Rafael

Odpowiedzi:

6

0<=y<=1

Greg Snow
źródło
6
Szukam 95% przedziału predykcji prognozy, która znajduje się w przestrzeni logarytmicznej. Później przekształcam to w przestrzeń prawdopodobieństwa. Przedział 100% przewidywania nigdy nie byłby interesujący dla żadnej procedury, prawda? Na przykład 100% przedział predykcji dla regresji liniowej obejmowałby -Inf do Inf ... W każdym razie, jak widać w moim kodzie, przedział przewidywania jest obliczany w przestrzeni logarytmicznej, która jest następnie przekształcana w przestrzeń prawdopodobieństwa później . Więc nie sądzę, aby moje pytanie było bezcelowe.
karbokon
2
Iloraz logarytmiczny można przeliczyć na prawdopodobieństwo i można obliczyć przedział ufności dla prawdopodobieństwa (lub logarytmicznego). Ale przedział przewidywania zależy od zmiennej odpowiedzi, która wynosi 0 lub 1. Jeśli twoim wynikiem jest przeżycie z 0 = martwy i 1 = żywy, możesz przewidzieć prawdopodobieństwo przeżycia dla danego zestawu zmiennych towarzyszących i obliczyć przedział ufności dla to prawdopodobieństwo. Ale wynik to 0/1, nie możesz mieć 62% żywego pacjenta, musi to być 0 lub 1, więc jedynymi możliwymi przedziałami prognozy są 0-0, 0-1 i 1-1 (co jest dlaczego większość ludzi trzyma się przedziałów ufności).
Greg Snow,
8
Jeśli masz sytuację, w której odpowiedź jest dwumianowa (która może być agregacją 0-1 w tych samych warunkach), wtedy przedział przewidywania może mieć sens.
Glen_b
7
Regresja logistyczna to regresja prawdopodobieństwa, próbująca modelować prawdopodobieństwo pewnego zdarzenia jako funkcję zmiennych regresora. Przedziały prognostyczne w tym ustawieniu są traktowane jako przedziały na skali prawdopodobieństwa lub skali logarytmiczno-szansowej, co czyni idealne senesy.
kjetil b halvorsen
2
@Cesar, formuła przedziału predykcji jest uzyskiwana przy założeniu, że Y jest normalnie rozmieszczone wokół linii, ale w regresji logistycznej nie mamy rozkładu normalnego, mamy Bernoulliego lub Dwumianowy. Zastosowanie formuł na tej stronie prowadziłoby do przedziału ufności (można to już zrobić) lub do sztucznie poszerzonego przedziału ufności, który nie spełnia definicji przedziału prognozowania (przewidywanie rzeczywistych wyników na oryginalnej skali wyników). Jak wspomniano Glen_b, przedział przewidywania może mieć sens, jeśli wynik jest naprawdę dwumianowy.
Greg Snow,