Quick overview of plot functions

In itsadug: Interpreting Time Series and Autocorrelated Data Using GAMMs

library(itsadug)
infoMessages("off")

The table present the different plot functions in the packages mgcv [@Wood_2006; @Wood_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()

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

²: 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')

Summary model m2:

gamtabs(m2, type='html')

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)

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

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

itsadug documentation built on June 17, 2022, 5:05 p.m.