如何将（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?

非常感谢！

Answer 1

总体思路是绘制模型的预测值，即您已经使用 plot.gam() plus 模型截距项可视化的平滑效果。您可以将模型截距添加到所有内容（请参阅 plot.gam 的 shift 参数并将其传递 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"))

然后剧情。（此示例在间隔 0,1 上具有 x2 统一。要根据您的示例将其转换为一年中的某一天，我只需将 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)

全部产生

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

how to change the value of y axis in (gam) plot from smoothed to actual values?

plot

r

smoothing

gam