ggplot2: stat_smooth per risultati logistici con facet_wrap che restituisce modelli 'full' o 'subset' glm

Question

Sto lavorando su un modello di regressione logistica con un predittore continuo e un predittore categoriale con diversi livelli. Voglio presentare i risultati utilizzando ggplot2 e sfruttando lo facet_wrap per mostrare le linee di regressione per ogni livello del predittore categoriale. Facendo questo ho notato che la curva adattata fornita da stat_smooth considera solo i dati in un particolare aspetto, non l'intero set di dati. Questa è una piccola differenza, ma notevole quando si guarda la trama rispetto ai valori previsti restituiti da predict.glm.

Ecco un esempio che ricrea il problema con l'immagine che segue il codice.

library(boot) # needed for inv.logit function 
library(ggplot2) # version 0.8.9 

set.seed(42) 
n <- 100 

df <- data.frame(location = rep(LETTERS[1:4], n), 
       score = sample(45:80, 4*n, replace = TRUE)) 

df$p <- inv.logit(0.075 * df$score + rep(c(-4.5, -5, -6, -2.8), n)) 
df$pass <- sapply(df$p, function(x){rbinom(1, 1, x)}) 

gplot <- ggplot(df, aes(x = score, y = pass)) + 
      geom_point() + 
      facet_wrap(~ location) + 
      stat_smooth(method = 'glm', family = 'binomial') 

# 'full' logistic model 
g <- glm(pass ~ location + score, data = df, family = 'binomial') 
summary(g) 

# new.data for predicting new observations 
new.data <- expand.grid(score = seq(46, 75, length = n), 
         location = LETTERS[1:4]) 

new.data$pred.full <- predict(g, newdata = new.data, type = 'response') 

pred.sub <- NULL 
for(i in LETTERS[1:4]){ 
    pred.sub <- c(pred.sub, 
    predict(update(g, formula = . ~ score, subset = location %in% i), 
      newdata = data.frame(score = seq(46, 75, length = n)), 
      type = 'response')) 
} 

new.data$pred.sub <- pred.sub 

gplot + 
    geom_line(data = new.data, aes(x = score, y = pred.full), color = 'green') + 
    geom_line(data = new.data, aes(x = score, y = pred.sub), color = 'red') 

Quello che ho notato e mi preoccupa è la facilità di vedere in sfaccettatura B. Le curve rosse sono i valori previsti dai modelli solo prendendo in considerazione una posizione, mentre le curve verdi sono previsioni utilizzando il pieno set di dati. I modelli basati sul sottoinsieme dei dati corrispondono alla trama da stat_smooth.

Vorrei tracciare, con sfumature errore standard, le curve verde via ggplot2. Sono sicuro che c'è un'opzione da qualche parte nel codice potrei usare che avrebbe fatto questo, ma devo ancora trovarla, o forse c'è un ordine o un passo dovrei seguire per ottenere le curve verde da una chiamata ggplot diverso. Ho riscontrato problemi simili quando ho tracciato tutto su un lato e utilizzando l'estetica del colore o del gruppo.

Qualsiasi suggerimento sarebbe molto apprezzato.

Answer 1

Hai ragione che il modo per farlo è quello di adattare il modello al di fuori di ggplot2 e quindi calcolare i valori e gli intervalli adattati come ti piace e passare quei dati separatamente.

Un modo per ottenere ciò che si descrive potrebbe essere qualcosa di simile:

preds <- predict(g, newdata = new.data, type = 'response',se = TRUE) 
new.data$pred.full <- preds$fit 

new.data$ymin <- new.data$pred.full - 2*preds$se.fit 
new.data$ymax <- new.data$pred.full + 2*preds$se.fit 

ggplot(df,aes(x = score, y = pass)) + 
    facet_wrap(~location) + 
    geom_point() + 
    geom_ribbon(data = new.data,aes(y = pred.full, ymin = ymin, ymax = ymax),alpha = 0.25) + 
    geom_line(data = new.data,aes(y = pred.full),colour = "blue") 

Questo avviene con i soliti avvertimenti circa intervalli di valori stimati: tocca a voi fare in modo che il l'intervallo che stai tramando è quello che vuoi veramente. Ci tende ad essere un sacco di confusione su "intervalli di previsione".