Как построить регрессию, преобразованную обратно в исходный масштаб с цветными полосами доверительного интервала?

Question

Я хотел бы построить линию и 95% доверительный интервал из линейной модели, где ответ был git преобразован обратно в исходную шкалу данных. Таким образом, результатом должна быть кривая линия, включающая доверительные интервалы в исходной шкале, где это будет прямая линия в шкале, преобразованной lo git. См. Код:

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

Что я пытаюсь здесь сделать, так это построить кривую линию регрессии и добавить доверительный интервал. В ggplot2 есть функция geom_smooth, в которой может быть указана линейная модель, но я не смог найти способ построения прогнозов из predict.lm(my model).

. Я также хотел бы знать, как добавить цветной многоугольник, который будет представлять доверительный интервал, как на рисунке ниже. Я знаю, что должен использовать функцию polygon и координаты, но я не знаю, как.

Answer 1

Вы можете использовать predict для возрастного диапазона, скажем, 1:100, укажите interval= для CI. Построение с type="l" сгладит красивую кривую. Затем можно добавить доверительные интервалы, используя 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)

Выход

---

Редактировать

Для получения закрашенных полос доверия используйте polygon. Поскольку вам нужны два уровня достоверности, вам, вероятно, нужно сделать один predict ион для каждого. Линия будет покрыта polygon с, поэтому лучше сначала сделать пустой plot, используя type="n" и нарисовать lines в конце. (Обратите внимание, что я также покажу вам некоторые подсказки для пользовательской маркировки осей.) Хитрость для polygons заключается в express значениях вперед и назад, используя 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)))

Выход