How to plot regression transformed back on original scale with colored confidence interval bands?

Question

I would like to plot the line and the 95% confidence interval from a linear model where the response has been logit transformed back on the original scale of the data. So the result should be a curved line including the confidence intervals on the original scale, where it would be a straight line on the logit transformed scale. See code:

# Data
dat <- data.frame(c(45,75,14,45,45,55,65,15,3,85),
                  c(.37, .45, .24, .16, .46, .89, .16, .24, .23, .49))
colnames(dat) <- c("age", "bil.")               


# Logit transformation
dat$bb_logit <- log(dat$bil./(1-dat$bil.))

# Model
modelbb <- lm(bb_logit ~ age + I(age^2), data=dat)
summary(modelbb)

# Backtranform
dat$bb_back <- exp(predict.lm(modelbb))/ (1 + exp(predict.lm(modelbb)))

# Plot
plot(dat$age, dat$bb_back)
abline(modelbb)

What do I try here is to plot the curved regression line and add the confidence interval. Within ggplot2 there is the geom_smooth function where the the linear model can be specified, but I could not find a way of plotting the predictions from the predict.lm(my model).

I would also like to know how to add a coloured polygon which will represent the confidence interval as in the image below. I know I have to use function polygon and coordinates but I do not know how.

Answer 1

You may use predict on an age range say 1:100, specify interval= option for the CIs. Plotting with type="l" will smooth a nice curve. Confidence intervals then can be added using lines.

p <- predict(modelbb, data.frame(age=1:100), interval="confidence")
# Backtransform
p.tr <- exp(p) / (1 + exp(p))

plot(1:100, p.tr[,1], type="l", ylim=range(p.tr), xlab="age", ylab="bil.")
sapply(2:3, function(i) lines(1:100, p.tr[,i], lty=2))
legend("topleft", legend=c("fit", "95%-CI"), lty=1:2)

Yields

---

Edit

To get shaded confidence bands use polygon. Since you want two confidence levels you probably need to make one prediction for each. The line will get covered by the polygons, so it's better to make an empty plot first using type="n" and draw the lines at the end. (Note that I'll also show you some hints for custom axis labeling.) The trick for the polygons is to express the values back and forth using rev.

p.95 <- predict(modelbb, data.frame(age=1:100), interval="confidence", level=.95)
p.99 <- predict(modelbb, data.frame(age=1:100), interval="confidence", level=.99)
# Backtransform
p.95.tr <- exp(p.95) / (1 + exp(p.95))
p.99.tr <- exp(p.99) / (1 + exp(p.99))

plot(1:100, p.99.tr[,1], type="n", ylim=range(p.99.tr), xlab="Age", ylab="",
     main="", yaxt="n")
mtext("Tree biomass production", 3, .5)
mtext("a", 2, 2, at=1.17, xpd=TRUE, las=2, cex=3)
axis(2, (1:5)*.2, labels=FALSE)
mtext((1:5)*2, 2, 1, at=(1:5)*.2, las=2)
mtext(bquote(Production ~(kg~m^-2~year^-1)), 2, 2)
# CIs
polygon(c(1:100, 100:1), c(p.99.tr[,2], rev(p.99.tr[,3])), col=rgb(.5, 1, .2),
        border=NA)
polygon(c(1:100, 100:1), c(p.95.tr[,2], rev(p.95.tr[,3])), col=rgb(0, .8, .5),
        border=NA)
# fit
lines(1:100, p.99.tr[,1], ylim=range(p.99.tr), lwd=2)
#legend
legend("topleft", legend=c("fit", "99%-CI", "95%-CI"), lty=c(1, NA, NA), lwd=2,
       pch=c(NA, 15, 15), bty="n",
       col=c("#000000", rgb(.5, 1, .2), rgb(0, .8, .5)))

Yields