gam can use isotropic smooths of any number of variables, specified via terms like
s(x,z,bs="tp",m=3) (or just
s(x,z) as this is the default basis). These terms are based on thin plate
m specifies the order of the derivatives in the thin plate spline penalty.
m is a vector of length 2 and the second element is zero, then the penalty null space of the smooth is not included in the smooth: this is useful if you need to test whether a smooth could be replaced by a linear term, or construct models with odd nesting structures.
Thin plate regression splines are constructed by starting with the basis and penalty for a full thin plate spline and then truncating this basis in an optimal manner, to obtain a low rank smoother. Details are given in Wood (2003). One key advantage of the approach is that it avoids the knot placement problems of conventional regression spline modelling, but it also has the advantage that smooths of lower rank are nested within smooths of higher rank, so that it is legitimate to use conventional hypothesis testing methods to compare models based on pure regression splines. Note that the basis truncation does not change the meaning of the thin plate spline penalty (it penalizes exactly what it would have penalized for a full thin plate spline).
The t.p.r.s. basis and penalties can become expensive to
calculate for large datasets. For this reason the default behaviour is to
max.knots unique data locations if there are more than
such, and to use the sub-sample for basis construction. The sampling is
always done with the same random seed to ensure repeatability (does not
reset R RNG).
max.knots is 2000, by default. Both seed and
max.knots can be modified using the
xt argument to
Alternatively the user can supply knots from which to construct a basis.
"ts" smooths are t.p.r.s. with the penalty modified so that the term is shrunk to zero
for high enough smoothing parameter, rather than being shrunk towards a function in the
penalty null space (see details).
1 2 3 4
a smooth specification object, usually generated by a term
a list containing just the data (including any
a list containing any knots supplied for basis setup — in same order and with same names as
The default basis dimension for this class is
M is the null space dimension
(dimension of unpenalized function space) and
k.def is 8 for dimension 1, 27 for dimension 2 and 100 for higher dimensions.
This is essentially arbitrary, and should be checked, but as with all penalized regression smoothers, results are statistically
insensitive to the exact choise, provided it is not so small that it forces oversmoothing (the smoother's
degrees of freedom are controlled primarily by its smoothing parameter).
The default is to set
m (the order of derivative in the thin plate spline penalty) to the smallest value satisfying
2m > d+1 where
d if the number of covariates of the term: this yields ‘visually smooth’ functions. In any case
2m>d must be satisfied.
The constructor is not normally called directly, but is rather used internally by
To use for basis setup it is recommended to use
For these classes the specification
object will contain
information on how to handle large datasets in their
xt field. The default is to randomly
subsample 2000 ‘knots’ from which to produce a tprs basis, if the number of
unique predictor variable combinations in excess of 2000. The default can be
modified via the
xt argument to
s. This is supplied as a
list with elements
seed containing a number
to use in place of 2000, and the random number seed to use (either can be
For these bases
knots has two uses. Firstly, as mentioned already, for large datasets
the calculation of the
tp basis can be time-consuming. The user can retain most of the advantages of the t.p.r.s.
approach by supplying a reduced set of covariate values from which to obtain the basis -
typically the number of covariate values used will be substantially
smaller than the number of data, and substantially larger than the basis dimension,
k. This approach is
the one taken automatically if the number of unique covariate values (combinations) exceeds
The second possibility
is to avoid the eigen-decomposition used to find the t.p.r.s. basis altogether and simply use
the basis implied by the chosen knots: this will happen if the number of knots supplied matches the
k. For a given basis dimension the second option is
faster, but gives poorer results (and the user must be quite careful in choosing knot locations).
The shrinkage version of the smooth, eigen-decomposes the wiggliness penalty matrix, and sets its zero eigenvalues to small multiples of the smallest strictly positive eigenvalue. The penalty is then set to the matrix with eigenvectors corresponding to those of the original penalty, but eigenvalues set to the peturbed versions. This penalty matrix has full rank and shrinks the curve to zero at high enough smoothing parameters.
An object of class
"ts.smooth". In addition to the usual elements of a
smooth class documented under
smooth.construct, this object will contain:
A record of the shift applied to each covariate in order to center it around zero and avoid any co-linearity problems that might otehrwise occur in the penalty null space basis of the term.
A matrix of the unique covariate combinations for this smooth (the basis is constructed by first stripping out duplicate locations).
The matrix mapping the t.p.r.s. parameters back to the parameters of a full thin plate spline.
The dimension of the space of functions that have zero wiggliness according to the wiggliness penalty for this term.
Simon N. Wood [email protected]
Wood, S.N. (2003) Thin plate regression splines. J.R.Statist.Soc.B 65(1):95-114
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
require(mgcv); n <- 100; set.seed(2) x <- runif(n); y <- x + x^2*.2 + rnorm(n) *.1 ## is smooth significantly different from straight line? summary(gam(y~s(x,m=c(2,0))+x,method="REML")) ## not quite ## is smooth significatly different from zero? summary(gam(y~s(x),method="REML")) ## yes! ## Fool bam(...,discrete=TRUE) into (strange) nested ## model fit... set.seed(2) ## simulate some data... dat <- gamSim(1,n=400,dist="normal",scale=2) dat$x1a <- dat$x1 ## copy x1 so bam allows 2 copies of x1 ## Following removes identifiability problem, by removing ## linear terms from second smooth, and then re-inserting ## the one that was not a duplicate (x2)... b <- bam(y~s(x0,x1)+s(x1a,x2,m=c(2,0))+x2,data=dat,discrete=TRUE) ## example of knot based tprs... k <- 10; m <- 2 y <- y[order(x)];x <- x[order(x)] b <- gam(y~s(x,k=k,m=m),method="REML", knots=list(x=seq(0,1,length=k))) X <- model.matrix(b) par(mfrow=c(1,2)) plot(x,X[,1],ylim=range(X),type="l") for (i in 2:ncol(X)) lines(x,X[,i],col=i) ## compare with eigen based (default) b1 <- gam(y~s(x,k=k,m=m),method="REML") X1 <- model.matrix(b1) plot(x,X1[,1],ylim=range(X1),type="l") for (i in 2:ncol(X1)) lines(x,X1[,i],col=i) ## see ?gam
Loading required package: nlme This is mgcv 1.8-20. For overview type 'help("mgcv-package")'. Family: gaussian Link function: identity Formula: y ~ s(x, m = c(2, 0)) + x Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.02133 0.05315 -0.401 0.689 x 1.18249 0.10564 11.193 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Approximate significance of smooth terms: edf Ref.df F p-value s(x) 0.9334 8 0.304 0.0781 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 R-sq.(adj) = 0.91 Deviance explained = 91.1% -REML = -70.567 Scale est. = 0.012767 n = 100 Family: gaussian Link function: identity Formula: y ~ s(x) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.5600 0.0113 49.56 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Approximate significance of smooth terms: edf Ref.df F p-value s(x) 1.933 2.417 413.1 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 R-sq.(adj) = 0.91 Deviance explained = 91.1% -REML = -69.353 Scale est. = 0.012767 n = 100 Gu & Wahba 4 term additive model
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.