Прогноз риска не соответствует ggplot geom_smooth (). Ошибка в моей формуле прогноза? - PullRequest
Новая
       48

Прогноз риска не соответствует ggplot geom_smooth (). Ошибка в моей формуле прогноза?

0 голосов
/ 08 января 2019

Я пытаюсь самостоятельно предсказать значения лесс, предоставляемые ggplot geom_smooth(). Я приложил ссылки на мои данные и график вывода прогнозов.

Данные можно найти здесь . Я следовал примеру из этого поста о прогнозировании лёсса, чтобы воспроизвести значения из ggplot, поэтому я думаю, что я на правильном пути, но что-то упустил. 

library("ggplot2")
load(file="data5a.RData")
lsmod = loess(Flux~DA_SQ_KM, data=data5a, control=loess.control(surface="direct"))
xrange <- max(data5a$DA_SQ_KM,na.rm=TRUE)
xseq <- c(0.01,0.05,0.1,0.2,0.3,0.5,seq(from=1, to=xrange, length=100))
pred = predict(lsmod,newdata=data.frame(DA_SQ_KM = xseq), se=TRUE)
y = pred$fit
ci <- pred$se.fit * qt(0.95 / 2 + .5, pred$df)
ymin = y - ci
ymax = y + ci
loess.DF <- data.frame(x = xseq, y, ymin, ymax, se = pred$se.fit)

ggplot(data5a, aes(DA_SQ_KM, Flux)) + 
  geom_point()+
  geom_smooth(method="loess")+
  geom_smooth(aes_auto(loess.DF), data=loess.DF, stat="identity",col="red")+
  geom_smooth(method="lm",se=FALSE,col="green")+
  theme(legend.position = "bottom")+
  scale_y_log10()+
  scale_x_log10()

Где ошибка в моем коде для воспроизведения данных в синей кривой, которая предсказывается geom_smooth()?

Вот изображение вывода в ggplot:

Output ggplot image

ОБНОВЛЕНИЕ 1:

Я включил сюда обновленный код на основе информации, предоставленной Роландом. Я изменил свой код для использования функции mgcv::gam, поскольку мои данные содержат более 1000 точек. Проблема все еще остается в том, что я не могу воспроизвести модель, созданную geom_smooth в ggplot. Также появилась новая проблема с доверительными интервалами. 

library("ggplot2")
library("mgcv")
load(file="data5a.RData")

#Attempt to re-create the gam model myself
gammod = mgcv::gam(Flux~s(DA_SQ_KM, bs = "cs"),data=data5a)

xrange <- max(data5a$DA_SQ_KM,na.rm=TRUE)
xseq <- c(0.001,0.01,0.05,0.1,0.2,0.3,0.5,seq(from=1, to=xrange, length=100))

pred = predict(gammod ,newdata=data.frame(DA_SQ_KM = xseq), se=TRUE)
y = pred$fit
ci <- pred$se.fit * qt(0.95 / 2 + .5, pred$df)
ymin = y - ci
ymax = y + ci
gam.DF <- data.frame(x = xseq, y, ymin, ymax, se = pred$se.fit)


ggplot(data5a, aes(DA_SQ_KM, Flux)) + 
  geom_point()+
  geom_smooth(aes_auto(gam.DF), data=gam.DF, stat="identity",col="red")+
  stat_smooth(method=mgcv::gam,formula = y ~ s(x, bs = "cs"),se=TRUE,col="purple")+
  theme(legend.position = "bottom")+
  scale_y_log10()+
  scale_x_log10()

Вот выход гаммы в ggplot:

Output2 ggplot image

1 Ответ

0 голосов
/ 09 января 2019

ggplot2 соответствует модели преобразованным переменным, если вы используете scale_* преобразования: 

DF <- data.frame(x = 1:3, y = c(10, 100, 1e3))

library(ggplot2)
p <- ggplot(DF, aes(x, y)) +
  geom_point() +
  scale_y_log10() +
  stat_smooth(method = "lm", n = 3)
g <- ggplot_build(p)
g[["data"]][[2]]
#  x y ymin ymax se PANEL group  colour   fill size linetype weight alpha
#1 1 1    1    1  0     1    -1 #3366FF grey60    1        1      1   0.4
#2 2 2    2    2  0     1    -1 #3366FF grey60    1        1      1   0.4
#3 3 3    3    3  0     1    -1 #3366FF grey60    1        1      1   0.4

Обратите внимание на нулевые SE, которые указывают на идеальную посадку. 

log10(predict(lm(y ~ x, data = DF)))
#  1        2        3 
#NaN 2.568202 2.937016 

predict(lm(log10(y) ~ x, data = DF))
#1 2 3 
#1 2 3
...