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

Question

I'm trying to plot three gam functions (same unit) in the same plot, with x-axis being dates (1 Jan to 31 Dec), y-axis being concentrations.

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

I'm not allowed to insert picture, but basically my y-axis range is now [-5,2]. My question is how can I change the y-axis from smoothed values to actual concentration values? in this case 1~98?

Many thanks in advance!

Answer 1

The general idea is to plot the predicted values from the model, i.e. the smooth effect you have already visualised with plot.gam() plus the model intercept term. You could just add the model intercept to everything (see the shift argument to plot.gam and pass it coef(mod)[1] to get the intercept.), but a more general solution is to predict from the model at a smooth set of time points, and then plot those.

This would also be easier with a single model, not three. I use a single model in the example below, but the ideas apply to separate models, you just need to predict from the three models separately and then combine (e.g. rbind()) the sets of predicted values.

Example data

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)

Now predict over the range of the covariate (x2 in my example, time in your's) for each level of the factor — you'll need to create data accordingly with a column for response with the stacked response values, and a fac (or other name) variable that codes for which type of response it is (macc, pre.pm10, pm10 would be the levels).

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

Then predict from the model at these observations

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

Then plot. (This example has x2 uniform on interval 0,1. To convert this to a day of year as per your example, I'll just multiply 365.25 by x2, but you can just use your time variable directly).

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

Which all produces