ggplot2: stat_smooth para resultados logísticos con facet_wrap modelos 'glm completos' o 'subconjuntos' glm

Question

Estoy trabajando en un modelo de regresión logística con un predictor continuo y un predictor categórico con varios niveles. Quiero presentar los resultados usando ggplot2 y explotándome el facet_wrap para mostrar las líneas de regresión para cada nivel del predictor categórico. Al hacer esto, noté que la curva ajustada proporcionada por stat_smooth solo considera los datos en una faceta particular, no el conjunto de datos completo. Esta es una pequeña diferencia, pero notable cuando se mira la trama frente a los valores predichos devueltos desde predict.glm.

Aquí hay un ejemplo que recrea el problema con el gráfico que sigue el código.

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

Lo que observó y me preocupa es la facilidad para ver en la faceta B. Las curvas rojas son los valores de predicción de los modelos considerando sólo un lugar, mientras que las curvas verdes son predicciones usando el pleno conjunto de datos. Los modelos basados ​​en el subconjunto de los datos coinciden con el diagrama de stat_smooth.

Me gustaría trazar, con sombreado de error estándar, las curvas verdes a través de ggplot2. Estoy seguro de que hay una opción en algún lugar del código que podría usar que haría esto, pero aún tengo que encontrarlo, o quizás hay un orden diferente o pasos que debo seguir para obtener las curvas verdes de una llamada ggplot. He encontrado problemas similares al trazar todo en una faceta y usar el color o la estética grupal.

Cualquier sugerencia sería muy apreciada.

Answer 1

Tiene razón en que la forma de hacerlo es ajustar el modelo fuera de ggplot2 y luego calcular los valores ajustados e intervalos como desee y pasar esos datos por separado.

Una forma de lograr lo que usted describe sería algo como esto:

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

Esto viene con las habituales advertencias acerca de los intervalos de valores ajustados: le toca a usted para asegurarse de que la el intervalo que está tramando es lo que realmente quiere. Tiende a haber mucha confusión sobre los "intervalos de predicción".