\HeaderA{smooth.construct}{Constructor functions for smooth terms in a GAM}{smooth.construct}
\aliasA{p.spline}{smooth.construct}{p.spline}
\aliasA{smooth.construct.cc.smooth.spec}{smooth.construct}{smooth.construct.cc.smooth.spec}
\aliasA{smooth.construct.cr.smooth.spec}{smooth.construct}{smooth.construct.cr.smooth.spec}
\aliasA{smooth.construct.cs.smooth.spec}{smooth.construct}{smooth.construct.cs.smooth.spec}
\aliasA{smooth.construct.tensor.smooth.spec}{smooth.construct}{smooth.construct.tensor.smooth.spec}
\aliasA{smooth.construct.tp.smooth.spec}{smooth.construct}{smooth.construct.tp.smooth.spec}
\aliasA{smooth.construct.ts.smooth.spec}{smooth.construct}{smooth.construct.ts.smooth.spec}
\keyword{models}{smooth.construct}
\keyword{smooth}{smooth.construct}
\keyword{regression}{smooth.construct}
\begin{Description}\relax
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{\LinkA{Predict.matrix}{Predict.matrix}} method function (see example code below).

In practice, \code{smooth.construct} is usually called via the wrapper
function \code{\LinkA{smoothCon}{smoothCon}}.
\end{Description}
\begin{Usage}
\begin{verbatim}
smooth.construct(object,data,knots)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{object}] is a smooth specification object, generated by an \code{\LinkA{s}{s}} or \code{\LinkA{te}{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:
\Itemize{
\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).
\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 0 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[full.call] The full version of the \code{s} term, with all defaults expanded explicitly.
\item[label] A suitable label for use with this term.
\item[null.space.dim] The dimension of the null space of the wiggliness
penalty.
\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).
}
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:
\Itemize{
\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[full.call] The full version of the \code{s} term, with all defaults expanded explicitly.
\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[\code{data}] a data frame in which the covariates and any \code{by} variable can be found.
\item[\code{knots}] an optional data frame specifying knot locations for each covariate. If it is null then the knot 
locations are generated automatically.
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
There are built in methods for objects with the following classes:
\code{tp.smooth.spec} (thin plate regression splines); 
\code{cr.smooth.spec} (cubic regression splines);
\code{ts.smooth.spec} (thin plate regression splines with shrinkage-to-zero);
\code{cs.smooth.spec} (cubic regression splines with shrinkage-to-zero);
\code{cc.smooth.spec} (cyclic cubic regression splines);
\code{tensor.smooth.spec} (tensor product smooths).

The returned objects for the built in smooth classes have the following extra elements.
\code{cr.smooth} objects (generated using \code{bs="cr"}) have an additional array \code{xp} giving the knot locations used to generate the basis.

\code{cyclic.smooth} objects (generated using \code{bs="cc"}) have an array \code{xp} of knot locations and a matrix 
\code{BD} used to define the basis (BD transforms function values at the knots to second derivatives at the knots).

\code{tprs.smooth} objects require several items to be stored in order to define the basis. These are:
\Itemize{ 
\item[shift] 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. 
\item[Xu] A matrix of the unique covariate combinations for this smooth (the basis is constructed by first stripping 
out duplicate locations).
\item[UZ] The matrix mapping the t.p.r.s. parameters back to the parameters of a full thin plate spline.
\item[null.space.dimension] The dimension of the space of functions that have zero wiggliness according to the 
wiggliness penalty for this term.
}

Again, these extra elements would be found in the elements of the \code{margin} list of \code{tensor.smooth} 
class object.


\code{tp.smooth.spec} (and \code{ts.smooth.spec}) objects will contain
information on how to handle large datasets in their \code{xt} field. The default is to randomly
subsample 3000 `knots' from which to produce a tprs basis, if the number of
unique predictor variable combinations in excess of 3000. The default can be
modified via the \code{xt} argument to \code{\LinkA{s}{s}}. This is supplied as a
list with elements \code{max.knots} and \code{seed} containing a number
to use in place of 3000, and the random number seed to use (either can be
missing). 


The classes \code{cs.smooth} and \code{ts.smooth} are modifications of the
classes \code{cr.smooth} and \code{tprs.smooth} respsectively: they differ
only in having an extra shrinkage component added to their penalties, so that
automatic model selection can effectively remove such terms from a model altogether.
\end{Details}
\begin{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:
\begin{ldescription}
\item[\code{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.
\item[\code{C}] The matrix defining any constraints on the term - usually a one row matrix giving the column sums of the  
model matrix, which defines the constraint that each term should sum to zero over the covariate values.
\item[\code{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[\code{rank}] an array giving the ranks of the penalties.
\item[\code{df}] the degrees of freedom associated with this term (at least when unpenalized).
\end{ldescription}


Usually the returned object will also include extra information required to define the basis, and used by 
\code{\LinkA{Predict.matrix}{Predict.matrix}} methods to make predictions using the basis. See the \code{Details} section for 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.
\end{Value}
\begin{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.
\end{Section}
\begin{Author}\relax
Simon N. Wood \email{simon.wood@r-project.org}
\end{Author}
\begin{References}\relax
Wood, S.N. (2000)  Modelling and Smoothing Parameter Estimation
with Multiple  Quadratic Penalties. J.R.Statist.Soc.B 62(2):413-428

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


Wood, S.N. (2004a) Stable and efficient multiple smoothing parameter estimation for
generalized additive models. J. Amer. Statist. Ass. 99:673-686

The p-spline code given in the example is based on:

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

\url{http://www.maths.bath.ac.uk/~sw283/}
\end{References}
\begin{SeeAlso}\relax
\code{get.var}, \code{\LinkA{gamm}{gamm}}, \code{\LinkA{gam}{gam}},
\code{\LinkA{Predict.matrix}{Predict.matrix}},
\code{\LinkA{smoothCon}{smoothCon}}, \code{\LinkA{PredictMat}{PredictMat}}
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}# adding "p-spline" classes and methods

smooth.construct.ps.smooth.spec<-function(object,data,knots)
# a p-spline constructor method function
{ require(splines)
  if (length(object$p.order)==1) m<-rep(object$p.order,2) 
  else m<-object$p.order  # m[1] - basis order, m[2] - penalty order
  nk<-object$bs.dim-m[1]  # number of interior knots
  if (nk<=0) stop("basis dimension too small for b-spline order")
  x <- get.var(object$term,data)  # find the data
  xl<-min(x);xu<-max(x);xr<-xu-xl # data limits and range
  xl<-xl-xr*0.001;xu<-xu+xr*0.001;dx<-(xu-xl)/(nk-1) 
  if (!is.null(knots)) k <- get.var(object$term,knots) 
  else k<-NULL
  if (is.null(k)) 
  k<-seq(min(x)-dx*(m[1]+1),max(x)+dx*(m[1]+1),length=nk+2*m[1]+2)   
  if (length(k)!=nk+2*m[1]+2) 
  stop(paste("there should be ",nk+2*m[1]+2," supplied knots"))
  object$X<-spline.des(k,x,m[1]+2,x*0)$design # get model matrix
  if (!object$fixed)       
  { S<-diag(object$bs.dim);if (m[2]) for (i in 1:m[2]) S<-diff(S)
    object$S<-list(t(S)%*%S)  # get penalty
    object$S[[1]] <- (object$S[[1]]+t(object$S[[1]]))/2 # exact symmetry
  }
  object$rank<-object$bs.dim-m[2]  # penalty rank 
  object$null.space.dim <- m[2]  # dimension of unpenalized space  
  object$knots<-k;object$m<-m      # store p-spline specific info.
  object$C<-matrix(colSums(object$X),1,object$bs.dim) #constraint
  object$df<-ncol(object$X)-1      # maximum DoF
  if (object$by!="NA")  # deal with "by" variable 
  { by <- get.var(object$by,data) # find by variable  
    if (is.null(by)) stop("Can't find by variable")
    object$X<-as.numeric(by)*object$X # form diag(by)%*%X
  }
  class(object)<-"pspline.smooth"  # Give object a class
  object
}

Predict.matrix.pspline.smooth<-function(object,data)
# prediction method function for the p.spline smooth class
{ require(splines)
  x <- get.var(object$term,data)
  X <- spline.des(object$knots,x,object$m[1]+2,x*0)$design
  if (object$by != "NA") { ## handle any by variables
        by <- get.var(object$by, data)
        if (is.null(by))
            stop("Can't find by variable")
        X <- as.numeric(by) * X
  }
  X
}

# an example, using the new class....
set.seed(0);n<-400;
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n)*2
y <- f + e
b<-gam(y~s(x0,bs="ps",m=2)+s(x1,bs="ps",m=c(1,3))+
         s(x2,bs="ps",m=2)+s(x3,bs="ps",m=2))
plot(b,pages=1)
# another example using tensor products of the new class

test1<-function(x,z,sx=0.3,sz=0.4)
{ (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
  0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
}
n<-400
x<-runif(n);z<-runif(n);
f <- test1(x,z)
y <- f + rnorm(n)*0.1
b <- gam(y~te(x,z,bs=c("ps","ps"),m=c(2,2)))
vis.gam(b)
\end{ExampleCode}
\end{Examples}

