\name{smooth.construct}
\alias{smooth.construct}
\alias{smooth.construct2}
\alias{user.defined.smooth}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Constructor functions for smooth terms in a GAM}
\description{Smooth terms in a GAM formula are turned into smooth specification objects of 
class \code{xx.smooth.spec} during processing of the formula. Each of these objects is
converted to a smooth object using an appropriate \code{smooth.construct} function. New smooth classes 
can be added by writing a new \code{smooth.construct} method function and a corresponding 
\code{\link{Predict.matrix}} method function (see example code below).

In practice, \code{smooth.construct} is usually called via \code{smooth.construct2} and the wrapper
function \code{\link{smoothCon}}, in order to handle \code{by} variables and
centering constraints (see the \code{\link{smoothCon}} documentation if 
you need to handle these things directly, for a user defined smooth class).

}

\usage{
smooth.construct(object,data,knots)
smooth.construct2(object,data,knots)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{object}{ is a smooth specification object, generated by an \code{\link{s}} or \code{\link{te}} term in a GAM 
formula. Objects generated by \code{s} terms have class \code{xx.smooth.spec} where \code{xx} is given by the 
\code{bs} argument of \code{s} (this convention allows the user to add their own smoothers). 
If \code{object} is not class \code{tensor.smooth.spec} it will have the following elements:

\describe{
\item{term}{The names of the covariates for this smooth, in an array.}
\item{bs.dim}{ Argument \code{k} of the \code{s} term generating the object. This is the dimension of the basis 
used to represent the term (or, arguably, 1 greater than the basis dimension for \code{cc} terms). 
\code{bs.dim<0} indicates that the constructor should set this to the default value.}
\item{fixed}{\code{TRUE} if the term is to be unpenalized, otherwise \code{FALSE}.}
\item{dim}{the number covariates of which this smooth is a function.}
\item{p.order}{the order of the smoothness penalty or \code{NA} for autoselection of this. This is argument 
\code{m} of the \code{s} term that generated \code{object}.}
\item{by}{the name of any \code{by} variable to multiply this term as supplied as an argument to \code{s}. 
\code{"NA"} if there is no such term.}
\item{label}{A suitable label for use with this term.}
\item{xt}{An object containing information that may be needed for basis setup
(used, e.g. by \code{"tp"} smooths to pass optional information on big dataset
handling).}
\item{id}{Any identity associated with this term --- used for linking bases
and smoothing parameters. \code{NULL} by default, indicating no linkage.}
\item{sp}{Smoothing parameters for the term. Any negative are estimated,
otherwise they are fixed at the supplied value. Unless \code{NULL} (default),
over-rides \code{sp} argument to \code{\link{gam}}.}
}

If \code{object} is of class \code{tensor.smooth.spec} then it was generated by a \code{te} term in the GAM formula, 
and specifies a smooth of several variables with a basis generated as a tensor product of lower dimensional bases. 
In this case the object will be different and will have the following elements:
\describe{
\item{margin}{is a list of smooth specification objects of the type listed above, defining the bases which have 
their tensor product formed in order to construct this term.}
\item{term}{is the array of names of the covariates that are arguments of the smooth.}
\item{by}{is the name of any \code{by} variable, or \code{"NA"}.}
\item{fx}{is an array, the elements of which indicate whether (\code{TRUE}) any of the margins in the 
tensor product should be unpenalized.}
\item{label}{A suitable label for use with this term.}
\item{dim}{is the number of covariates of which this smooth is a function.}
\item{mp}{\code{TRUE} if multiple penalties are to be used.}
\item{np}{\code{TRUE} if 1-D marginal smooths are to be re-parameterized in terms of
function values.}
\item{id}{Any identity associated with this term --- used for linking bases
and smoothing parameters. \code{NULL} by default, indicating no linkage.}
\item{sp}{Smoothing parameters for the term. Any negative are estimated,
otherwise they are fixed at the supplied value. Unless \code{NULL} (default),
over-rides \code{sp} argument to \code{\link{gam}}.}
}}

\item{data}{For \code{smooth.construct} a data frame or list containing the evaluation of the elements of \code{object$term},
with names given by \code{object$term}. The last entry will be the \code{by} variable, if \code{object$by}
is not \code{"NA"}. For \code{smooth.construct2} \code{data} need only be an object within which \code{object$term} 
can be evaluated, the variables can be in any order, and there can be irrelevant variables present as well. }

\item{knots}{an optional data frame or list containing the knots relating to \code{object$term}. 
If it is \code{NULL} then the knot locations are generated automatically. The structure of \code{knots} should
be as for \code{data}, depending on whether \code{smooth.construct} or \code{smooth.construct2} is used.}
}

\value{
The input argument \code{object}, assigned a new class to indicate what type of smooth it is and with at least the 
following items added:
\item{X}{The model matrix from this term. This may have an \code{"offset"}
attribute:  a vector of length \code{nrow(X)} containing any contribution of
the smooth to the model offset term. \code{by} variables do not need to be
dealt with here, but if they are then an item \code{by.done} must be added to
the \code{object}.}

\item{S}{A list of positive semi-definite penalty matrices that apply to this term. The list will be empty 
if the term is to be left un-penalized.}

\item{rank}{An array giving the ranks of the penalties.}

\item{null.space.dim}{The dimension of the penalty null space (before centering).}

The following items may be added:

\item{C}{The matrix defining any identifiability constraints on the term, for use when fitting. If this is \code{NULL} then
\code{smoothCon} will add an identifiability constraint that each term should
sum to zero over the covariate values. Set to a zero row matrix if no
constraints are required. If a supplied \code{C} has an attribute \code{"always.apply"} then it is never ignored, even if any
\code{by} variables of a smooth imply that no constraint is actually needed. Code for creating \code{C} should check whether 
the specification object already contains a zero row matrix, and leave this unchanged if it is (since this signifies 
no constraint should be produced). }

\item{Cp}{An optional matrix supplying alternative identifiability constraints for use when predicting. By default the 
fitting constrants are used. This option is useful when some sort of simple sparse constraint is required for fitting, but the 
usual sum-to-zero constraint is required for prediction so that, e.g. the CIs for model components are as narrow as possible. }

\item{no.rescale}{if this is non-NULL then the penalty coefficient matrix of the smooth will not be 
rescaled for enhaced numerical stability (rescaling is the default, because \code{\link{gamm}} requires it). 
Turning off rescaling is useful if the values of the smoothing parameters should be interpretable in a model, 
for example because they are inverse variance components.  }

\item{df}{the degrees of freedom associated with this term (when
unpenalized and unconstrained). If this is null then \code{smoothCon} will set it to the basis 
dimension. \code{smoothCon} will reduce this by the number of constraints.}

\item{te.ok}{\code{0} if this term should not be used as a tensor product marginal, \code{1} if 
it can be used and plotted, and \code{2} is it can be used but not plotted. Set to \code{1} if \code{NULL}.}

\item{plot.me}{Set to \code{FALSE} if this smooth should not be plotted by \code{\link{plot.gam}}.  Set to \code{TRUE} if \code{NULL}.}

\item{side.constrain}{Set to \code{FALSE} to ensure that the smooth is never subject to side constraints as a result of nesting. }

\item{L}{smooths may depend on fewer `underlying' smoothing parameters than there are elements of
\code{S}. In this case \code{L} is the matrix mapping the vector of underlying log smoothing 
parameters to the vector of logs of the smoothing parameters actually multiplying the \code{S[[i]]}. 
\code{L=NULL} signifies that there is one smoothing parameter per \code{S[[i]]}. }

Usually the returned object will also include extra information required to define the basis, and used by 
\code{\link{Predict.matrix}} methods to make predictions using the basis. See
the \code{Details} section for links to the information included for the built in smooth classes. 

\code{tensor.smooth} returned objects will additionally have each element of
the \code{margin} list updated in the same way. \code{tensor.smooths} also
have a list, \code{XP}, containing re-parameterization matrices for any 1-D marginal terms
re-parameterized in terms of function values. This list will have \code{NULL}
entries for marginal smooths that are not re-parameterized, and is only long
enough to reach the last re-parameterized marginal in the list.  

}

\details{ There are built in methods for objects with the following classes:
\code{tp.smooth.spec} (thin plate regression splines: see \code{\link{tprs}}); 
\code{ts.smooth.spec} (thin plate regression splines with shrinkage-to-zero);
\code{cr.smooth.spec} (cubic regression splines: see \code{\link{cubic.regression.spline}};
\code{cs.smooth.spec} (cubic regression splines with shrinkage-to-zero);
\code{cc.smooth.spec} (cyclic cubic regression splines);
\code{ps.smooth.spec} (Eilers and Marx (1986) style P-splines: see \code{\link{p.spline}});
\code{cp.smooth.spec} (cyclic P-splines);
\code{ad.smooth.spec} (adaptive smooths of 1 or 2 variables: see \code{\link{adaptive.smooth}});
\code{re.smooth.spec} (simple random effect terms);
\code{mrf.smooth.spec} (Markov random field smoothers for smoothing over discrete districts);
\code{tensor.smooth.spec} (tensor product smooths).

There is an implicit assumption that the basis only depends on the knots and/or the set of unique 
covariate combinations; i.e. that the basis is the same whether generated from
the full set of covariates, or just the unique combinations of covariates. 

Plotting of smooths is handled by plot methods for smooth objects. A default \code{mgcv.smooth} method 
is used if there is no more specific method available. Plot methods can be added for specific smooth classes, see 
source code for \code{mgcv:::plot.sos.smooth}, \code{mgcv:::plot.random.effect}, \code{mgcv:::plot.mgcv.smooth} 
for example code.

}


\references{ 

Wood, S.N. (2003) Thin plate regression splines. J.R.Statist.Soc.B 65(1):95-114

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for
generalized additive mixed models. Biometrics 62(4):1025-1036

The code given in the example is based on the smooths advocated in:

Ruppert, D., M.P. Wand and R.J. Carroll (2003) Semiparametric Regression. Cambridge 
University Press.

However if you want p-splines, rather than splines with derivative based penalties,
then the built in "ps" class is probably a marginally better bet. It's based on

Eilers, P.H.C. and B.D. Marx (1996) Flexible Smoothing with B-splines and Penalties. 
Statistical Science, 11(2):89-121

}

\author{Simon N. Wood \email{simon.wood@r-project.org}}

\seealso{ \code{\link{s}},\code{\link{get.var}}, \code{\link{gamm}}, \code{\link{gam}},
 \code{\link{Predict.matrix}},
\code{\link{smoothCon}}, \code{\link{PredictMat}}
 }

\section{WARNING}{User defined smooth objects should avoid having attributes names
\code{"qrc"} or \code{"nCons"} as these are used internally to provide
constraint free parameterizations.}

\examples{
## Adding a penalized truncated power basis class and methods
## as favoured by Ruppert, Wand and Carroll (2003) 
## Semiparametric regression CUP. (No advantage to actually
## using this, since mgcv can happily handle non-identity 
## penalties.)

smooth.construct.tr.smooth.spec<-function(object,data,knots) {
## a truncated power spline constructor method function
## object$p.order = null space dimension
  m <- object$p.order[1]
  if (is.na(m)) m <- 2 ## default 
  if (m<1) stop("silly m supplied")
  if (object$bs.dim<0) object$bs.dim <- 10 ## default
  nk<-object$bs.dim-m-1 ## number of knots
  if (nk<=0) stop("k too small for m")
  x <- data[[object$term]]  ## the data
  x.shift <- mean(x) # shift used to enhance stability
  k <- knots[[object$term]] ## will be NULL if none supplied
  if (is.null(k)) # space knots through data
  { n<-length(x)
    k<-quantile(x[2:(n-1)],seq(0,1,length=nk+2))[2:(nk+1)]
  }
  if (length(k)!=nk) # right number of knots?
  stop(paste("there should be ",nk," supplied knots"))
  x <- x - x.shift # basis stabilizing shift
  k <- k - x.shift # knots treated the same!
  X<-matrix(0,length(x),object$bs.dim)
  for (i in 1:(m+1)) X[,i] <- x^(i-1)
  for (i in 1:nk) X[,i+m+1]<-(x-k[i])^m*as.numeric(x>k[i])
  object$X<-X # the finished model matrix
  if (!object$fixed) # create the penalty matrix
  { object$S[[1]]<-diag(c(rep(0,m+1),rep(1,nk)))
  }
  object$rank<-nk  # penalty rank
  object$null.space.dim <- m+1  # dim. of unpenalized space
  ## store "tr" specific stuff ...
  object$knots<-k;object$m<-m;object$x.shift <- x.shift
 
  object$df<-ncol(object$X)     # maximum DoF (if unconstrained)
 
  class(object)<-"tr.smooth"  # Give object a class
  object
}

Predict.matrix.tr.smooth<-function(object,data) {
## prediction method function for the `tr' smooth class
  x <- data[[object$term]]
  x <- x - object$x.shift # stabilizing shift
  m <- object$m;     # spline order (3=cubic)
  k<-object$knots    # knot locations
  nk<-length(k)      # number of knots
  X<-matrix(0,length(x),object$bs.dim)
  for (i in 1:(m+1)) X[,i] <- x^(i-1)
  for (i in 1:nk) X[,i+m+1] <- (x-k[i])^m*as.numeric(x>k[i])
  X # return the prediction matrix
}

# an example, using the new class....
require(mgcv)
set.seed(100)
dat <- gamSim(1,n=400,scale=2)
b<-gam(y~s(x0,bs="tr",m=2)+s(x1,bs="ps",m=c(1,3))+
         s(x2,bs="tr",m=3)+s(x3,bs="tr",m=2),data=dat)
plot(b,pages=1)
b<-gamm(y~s(x0,bs="tr",m=2)+s(x1,bs="ps",m=c(1,3))+
         s(x2,bs="tr",m=3)+s(x3,bs="tr",m=2),data=dat)
plot(b$gam,pages=1)
# another example using tensor products of the new class
dat <- gamSim(2,n=400,scale=.1)$data
b <- gam(y~te(x,z,bs=c("tr","tr"),m=c(2,2)),data=dat)
vis.gam(b)
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ...