shaded confident interval for two model in the same plot in Rbase

Question

I have a dataset with 142 observations about probability of individual crossing (nRF01) by different width gap (LargDiscon) at four sites (Site), like that :

  LargDiscon nRF01            Site
1         90     0 Azay-sur-Thouet
2         79     0 Azay-sur-Thouet
3         55     0 Azay-sur-Thouet
4        125     0 Azay-sur-Thouet
5         40     0 Azay-sur-Thouet
6         28     0 Azay-sur-Thouet

I performed a GLM with a binomial error distribution with width gap, site (colored slopes in the plot) and interaction between the two as explanatory variables. I also need the model overall (black slope in the plot) because we used this model for an other analysis.

gt <- glm(nRF01 ~ LargDiscon*Site, data=p, family="binomial") #model with four sites
gt2 <- glm(nRF01 ~ LargDiscon, data=p, family="binomial") #model overall

After that, I plot the model with sites but I also added the overall slope with this following code :

# Plot
par(cex=0.8, cex.lab = 1.5, cex.sub=1) # increase size font
par(cex.lab=1.7, cex.axis=1.9, cex.main=1.9,
    mar = c(5, 5, 5, 5))
plot(jitter(nRF01, 0.1) ~ LargDiscon, 
     data=p, t="n", 
     xlim=c(0, 170), 
     xaxs="i", 
     ylim=c(-0.05, 1.05), 
     xlab="Gap width (m)", 
     ylab="Probability of crossing", 
     main=" ")


lg <- seq(1, 170, 1)
nd <- data.frame(LargDiscon=rep(lg,times=4), Site=rep(levels(p$Site), each=length(lg))) #predictions
nd$pred <- as.vector(predict(gt, nd, type="response"))  # using model to predict values of points 
pred <- as.vector(predict(gt2,  list(LargDiscon=lg), type="response")) #pred of overall model

# add confidence interval
# how ?

abline(h=c(0,1), col="darkgrey")
abline(v=0, col="darkgrey")
col <- c("#FF6101", "#D90000", "#0093A5", "#902E98")
lty <- c("dotted", "dotdash", "twodash", "longdash")

for(i in 1:4) { 
    points(pred ~ LargDiscon, data=nd[nd$Site==levels(p$Site)[i], ], 
           lwd=2, 
           col=col[i],
           lty = lty[i],
           t="l")
    }
points(lg, pred, lwd=2, col="black", t="l")
points(jitter(nRF01, 0.08) ~ LargDiscon, data=p, pch=20, bg="black", col = c("#FF6101", "#D90000", "#0093A5", "#902E98"))
abline(h=0.5, col="darkgrey", lty=2)

But I really need to add confident interval for each model, and I failed. I try with ggplot but I failed too. Could someone help me add the confidence intervals to this graph?

Answer 1

I succeeded with ggplot. I've included the script below. This allows you to have your model + the overall model on the same graph.

 # Create dataframe with fittedvalues
 lg <- seq(1, 170, 1) #width of my explcative variable
    nd <- data.frame(LargDiscon=rep(lg,times=4), Site=rep(levels(p$Site), each=length(lg))) # predictions
    ## add the fitted values by predicting from the model for the new data
    nd <- add_column(nd, fit = predict(gt, newdata = nd, type = 'response'))
    nd <- add_column(nd, wrong_se = predict(gt, newdata = nd, type = 'response',
                                            se.fit = TRUE)$se.fit)
    nd <- mutate(nd, wrong_upr = fit + (2 * wrong_se), wrong_lwr = fit - (2 * wrong_se))
    
    
    # Same for the overall model
    lg <- seq(1, 170, 1) #width of my explcative variable
    nd2 <- data.frame(LargDiscon = lg, fit = as.vector(predict(gt2,  list(LargDiscon=lg), type="response")))
    nd2 <- add_column(nd2, fit = predict(gt2, newdata = nd2, type = 'response'))
    nd2 <- add_column(nd2, wrong_se = predict(gt2, newdata = nd2, type = 'response',
                                            se.fit = TRUE)$se.fit)
    nd2 <- mutate(nd2, wrong_upr = fit + (2 * wrong_se), wrong_lwr = fit - (2 * wrong_se))
    nd2$Site <- "Overall"
    
    
    ## plot it
    
    ggplot() +
      # hline
      geom_hline(yintercept=1, linetype="dashed", color = "darkgrey")+
      geom_hline(yintercept=0, linetype="dashed", color = "darkgrey")+
      geom_hline(yintercept=0.5, linetype="dashed", color = "darkgrey")+
      #raw data
      geom_jitter(data = p, aes(x = LargDiscon, y = nRF01, fill = Site, color = Site), width = 0.005, height = 0.01, alpha = 0.5) +
      # model site
      geom_line(data = nd, aes(x = LargDiscon, y = fit, colour = Site), size = 1) +
      geom_ribbon(data = nd, aes(x = LargDiscon, y = fit, ymin = wrong_lwr, ymax = wrong_upr, fill = Site),
                  alpha = 0.1,
                  colour = NA) +
      # #overall model
      geom_line(data=nd2, aes(x=LargDiscon, y=fit, color = Site), size=1) +
      geom_ribbon(data = nd2, aes(x = LargDiscon, y = fit, ymin = wrong_lwr, ymax = wrong_upr),
                  alpha = 0.1,
                  colour = NA) +
      # geom_ribbon(data = nd2, aes(ymin = wrong_lwr, ymax = wrong_upr),
      #             alpha = 0.1,
      #             colour = NA) +
      # geom_line(nd2, aes(x = LargDiscon, y = fit), size = 1) +
      #scale_colour_discrete(name = 'Site') +
      scale_fill_manual(values = c("#FF6101", "#D90000", "#0093A5", "#902E98", "black"), aesthetics = c("color", "fill")) +
      labs(x = 'Gap width (m)', y = 'Probability of gap crossing') +
      theme_classic()