gam
object produced by gam()
and plots the
component smooth functions that make it up, on the scale of the linear
predictor. Optionally produces term plots for parametric model components
as well.## S3 method for class 'gam':
plot(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scale=-1,
n=100,n2=40,pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,
ylab=NULL,main=NULL,ylim=NULL,xlim=NULL,too.far=0.1,
all.terms=FALSE,shade=FALSE,shade.col="gray80",shift=0,
trans=I,seWithMean=FALSE,unconditional=FALSE,by.resids=FALSE,
scheme=0,...)
gam
object as produced by gam()
.TRUE
then partial residuals are added to plots of 1-D smooths. If FALSE
then no residuals are added. If this is an array of the correct length then it is used as the array of
residuals to be used for producing partial residupages=1
then all terms will be plotted on one page with the layout performed automatically.
Set to 0 to have the routine leave all graphics settings as they arselect=2
.ylim
supplied.TRUE
if you want perspective plots for 2-d
terms.TRUE
then the partial effects of parametric
model components are also plotted, via a call to termplot
. Only
terms of order 1 can be plotted in this way.TRUE
to produce shaded regions as confidence bands
for smooths (not avaliable for parametric terms, which are plotted using termplot
).trans
.shift
and trans
are occasionally useful as a means for
getting plots on the response scale, when the model consists only of a single smooth.TRUE
the component smooths are shown with confidence
intervals that include the uncertainty about the overall mean. If FALSE
then the
uncertainty relates purely to the centred smooth itself. Marra and Wood (2012) suggests
tTRUE
then the smoothing parameter uncertainty corrected
covariance matrix is used to compute uncertainty bands, if available. Otherwise the bands
treat the smoothing parameters as fixed.by
variables?
Usually the answer is no, they would be meaningless.plot.gam()
in S-PLUS.Plotting can be slow for models fitted to large datasets. Set rug=FALSE
to improve matters.
If it's still too slow set too.far=0
, but then take care not to overinterpret smooths away from
supporting data.
Plots of 2-D smooths with standard error contours shown can not easily be customized.
The function can not deal with smooths of more than 2 variables!
termplot
.For smooth terms plot.gam
actually calls plot method functions depending on the
class of the smooth. Currently random.effects
, Markov random fields (mrf
),
Spherical.Spline
and
factor.smooth.interaction
terms have special methods (documented in their help files),
the rest use the defaults described below.
For plots of 1-d smooths, the x axis of each plot is labelled
with the covariate name, while the y axis is labelled s(cov,edf)
where cov
is the covariate name, and edf
the estimated (or user defined for regression splines)
degrees of freedom of the smooth. scheme == 0
produces a smooth curve with dashed curves
indicating 2 standard error bounds. scheme == 1
illustrates the error bounds using a shaded
region.
For scheme==0
, contour plots are produced for 2-d smooths with the x-axes labelled with the first covariate
name and the y axis with the second covariate name. The main title of
the plot is something like s(var1,var2,edf)
, indicating the
variables of which the term is a function, and the estimated degrees of
freedom for the term. When se=TRUE
, estimator variability is shown by overlaying
contour plots at plus and minus 1 s.e. relative to the main
estimate. If se
is a positive number then contour plots are at plus or minus se
multiplied
by the s.e. Contour levels are chosen to try and ensure reasonable
separation of the contours of the different plots, but this is not
always easy to achieve. Note that these plots can not be modified to the same extent as the other plot.
For 2-d smooths scheme==1
produces a perspective plot, while scheme==2
produces a heatmap,
with overlaid contours.
Smooths of more than 2 variables are not plotted, but see vis.gam
.
Fine control of plots for parametric terms can be obtained by calling
termplot
directly, taking care to use its terms
argument.
Note that, if seWithMean=TRUE
, the confidence bands include the uncertainty about the overall mean. In other words
although each smooth is shown centred, the confidence bands are obtained as if every other term in the model was
constrained to have average 0, (average taken over the covariate values), except for the smooth concerned. This seems to correspond more closely to how most users interpret componentwise intervals in practice, and also results in intervals with
close to nominal (frequentist) coverage probabilities by an extension of Nychka's (1988) results presented in Marra and Wood (2012).
Sometimes you may want a small change to a default plot, and the arguments to plot.gam
just won't let you do it.
In this case, the quickest option is sometimes to clone the smooth.construct
and Predict.matrix
methods for
the smooth concerned, modifying only the returned smoother class (e.g. to foo.smooth
).
Then copy the plot method function for the original class (e.g. mgcv:::plot.mgcv.smooth
), modify the source code to plot exactly as you want and rename the plot method function (e.g. plot.foo.smooth
). You can then use the cloned
smooth in models (e.g. s(x,bs="foo")
), and it will automatically plot using the modified plotting function.
Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals for Generalized Additive Model Components. Scandinavian Journal of Statistics.
Nychka (1988) Bayesian Confidence Intervals for Smoothing Splines. Journal of the American Statistical Association 83:1134-1143.
Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.
gam
, predict.gam
, vis.gam
library(mgcv)
set.seed(0)
## fake some data...
f1 <- function(x) {exp(2 * x)}
f2 <- function(x) {
0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
}
f3 <- function(x) {x*0}
n<-200
sig2<-4
x0 <- rep(1:4,50)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
e <- rnorm(n, 0, sqrt(sig2))
y <- 2*x0 + f1(x1) + f2(x2) + f3(x3) + e
x0 <- factor(x0)
## fit and plot...
b<-gam(y~x0+s(x1)+s(x2)+s(x3))
plot(b,pages=1,residuals=TRUE,all.terms=TRUE,shade=TRUE,shade.col=2)
plot(b,pages=1,seWithMean=TRUE) ## better coverage intervals
## just parametric term alone...
termplot(b,terms="x0",se=TRUE)
## more use of color...
op <- par(mfrow=c(2,2),bg="blue")
x <- 0:1000/1000
for (i in 1:3) {
plot(b,select=i,rug=FALSE,col="green",
col.axis="white",col.lab="white",all.terms=TRUE)
for (j in 1:2) axis(j,col="white",labels=FALSE)
box(col="white")
eval(parse(text=paste("fx <- f",i,"(x)",sep="")))
fx <- fx-mean(fx)
lines(x,fx,col=2) ## overlay `truth' in red
}
par(op)
## example with 2-d plots, and use of schemes...
b1 <- gam(y~x0+s(x1,x2)+s(x3))
op <- par(mfrow=c(2,2))
plot(b1,all.terms=TRUE)
par(op)
op <- par(mfrow=c(2,2))
plot(b1,all.terms=TRUE,scheme=1)
par(op)
op <- par(mfrow=c(2,2))
plot(b1,all.terms=TRUE,scheme=c(2,1))
par(op)
Run the code above in your browser using DataLab