R 黄土预测与 ggplot geom_smooth() 不匹配。我的预测公式有误?

R loess prediction does not match ggplot geom_smooth(). Error in my prediction formula?

我正在尝试自行预测 ggplot geom_smooth() 提供的黄土值。我附上了指向我的数据和预测输出图的链接。

可以找到有关黄土预测的数据 here. I followed an example provided from this post 以重现 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 中的输出图像:

更新 1:

我根据 Roland 提供的输入在此处包含了更新的代码。我修改了我的代码以使用 mgcv::gam 函数,因为我的数据包含超过 1000 个点。问题仍然存在,我无法在 ggplot 中重现 geom_smooth 创建的模型。置信区间也出现了一个新问题。

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 中的 gam 输出:

如果您使用 scale_* 转换,

ggplot2 使模型适合转换后的变量:

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