Plot predicted probabilities and confidence intervals in R

Question

This is a plot I did, I want the confidence intervals for the plot, both upper and lower. I have come so far that I have produced both the upper and lower range but I have problems with the plot that includes the confidence interval.

Here are a few lines of my data, gdk is my binary response and the second variable is the age

               gdk    age prog      calender
29            FALSE    59 NASTK       11
30            FALSE    59 NASTK       10
91             TRUE    49 NMATK        9
129            TRUE    47 NFYSK        8
227           FALSE    46 LARAA       13
244            TRUE    44 LARAA       11
256            TRUE    41 LARAA        9
311           FALSE    38 NMATK        7
323           FALSE    42 NSFYY       11
393            TRUE    40 LARAA       11
449           FALSE    37 NSFYY        9
450           FALSE    38 NSFYY       10

This is the code for my first plot:

prop<-numeric()
for (i in 18:60){prop[i-17]<-mean(both$gdk[both$age==i],na.rm=TRUE)}

mod.red.fin<-glm(respons ~prog+age+calender, family=binomial,data=both)
newdata<-data.frame(prog="NMATK",calender=7, age=18:60)
plot(18:60, predict(mod.red.fin, newdata, type="respons"))

to bring up my confidence, I used the code:

newdata<-data.frame(prog="NMATK",calender=7, age=18:60)
newdata2<-cbind(newdata, predict(mod.red.fin, newdata, type="link", se=TRUE))
newdata2<-within(newdata2, {PredictedProb<-plogis(fit)
                            LL<-plogis(fit-(1.96*se.fit))
                            UL<-plogis(fit+(1.96*se.fit))})

head(newdata2)
   prog    calender age    fit    se.fit  residual.scale        UL        LL   PredictedProb
1 NMATK        7    18 1.637162 0.2128354              1 0.8863833 0.7720644     0.8371484
2 NMATK        7    19 1.569661 0.2072370              1 0.8782376 0.7619639     0.8277353
3 NMATK        7    20 1.502160 0.2032196              1 0.8699448 0.7509808     0.8178965
4 NMATK        7    21 1.434660 0.2008779              1 0.8615687 0.7390311     0.8076263
5 NMATK        7    22 1.367159 0.2002708              1 0.8531708 0.7260410     0.7969207
6 NMATK        7    23 1.299658 0.2014139              1 0.8448057 0.7119527     0.7857774

how do I then plot the confidence interval? Need help with the code. checked a bit on the library(ggplot2) but did not come up with anything.

Answer 1

If you want to use ggplot (probably the easiest way to create your desired plots), use the stat_smooth() geom.

However, you have a problem with your desired plot. You can only have 1 x variable plotted at a time with ggplot. That being said, here's some example code that should get you stated:

d = read.table(header = TRUE, text =          
"              gdk    age prog      calender
             FALSE    59 NASTK       11
             FALSE    59 NASTK       10
              TRUE    49 NMATK        9
              TRUE    47 NFYSK        8
             FALSE    46 LARAA       13
              TRUE    44 LARAA       11
              TRUE    41 LARAA        9
             FALSE    38 NMATK        7
             FALSE    42 NSFYY       11
              TRUE    40 LARAA       11
             FALSE    37 NSFYY        9
             FALSE    38 NSFYY       10")
## Convert gkk from T/F to 1/0
d$gdk2 <- as.numeric(d$gdk)
library(ggplot2)

plot1 <- ggplot(data = d, aes(x = age, y = gdk2)) + 
            stat_smooth(method = 'glm', family = 'binomial') +
            theme_bw()

ggsave('plot1.jpg', plot1, width = 6, height = 4)

Which gives you this plot:

plot2 <- ggplot(data = d, aes(x = calender, y = gdk2)) + stat_smooth(method = 'glm', family = 'binomial') + theme_bw()

ggsave('plot2.jpg', plot2, width = 6, height = 4)   

Which gives you this figure.

BTW, I know ggplot2 can be hard to learn. I would suggest checking out this page for more information. Also, the first Google hit for "confidence ggplot2" was the offical ggplot2 documentation for plotting confidence intervals.