r - 如何将（gam）图中的y轴值从平滑值更改为实际值？

Question

我试图在同一个图中绘制三个 gam 函数（相同的单位），x 轴是日期（1 月 1 日至 12 月 31 日），y 轴是浓度。

## pm, macc and pred in a same plot

gam.pre.pm10.time<-mgcv::gam(pre.pm10~s(time),data=mypred1)
plot(gam.pre.pm10.time,shade=T,xaxt="n",scale=-1,lty=3)
axis(1,labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"),at=seq(1,365,31),las=1)

plot(gam.pm10.time,shade=T,shade.col = "blue", xaxt="n",yaxt="n",xlab="",ylab="PM10", scale = -1)
par(new=TRUE)
plot(gam.macc.time,shade=T,shade.col = "green", xaxt="n",yaxt="n",lty=2,xlab="",ylab="", scale = -1)
par(new=TRUE)
plot(gam.pre.pm10.time,shade=T,shade.col="grey", xaxt="n",yaxt="n",xlab="",ylab="",scale=-1,lty=3)
axis(1,labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"),at=seq(1,365,31),las=1)
legend(x="bottomleft",y=8,bg='transparent',
       legend=c("PM10","MACC","PRED"),
       lty=1:3,cex=0.8)
par(new=FALSE)
#

我不允许插入图片，但基本上我的 y 轴范围现在是 [-5,2]。我的问题是如何将 y 轴从平滑值更改为实际浓度值？在这种情况下1~98？

提前谢谢了！

标签： rplotsmoothinggam

Answer 1

总体思路是绘制模型的预测值，即您已经可视化的平滑效果plot.gam() 加上模型截距项。您可以将模型截距添加到所有内容（请参阅shift参数plot.gam并将其传递coef(mod)[1]以获取截距。），但更通用的解决方案是在一组平滑的时间点从模型进行预测，然后绘制它们。

使用单个模型而不是三个模型也会更容易。我在下面的示例中使用单个模型，但这些想法适用于单独的模型，您只需要分别从三个模型中进行预测，然后组合（例如rbind()）预测值集。

示例数据

library('mgcv')

## Factor `by' variable example (with a spurious covariate x0)
## simulate data...
dat <- gamSim(4)

## fit model...
b <- gam(y ~ fac +s(x2, by = fac), data = dat)

现在预测因子的每个级别的协变量范围（x2在我的示例中，time在您的范围内） - 您需要相应地创建数据，其中包含一个具有堆叠响应值的响应列和一个fac（或其他名称）变量它是哪种类型的响应的代码（macc, pre.pm10,pm10将是级别）。

pdat <- with(dat, expand.grid(fac = levels(fac),
                              x2  = seq(min(x2), max(x2), length = 200)
                             )
            )

然后根据这些观察值从模型中预测

pdat <- transform(pdat, pred = predict(b, newdata = pdat, type = "response"))

然后情节。（这个例子x2在间隔 0,1 上是统一的。要根据你的例子将它转换为一年中的一天，我只需将 365.25 乘以x2，但你可以直接使用你的时间变量）。

## create a `time` variable for plotting
dat  <- transform(dat,  time = 365.25 * x2)
pdat <- transform(pdat, time = 365.25 * x2)

## ylims for plot, contain data
ylims <- with(dat, range(y))

## draw base plot
plot(y ~ time, data = dat, xaxt = 'n')
levs <- levels(dat[['fac']])
cols <- c('red', 'green', 'blue')

## add the fitted lines
for (l in seq_along(levs)) {
    dd <- subset(pdat, fac == levs[l])
    lines(pred ~ time, data = dd, col = cols[[l]])
}

## using *your* code add axis
##  --- this gets the wrong days of year for months bc not all have 31 days!
axis(1, labels = month.abb, at = seq(1, 365, 31), las = 1)

这一切都产生