Wie füge ich Konfidenzintervalle zum glm-Modell in ggplot hinzu?

Question

Hier ist ein Beispiel dafür, was wie meine Daten aussehen:

DATA <- data.frame(
TotalAbund = sample(1:10), 
TotalHab = sample(0:1), 
TotalInv = sample(c("yes", "no"), 20, replace = TRUE) 
) 
DATA$TotalHab<-as.factor(DATA$TotalHab) 
DATA 

Hier ist mein Modell:

MOD.1<-glm(TotalAbund~TotalInv+TotalHab, family=quasipoisson, data=DATA) 

Hier ist mein Grundstück:

NEWDATA <- with(DATA, 
       expand.grid(TotalInv=unique(TotalInv), 
         TotalHab=unique(TotalHab))) 

pred <- predict(MOD.1,newdata= NEWDATA,se.fit=TRUE) 
gg1 <- ggplot(NEWDATA, aes(x=factor(TotalHab), y=TotalAbund,colour=TotalInv)) 

ich folgende Fehlermeldung erhalten. ..

Error in eval(expr, envir, enclos) : object 'TotalAbund' not found 

... wenn man versucht, die letzte Zeile Code auszuführen:

gg1 + geom_point(data=pframe,size=8,shape=17,alpha=0.7, 
      position=position_dodge(width=0.75)) 

Kann jemand helfen? Wie füge ich 95% Konfidenzintervalle zu meinen Punkten hinzu? Vielen Dank.

Answer 1

Sie müssen die 95% -Konfidenzintervalle selbst berechnen. Sie waren auf dem richtigen Weg mit predict und fragen nach der se.fit. Wir werden zuerst nach den Vorhersagen auf der Verbindungsskala fragen, 95% -Konfidenzintervalle berechnen und sie dann in den realen Maßstab für das Plotten umwandeln. Hier ist eine Komfortfunktion zur Berechnung Ihrer CI für den Log-Link (den Sie im Modell verwendet haben).

# get your prediction 
pred <- predict(MOD.1,newdata= NEWDATA,se.fit=TRUE, 
      type = "link") 

# CI function 
make_ci <- function(pred, data){ 

# fit, lower, and upper CI 
fit <- pred$fit 
lower <- fit - 1.96*pred$se.fit 
upper <- fit + 1.96*pred$se.fit 

return(data.frame(exp(fit), exp(lower), exp(upper), data)) 
} 

my_pred <- make_ci(pred, NEWDATA) 

# to be used in geom_errorbar 
limits <- aes(x = factor(TotalHab), ymax = my_pred$exp.upper., ymin = my_pred$exp.lower., 
        group = TotalInv) 

wir es dann plotten aus, ich werde die endgültige Zwicken Sie verlassen, um die Figur zu machen, wie Sie es wollen.

ggplot(my_pred, aes(x = factor(TotalHab), y = exp.fit., color = TotalInv))+ 
geom_errorbar(limits, position = position_dodge(width = 0.75), 
      color = "black")+ 
geom_point(size = 8, position = position_dodge(width = 0.75), shape = 16)+ 
ylim(c(0,15))+ 
geom_point(data = DATA, aes(x = factor(TotalHab), y = TotalAbund, colour = TotalInv), 
     size = 8, shape = 17, alpha = 0.7, 
     position = position_dodge(width = 0.75))