The table present the different plot functions in the packages mgcv (Wood 2006, 2011) and itsadug for visualizing GAMM models.


Partial effect Sum of all effects Sum of “fixed” effects1
surface plot.gam() vis.gam()
pvisgam() fvisgam()
smooth plot.gam() plot_smooth()
group estimates plot.gam()2 plot_parametric()
random smooths get_random(), inspect_random()

1: include rm.ranef=TRUE to zero all random effects.

2: include all.terms=TRUE to visualize parametric terms.


Examples

library(itsadug)
library(mgcv)
data(simdat)

## Not run: 
# Model with random effect and interactions:
m1 <- bam(Y ~ Group + te(Time, Trial, by=Group)
    + s(Time, Subject, bs='fs', m=1, k=5),
    data=simdat)
# Simple model with smooth:
m2 <-  bam(Y ~ Group + s(Time, by=Group)
    + s(Subject, bs='re')
    + s(Subject, Time, bs='re'),
    data=simdat)

Summary model m1:

gamtabs(m1, type='html')
A. parametric coefficients Estimate Std. Error t-value p-value
(Intercept) 2.0505 0.6881 2.9799 0.0029
GroupAdults 3.1204 0.9731 3.2065 0.0013
B. smooth terms edf Ref.df F-value p-value
te(Time,Trial):GroupChildren 19.7347 21.9398 1651.5571 < 0.0001
te(Time,Trial):GroupAdults 19.6302 21.8781 1466.8432 < 0.0001
s(Time,Subject) 175.2256 178.0000 5803.4508 < 0.0001

Summary model m2:

gamtabs(m2, type='html')
A. parametric coefficients Estimate Std. Error t-value p-value
(Intercept) 2.0575 2.9497 0.6975 0.4855
GroupAdults 3.1264 4.1715 0.7495 0.4536
B. smooth terms edf Ref.df F-value p-value
s(Time):GroupChildren 8.6203 8.9584 1630.0297 < 0.0001
s(Time):GroupAdults 8.8372 8.9921 8798.1199 < 0.0001
s(Subject) 33.9912 34.0000 125335342.8798 < 0.0001
s(Subject,Time) 33.9878 34.0000 125318454.9121 < 0.0001

a. Surfaces

Plotting the partial effects of te(Time,Trial):GroupAdults and te(Time,Trial):GroupChildren.

pvisgam()

par(mfrow=c(1,2), cex=1.1)

pvisgam(m1, view=c("Time", "Trial"), select=1,
        main="Children", zlim=c(-12,12))
pvisgam(m1, view=c("Time", "Trial"), select=2,
        main="Adults", zlim=c(-12,12))

Notes:

  • Plots same data as plot(m1, select=1): partial effects plot, i.e., the plot does not include intercept or any other effects.

  • Make sure to set the zlim values the same when comparing surfaces

  • Use the argument cond to specify the value of other predictors in a more complex interaction. For example, for plotting a modelterm te(A,B,C) use something like pvisgam(model, view=c("A", "B"), select=1, cond=list(C=5)).

fvisgam()

Plotting the fitted effects of te(Time,Trial):GroupAdults and te(Time,Trial):GroupChildren.

par(mfrow=c(1,2), cex=1.1)

fvisgam(m1, view=c("Time", "Trial"), cond=list(Group="Children"),
        main="Children", zlim=c(-12,12), rm.ranef=TRUE)
fvisgam(m1, view=c("Time", "Trial"), cond=list(Group="Adults"),
        main="Adults", zlim=c(-12,12), rm.ranef=TRUE)

Notes:

  • Plots the fitted effects, i.e., the plot includes intercept and estimates for the other modelterms.

  • Make sure to set the zlim values the same when comparing surfaces

  • Use the argument cond to specify the value of other predictors in the model.

  • The argument rm.ranef cancels the contribution of random effects.

  • The argument transform accepts a function for transforming the fitted values into the original scale.

b. Smooths

plot.gam()

Plotting the partial effects of s(Time):GroupAdults and s(Time):GroupChildren.

par(mfrow=c(1,2), cex=1.1)

plot(m2, select=1, shade=TRUE, rug=FALSE, ylim=c(-15,10))
abline(h=0)
plot(m2, select=2, shade=TRUE, rug=FALSE, ylim=c(-15,10))
abline(h=0)

Notes:

  • Not possible to combine different smooth terms in a single plot.

Alternatively:

par(mfrow=c(1,1), cex=1.1)

# Get model term data:
st1 <- get_modelterm(m2, select=1)
st2 <- get_modelterm(m2, select=2)

# plot model terms:
emptyPlot(2000, c(-15,10), h=0,
    main='s(Time)', 
    xmark = TRUE, ymark = TRUE, las=1)
plot_error(st1$Time, st1$fit, st1$se.fit, shade=TRUE)
plot_error(st2$Time, st2$fit, st2$se.fit, shade=TRUE, col='red', lty=4, lwd=2)

# add legend:
legend('bottomleft',
  legend=c('Children', 'Adults'),
  fill=c(alpha('black'), alpha('red')),
  bty='n')

plot_smooth()

Plotting the fitted effects of te(Time,Trial):GroupAdults and te(Time,Trial):GroupChildren i.e., the plot includes intercept and estimates for the other modelterms.

par(mfrow=c(1,2), cex=1.1)

plot_smooth(m1, view="Time", cond=list(Group="Children"),
    rm.ranef=TRUE, ylim=c(-6,10))
plot_smooth(m1, view="Time", cond=list(Group="Adults"),
    col="red", rug=FALSE, add=TRUE,
    rm.ranef=TRUE)

# or alternatively:
plot_smooth(m1, view="Time", plot_all="Group",
    rm.ranef=TRUE)

Notes:

  • Use the argument cond to specify the value of other predictors in the model.

  • The argument rm.ranef cancels the contribution of random effects.

  • The argument transform accepts a function for transforming the fitted values into the original scale.

  • The argument plot_all plots all levels for the given predictor(s).

c. Group estimates

plot.gam()

Plotting the partial effect of grouping predictors such as Group:

par(mfrow=c(1,1), cex=1.1)

plot.gam(m1, select=4, all.terms=TRUE, rug=FALSE)

Alternatively, use get_coefs() to extract the coefficients and plot these:

coefs <- get_coefs(m1)
coefs

par(mfrow=c(1,2), cex=1.1)

b <- barplot(coefs[,1], beside=TRUE, 
             main="Parametric terms",
             ylim=c(0,5))
errorBars(b, coefs[,1], coefs[,2], xpd=TRUE)

# Note that the effect of Group is a *difference* estimate
# between intercept (=GroupChildren) and Group Adults

b2 <- barplot(coefs[1,1], beside=TRUE, 
             main="Estimate for Group",
             ylim=c(0,5), xlim=c(0.1,2.5))
mtext(row.names(coefs), at=b, side=1, line=1)
abline(h=coefs[1,1], lty=2)
rect(b[2]-.4, coefs[1,1], b[2]+.4, coefs[1,1]+coefs[2,1],
     col='gray')
errorBars(b, coefs[,1]+c(0,coefs[1,1]), coefs[,2], xpd=TRUE)

##             Estimate Std. Error
## (Intercept) 2.050539  0.6881195
## GroupAdults 3.120351  0.9731479

Notes:

  • For large models get_coefs() is faster than summary(model).

plot_parametric()

Plotting the fitted effects of grouping predictors such as Group:

pp <- plot_parametric(m1, pred=list(Group=c("Children", "Adults")) )

pp
## $fv
##      Group    Time Trial Subject       fit       CI        rm.ranef
## 1 Children 989.899     0     a01  4.182214 1.705797 s(Time,Subject)
## 2   Adults 989.899     0     a01 11.216176 1.789386 s(Time,Subject)

Random effects

get_random()

For extracting the random effects coefficients (random adjustments of intercept and slopes):

par(mfrow=c(1,2), cex=1.1)
plot(m2, select=3)
plot(m2, select=4)

Or alternatively:

pp <- get_random(m2)

emptyPlot(range(pp[[1]]), range(pp[[2]]), h=0,v=0,
     xlab='Random intercepts', ylab='Random slopes',
     main='Correlation')

text(pp[[1]], pp[[2]], labels=names(pp[[1]]), 
     col='steelblue', xpd=TRUE)

inspect_random()

For plotting and extracting the random smooths:

par(mfrow=c(1,2), cex=1.1)

inspect_random(m1, select=3, main='s(Time, Subject)')

children <- unique(simdat[simdat$Group=="Children", "Subject"])
adults   <- unique(simdat[simdat$Group=="Adults", "Subject"])

inspect_random(m1, select=3, main='Averages', 
      fun=mean, 
      cond=list(Subject=children))
inspect_random(m1, select=3, 
      fun=mean, cond=list(Subject=adults),
      add=TRUE, col='red', lty=5)

# add legend:
legend('bottomleft',
  legend=c('Children', 'Adults'),
  col=c('black', 'red'), lty=c(1,5),
  bty='n')

References

Wood, Simon N. 2006. Generalized Additive Models: An Introduction with R. Chapman; Hall/CRC.
———. 2011. “Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models.” Journal of the Royal Statistical Society (B) 73 (1): 3–36.