How2plot nicer GAM curves


Generalized additive models visualize potential non-linear associations between a predictor and a response variable.

The default plotting method produces clean plots for all covariates in the model (or a selection) but: They do not have presentation quality by any means in terms of colors, understandable axes-labels, scaling, etc.

This example shows a customization variant for a additive regression model with one covariate. The goal was to display the absolute value of the response variable on the y-axis and not the “difference from intercept” which is default.

This is only meaningful for a single covariate in the model.

1 The default plot.gam() method

MyGAM1<- with(MyData[MyData$Strata==1,], gam(Y ~ s(Covariate)))
MyGAM0<- with(MyData[MyData$Strata==0,], gam(Y ~ s(Covariate)))

par( mfcol=c(1,2))
plot(MyGAM0)
plot(MyGAM1)

This is the resulting plot:

Stratified Additive Regression Model
Stratified Additive Regression Model

2 The fancy way

1. Extract the values of the model response from the GAM object:

response1 <- predict(MyGAM1, type="response", se.fit=T)
response0 <- predict(MyGAM0, type="response", se.fit=T)

2. Print the response values against the covariate (note: this works just with one covariate)

par(mfcol=c(1,1))

plot(0, type="n", bty="n", main="Fancy GAM plot", xlab="MyCovariate", ylab="MyResponse", lwd=3,ylim=c(0,60), xlim=c(0,200))
legend("bottomright", bty="n", lwd=5, col=c("green","red"), legend=c("Strata = 0", "Strata = 1"))

lines(sm.spline(MyGAM1$model$Covariate , response1$fit) , lwd = 3 , col = "red")
lines(sm.spline(MyGAM1$model$Covariate , response1$fit+1.96*response1$se) , lty = 3 , lwd = 2 , col = "red")
lines(sm.spline(MyGAM1$model$Covariate , response1$fit-1.96*response1$se) , lty = 3 , lwd = 2 , col = "red")

lines(sm.spline(MyGAM0$model$Covariate , response0$fit) , lwd = 3 , col = "green")
lines(sm.spline(MyGAM0$model$Covariate, response0$fit + 1.96 * response0$se) , lty = 3 , lwd = 2, col = "green")
lines(sm.spline(MyGAM0$model$Covariate, response0$fit - 1.96 * response0$se) , lty = 3 , lwd = 2 , col = "green")

abline(h=gam.dm1$coefficients[1], lty=2, lwd=1, col="red")
abline(h=gam.dm0$coefficients[1], lty=2, lwd=1, col="green")

Stratified Additive Regression Model on Response Scale
Stratified Additive Regression Model on Response Scale

 

Advertisement

2 thoughts on “How2plot nicer GAM curves

  1. This produces really neat graphs, but it does not reproduce a plot made with plot.gam(). plotgam() plots confidence bands for the parameter estimates while this code produces confidence intervals for fitted values (which are necessarily wider). This doesn’t make it a bad solution, but one should keep it in mind when interpreting the result.
    See this example code:
    library(mgcv)
    set.seed(0)
    dat <- gamSim(1, n=400, dist="normal", scale=2)
    gam1 <- gam(y~s(x0), data=dat)
    plot(gam1)
    newdat <- data.frame(x0=seq(from=0, to=1, by=0.01))
    fits1 <- predict(gam1, newdata=newdat, type="response", se.fit=TRUE)
    lines(newdat$x0, fits1$fit-mean(predict(gam1, newdata=dat)), col="red")
    lines(newdat$x0, fits1$fit-mean(predict(gam1, newdata=dat)) + 1.96*fits1$se.fit, col="red")
    lines(newdat$x0, fits1$fit-mean(predict(gam1, newdata=dat)) – 1.96*fits1$se.fit, col="red")

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s