Adaptive smooths in GAMs
gam can use adaptive smooths of one or two variables, specified
via terms like s(...,bs="ad",...). (gamm can not use such terms — check out
package AdaptFit if this is a problem.) The basis for such a term is a (tensor product of)
p-spline(s) or cubic regression spline(s). Discrete P-spline type penalties are applied directly to the coefficients
of the basis, but the penalties themselves have a basis representation, allowing the strength of the
penalty to vary with the covariates. The coefficients of the penalty basis are the smoothing parameters.
When invoking an adaptive smoother the k argument specifies the dimension of the smoothing basis (default 40 in 1D, 15 in 2D),
while the m argument specifies the dimension of the penalty basis (default 5 in 1D, 3 in 2D). For an adaptive smooth of two
variables k is taken as the dimension of both marginal bases: different marginal basis dimensions can be
specified by making k a two element vector. Similarly, in the two dimensional case m is the
dimension of both marginal bases for the penalties, unless it is a two element vector, which specifies different
basis dimensions for each marginal (If the penalty basis is based on a thin plate spline then m specifies
its dimension directly).
By default, P-splines are used for the smoothing and penalty bases, but this can be modified by supplying a list as
argument xt with a character vector xt$bs specifying the smoothing basis type.
Only "ps", "cp", "cc" and "cr" may be used for the smoothing basis. The penalty basis is
always a B-spline, or a cyclic B-spline for cyclic bases.
The total number of smoothing parameters to be estimated for the term will be the dimension of the penalty basis.
Bear in mind that adaptive smoothing places quite severe demands on the data. For example, setting m=10 for a
univariate smooth of 200 data is rather like estimating 10 smoothing parameters, each from a data series of length 20.
The problem is particularly serious for smooths of 2 variables, where the number of smoothing parameters required to
get reasonable flexibility in the penalty can grow rather fast, but it often requires a very large smoothing basis
dimension to make good use of this flexibility. In short, adaptive smooths should be used sparingly and with care.
In practice it is often as effective to simply transform the smoothing covariate as it is to use an adaptive smooth.
## S3 method for class 'ad.smooth.spec' smooth.construct(object, data, knots)
object |
a smooth specification object, usually generated by a term |
data |
a list containing just the data (including any |
knots |
a list containing any knots supplied for basis setup — in same order and with same names as |
The constructor is not normally called directly, but is rather used internally by gam.
To use for basis setup it is recommended to use smooth.construct2.
This class can not be used as a marginal basis in a tensor product smooth, nor by gamm.
An object of class "pspline.smooth" in the 1D case or "tensor.smooth" in the 2D case.
Simon N. Wood simon.wood@r-project.org
## Comparison using an example taken from AdaptFit
## library(AdaptFit)
require(mgcv)
set.seed(0)
x <- 1:1000/1000
mu <- exp(-400*(x-.6)^2)+5*exp(-500*(x-.75)^2)/3+2*exp(-500*(x-.9)^2)
y <- mu+0.5*rnorm(1000)
##fit with default knots
## y.fit <- asp(y~f(x))
par(mfrow=c(2,2))
## plot(y.fit,main=round(cor(fitted(y.fit),mu),digits=4))
## lines(x,mu,col=2)
b <- gam(y~s(x,bs="ad",k=40,m=5)) ## adaptive
plot(b,shade=TRUE,main=round(cor(fitted(b),mu),digits=4))
lines(x,mu-mean(mu),col=2)
b <- gam(y~s(x,k=40)) ## non-adaptive
plot(b,shade=TRUE,main=round(cor(fitted(b),mu),digits=4))
lines(x,mu-mean(mu),col=2)
b <- gam(y~s(x,bs="ad",k=40,m=5,xt=list(bs="cr")))
plot(b,shade=TRUE,main=round(cor(fitted(b),mu),digits=4))
lines(x,mu-mean(mu),col=2)
## A 2D example (marked, 'Not run' purely to reduce
## checking load on CRAN).
par(mfrow=c(2,2),mar=c(1,1,1,1))
x <- seq(-.5, 1.5, length= 60)
z <- x
f3 <- function(x,z,k=15) { r<-sqrt(x^2+z^2);f<-exp(-r^2*k);f}
f <- outer(x, z, f3)
op <- par(bg = "white")
## Plot truth....
persp(x,z,f,theta=30,phi=30,col="lightblue",ticktype="detailed")
n <- 2000
x <- runif(n)*2-.5
z <- runif(n)*2-.5
f <- f3(x,z)
y <- f + rnorm(n)*.1
## Try tprs for comparison...
b0 <- gam(y~s(x,z,k=150))
vis.gam(b0,theta=30,phi=30,ticktype="detailed")
## Tensor product with non-adaptive version of adaptive penalty
b1 <- gam(y~s(x,z,bs="ad",k=15,m=1),gamma=1.4)
vis.gam(b1,theta=30,phi=30,ticktype="detailed")
## Now adaptive...
b <- gam(y~s(x,z,bs="ad",k=15,m=3),gamma=1.4)
vis.gam(b,theta=30,phi=30,ticktype="detailed")
cor(fitted(b0),f);cor(fitted(b),f)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.