\name{GAMsetup}
\alias{GAMsetup}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{ Set up GAM using penalized regression splines} 
\description{ Sets up design matrix \eqn{\bf X}{X}, penalty matrices \eqn{ {\bf
S}_i }{S_i} and linear equality constraint matrix \eqn{\bf C}{C} for a GAM defined in terms of 
penalized regression splines. Various other information characterising the bases used is also returned.
 The output is such that the model can be fitted and 
smoothing parameters estimated by the method of Wood (2000) as implemented in routine
\code{mgcv()}. This is usually called by \code{gam}.

}
\usage{
GAMsetup(G)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{G}{ is the single argument to this function: it is a list
  containing at least the elements listed below:
\itemize{
\item{m}{the number of smooth terms in the model}
\item{df}{an array of \code{G$m} integers specifying the maximum d.f. for each spline 
  term.}
\item{n}{the number of data to be modelled}
\item{nsdf}{the number of user supplied columns of the design matrix for any parametric 
model parts}
\item{dim}{An array of dimensions for the smooths. \code{dim[i]} is the
  number of covariates that smooth \code{i} is a function of.}
\item{sp}{array of supplied smoothing parameters. If \code{fit.method} is 
\code{"magic"} then this may be a mixture of positive numbers, which are 
used as the smoothing parameters, and negative to indicate that the 
parameters are to be estimated. With \code{"mgcv"} this is unused.}
\item{fix}{An array of logicals indicating whether each smooth term has fixed degrees of freedom or not.}
\item{s.type}{An array giving the type of basis used for each term. 0 for
  cubic regression spline, 1 for t.p.r.s }
\item{p.order}{An array giving the order of the penalty for each term. 0
for auto selection.}
\item{x}{an array of \code{G$n} element arrays of data and (optionally) design matrix 
columns. The first \code{G$nsdf} elements of \code{G$x} should contain the elements of 
the columns of the design matrix corresponding to the parametric part of the model. The 
remaining \code{G$m} elements of \code{G$x} are the values of the covariates that are 
arguments of the spline terms. Note that the smooths will be centred and no intercept term 
will be added unless an array of 1's is supplied as part of in
\code{G$x}}
\item{vnames}{Array of variable names, including the constant, if present.}
\item{w}{prior weights on response data.}
\item{by}{a 2-d array of \code{by} variables (i.e. covariates that multiply a
smooth term) \code{by[i,j]} is the jth value for the ith \code{by}
variable. There are only as many rows of this array as there are
\code{by} variables in the model (often 0). The rownames of \code{by}
give the \code{by} variable names.}
\item{by.exists}{an array of logicals: \code{by.exists[i]} is
\code{TRUE} if the ith smooth has a \code{by} variable associated with
it, \code{FALSE} otherwise.}
\item{knots}{ a compact array of user supplied knot locations for each smooth, in the order corresponding 
to the  row order in \code{G$x}. There are \code{G$dim[i]} arrays of length \code{G$n.knots[i]} for the ith
         smooth - all these arrays are packed end to end in 1-d array \code{G$knots} - 
        zero length 1 for no knots.}
\item{n.knots}{ array giving number of user supplied knots of basis for each smooth term 0's for none supplied.}
\item{fit.method}{one of \code{"mgcv"} for the Wood (2000) method or 
\code{"magic"} for a more recent and in principle more stable 
alternative.} 
\code{min.sp}{lower bounds on the smoothing parameters: only possible if 
fit method is \code{"magic"}}
\code{H}{the offset penalty matrix, \code{NULL} for none. This is the 
coefficient matrix of any user supplied fixed penalty.}
} % end of itemize
} % end of G
}
\value{
  A list \code{H}, containing the elements of \code{G} (the input list) plus the 
following:  
\item{X}{the full design matrix.}
\item{S}{ If \code{fit.method} is \code{"magic"} then this is a one 
dimensional array containing the non-zero elements of the
   penalty matrices. Let \code{start[k+1]<-start[k]+H$df[1:(k-1)]^2} and
   \code{start[1]<-0}. Then penalty matrix \code{k} has
   \code{H$S[start[k]+i+H$df[i]*(j-1)} on its ith row and jth column.
   To get the kth full penalty matrix the matrix so obtained would be
   inserted into a full matrix of zeroes with it's 1,1 element at \code{H$off[k],H$off[k]}. 
   If \code{fit.method} is \code{"mgcv"} then this is a list of penalty 
matrices, again stored as smallest matrices including all the non-zero 
elements of the penalty matrix concerned. }
\item{off}{is an array of offsets, used to facilitate efficient storage of the penalty 
matrices and to indicate where in the overall parameter vector the parameters of the ith 
spline reside (e.g. first parameter of ith spline is at \code{p[off[i]+1]}).}
\item{C}{ a matrix defining the linear equality constraints on the parameters used to define the
  the model (i.e. \eqn{\bf C}{C} in \eqn{ {\bf Cp } ={\bf 0}}{Cp=0}). }
\item{UZ}{Array containing matrices, which transform from a t.p.r.s. basis to the
  equivalent t.p.s. basis (for t.p.r.s. terms only). The packing method
  is as follows: \cr
set \code{start[1]<-0} and
  \code{start[k+1]<-start[k]+(M[k]+n)*tp.bs[k]} where \code{n} is number
  of data, \code{M[k]} is penalty null space dimension and
  \code{tp.bs[k]} is zero for a cubic regression spline and the basis
  dimension for a t.p.r.s. Then element \code{i,j} of the UZ matrix for
  model term \code{k} is:\cr
 \code{UZ[start[k]+i+(j=1)*(M[k]+n)]}.}
\item{Xu}{Set of unique covariate combinations for each term.  The packing method
  is as follows: \cr
set \code{start[1]<-0} and
  \code{start[k+1]<-start[k]+(xu.length[k])*tp.dim[k]} where \code{xu.length[k]} is number
  of unique covariate combinations and \code{tp.dim[k]} is zero for a
  cubic regression spline
  and the dimension of the smooth (i.e. number of covariates it is a
  function of) for a t.p.r.s. Then element \code{i,j} of the Xu matrix for
  model term \code{k} is:\cr
\code{Xu[start[k]+i+(j=1)*(xu.length[k])]}.}
\item{xu.length}{Number of unique covariate combinations for each
  t.p.r.s. term.}
\item{covariate.shift}{All covariates are centred around zero before
  bases are constructed - this is an array of the applied shifts.}
\item{xp}{matrix whose rows contain the covariate values corresponding to the 
parameters  of each cubic regression spline - the cubic regression splines are parameterized using their \eqn{y}{y}- 
values at a series of \eqn{x}{x} values - these vectors contain those \eqn{x}{x} 
values! Note that these will be covariate shifted.}
\item{rank}{an array giving the ranks of the penalty matrices.}
\item{m.free}{this is only for use with \code{"magic"} and is the number 
of smoothing parameters that must be estimated.}
\item{m.off}{again only for \code{"magic"}: the offests for the penalty 
matrices for the penalties with smoothing parameters that must be 
estimated.} }
\references{ 

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

\url{http://www.stats.gla.ac.uk/~simon/}
}

\author{Simon N. Wood \email{simon@stats.gla.ac.uk}}

\seealso{ \code{\link{mgcv}} \code{\link{gam}} }

\examples{
set.seed(0)
n<-100 # number of observations to simulate
x <- runif(5 * n, 0, 1) # simulate covariates
x <- array(x, dim = c(5, n)) # put into array for passing to GAMsetup
pi <- asin(1) * 2  # begin simulating some data
y <- 2 * sin(pi * x[2, ])
y <- y + exp(2 * x[3, ]) - 3.75887
y <- y + 0.2 * x[4, ]^11 * (10 * (1 - x[4, ]))^6 + 10 * (10 * 
     x[4, ])^3 * (1 - x[4, ])^10 - 1.396
sig2<- -1    # set magnitude of variance 
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- y + e          # simulated data
w <- matrix(1, n, 1) # weight matrix
par(mfrow = c(2, 2)) # scatter plots of simulated data
plot(x[2, ], y);plot(x[3, ], y);plot(x[4, ], y);plot(x[5, ], y)
x[1,]<-1
# create list for passing to GAMsetup....
G <- list(m = 4, n = n, nsdf = 0, df = c(15, 15, 15, 15),dim=c(1,1,1,1),
     s.type=c(0,0,0,0),by=0,by.exists=c(FALSE,FALSE,FALSE,FALSE),
     p.order=c(0,0,0,0),x = x,n.knots=rep(0,4),fit.method="mgcv")
H <- GAMsetup(G)
H$y <- y    # add data to H
H$sig2 <- sig2  # add variance (signalling GCV use in this case) to H
H$w <- w       # add weights to H
H$sp<-array(-1,H$m)
H$fix<-array(FALSE,H$m)
H$conv.tol<-1e-6;H$max.half<-15
H$min.edf<-5;H$fixed.sp<-0
H <- mgcv(H)  # select smoothing parameters and fit model    
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ...


\eof
\name{QT}
\alias{QT}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{ QT factorisation of a matrix }
\description{ If \eqn{\bf A}{A} is an \eqn{r \times c}{r by c} matrix, 
where \eqn{r<c}{r<c}, then it can be factorised: 
\eqn{ {\bf AQ}=[{\bf 0},{\bf T}] }{ AQ=[0,T] }, 
where \eqn{\bf Q}{Q} is an orthogonal matrix and \eqn{\bf T}{T}
is reverse lower triangular. The first \eqn{c-r}{c-r} columns of \eqn{\bf Q}{Q} 
from an orthogonal basis for the null space of \eqn{\bf A}{A}. \eqn{\bf Q}{Q} is 
returned as a sequence of \eqn{r}{r} Householder rotations (from the right), as 
described below.
}
\usage{
QT(A)
}
%- maybe also `usage' for other objects documented here.
\arguments{
 \item{A}{ is an \R matrix having more columns than rows }
}
\details{ This function is primarily useful for providing the null space of the 
linear constraint matrix \eqn{\bf C}{C}, from the linear constraints 
\eqn{ {\bf Cp} = {\bf 0} }{ Cp=0 }, as a series of Householder rotations of the form 
used internally by \code{mgcv()}. It does not need to be called to set up a problem 
for solution by mgcv. 

   
}
\value{

The returned matrix \eqn{\bf R}{R}, say, is of the same dimension as the input 
matrix \eqn{\bf A}{A}. Each of its rows contains a vector, \eqn{ {\bf u}_i }{ u_i }, 
defining one Householder rotation, \eqn{ {\bf H}_i =({\bf I} - {\bf u}_i {\bf 
u}_i^\prime) }{ H_i = (I - u_i u_i') }. The orthogonal matrix \eqn{\bf Q}{Q} is 
defined by: \eqn{ {\bf Q}={\bf H}_1 {\bf H}_2 \ldots {\bf H}_r }{ Q=H_1 H_2 ... 
H_r }.
 
}
\references{Gill, Murray and Wright (1981) Practical Optimization, Academic %%@
Press}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}


\examples{
set.seed(0)
y<-rnorm(12,1,1)
A<-matrix(y,2,6)
B<-QT(A)
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..










\eof
\name{exclude.too.far}
\alias{exclude.too.far}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Exclude prediction grid points too far from data}
\description{ Takes two arrays defining the nodes of a grid over a 2D covariate space and two arrays 
defining the location of data in that space, and returns a logical vector with elements \code{TRUE} if 
corresponding node is too far from data and \code{FALSE} otherwise. Basically a service routine for 
\code{vis.gam}.
}
\usage{
exclude.too.far(g1,g2,d1,d2,dist)
}
%- maybe also `usage' for other objects documented here.
\arguments{ 
\item{g1}{co-ordinates of grid relative to first axis.}
\item{g2}{co-ordinates of grid relative to second axis.}
\item{d1}{co-ordinates of data relative to first axis.}
\item{d2}{co-ordinates of data relative to second axis.}
\item{dist}{how far away counts as two far. Grid and data are first scaled so that the grid lies exactly 
in the unit square, and \code{dist} is a distance within this unit square.} 
}
\details{ Linear scalings of the axes are first determined so that the grid defined by the nodes in 
\code{g1} and \code{g2} lies exactly in the unit square (i.e. on [0,1] by [0,1]). These scalings are 
applied to \code{g1}, \code{g2}, \code{d1} and \code{d2}. The minimum Euclidean 
distance from each node to a datum is then determined and if it is greater than \code{dist} the 
corresponding entry in the returned array is set to \code{TRUE} (otherwise to \code{FALSE}). The 
calculations are performed by looping through the nodes (although they are vectorized with respect to the 
data): unless there are very few data this is not overly ineffeicient, and avoids using excessive amounts of 
memory. 
}

\value{An logical array with \code{TRUE} indicating a node in the grid defined by \code{g1}, \code{g2} that 
is too far from any datum to be included. 
}

\references{
\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\section{WARNING }{ Routine is vectorized w.r.t. the data and not the grid: it will be very inefficient if 
\code{d1} and \code{d2} are very short while \code{g1} and \code{g2} are very long.
} 

\seealso{ \code{\link{vis.gam}} }

\examples{
library(mgcv)
x<-rnorm(100);y<-rnorm(100) # some "data"
n<-20 # generate a grid....
mx<-seq(min(x),max(x),length=n)
my<-seq(min(y),max(y),length=n)
gx<-rep(mx,n);gy<-rep(my,rep(n,n))
tf<-exclude.too.far(gx,gy,x,y,0.1)
plot(gx[!tf],gy[!tf],pch=".");points(x,y,col=2)
}
\keyword{hplot}%-- one or more ...







\eof
\name{full.score}
\alias{full.score}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{GCV/UBRE score for use withn nlm}
\description{ \code{gam} allows the option to `polish' smoothing parameter 
estimates by minimising the GCV/UBRE score calculated at convergence of the IRLS 
algorithm, given a set of smoothing parameters, rather than using the `performance 
iteration' method which estimates smoothing parameters within the IRLS loop. The 
estimates are often slightly different, since the power iteration effectively 
neglects the depedence of the iterative weights on the smoothing parameters.

The `polishing' optimization is fairly crude and numerically costly. 
\code{\link{nlm}} is used to minimize the scores with respect to the smoothin 
parameters, and this routine is designed to be supplied to it as an argument. 

This is basically a service routine for \code{\link{gam}}, and is not usually 
called directly by users.
}
\usage{
full.score(sp,G,family,control,gamma)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{sp}{The logs of the smoothing parameters}

\item{G}{a list reurned by \code{\link{GAMsetup}}}

\item{family}{The family object for the GAM.}

\item{control}{a list returned be \code{\link{gam.control}}}

\item{gamma}{the degrees of freedom inflation factor (usually 1).}
}
\value{ The value of the GCV/UBRE score, with attribute \code{"full.gam.object"}
which is the full object returned by \code{\link{gam.fit}}.
}
 
}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}


\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..










\eof
\name{gam}
\alias{gam}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Additive Models using penalized regression splines and 
GCV}
\description{ Fits the specified  generalized additive model (GAM) to
  data. \code{gam()} is not a clone of what Splus provides.
Smooth terms are represented using penalized regression splines
with smoothing parameters selected by GCV/UBRE or by regression splines with
fixed degrees of freedom (mixtures of the two are
permitted). Multi-dimensional smooths are available using penalized thin plate
regression splines, but the user must make sure that covariates are sensibly scaled
relative to each other when using such terms. For a general overview see
Wood (2001). For more on specifying models see \code{\link{gam.models}}. For more on model 
selection see \code{\link{gam.selection}}.
}
\usage{

gam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
    na.action,control=gam.control(),scale=0,knots=NULL,sp=NULL,
    min.sp=NULL,H=NULL,gamma=1,...)
}
%- maybe also `usage' for other objects documented here.
\details{ 
A generalized additive model (GAM) is a generalized linear model (GLM) in which the linear 
predictor is given by a user specified sum of smooth functions of the covariates plus a 
conventional parametric component of the linear predictor. A simple example is:
\deqn{\log(E(y_i)) = f_1(x_{1i})+f_2(x_{2i})}{log(E(y_i))=f_1(x_1i)+f_2(x_2i)}
where the (idependent) response variables \eqn{y_i \sim {\rm Poi }}{y_i~Poi}, and
\eqn{f_1}{f_1} and \eqn{f_2}{f_2} are smooth functions of covariates \eqn{x_1}{x_1} and 
\eqn{x_2}{x_2}. The log is an example of a link function. 

If absolutely any smooth functions were allowed in model fitting then maximum likelihood 
estimation of such models would invariably result in complex overfitting estimates of 
\eqn{f_1}{f_1}  and \eqn{f_2}{f_2}. For this reason the models are usually fit by 
penalized likelihood 
maximization, in which the model (negative log) likelihood is modified by the addition of 
a penalty for each smooth function, penalizing its `wiggliness'. To control the tradeoff 
between penalizing wiggliness and penalizing badness of fit each penalty is multiplied by 
an associated smoothing parameter: how to estimate these parameters, and 
how to practically represent the smooth functions are the main statistical questions 
introduced by moving from GLMs to GAMs. 

The \code{mgcv} implementation of \code{gam} represents the smooth functions using 
penalized regression splines, and by default uses basis functions for these splines that 
are designed to be optimal, given the number basis functions used. The smooth terms can be 
functions of any number of covariates and the user has some control over how smoothness of 
the functions is measured. 

\code{gam} in \code{mgcv} solves the smoothing parameter estimation problem by using the 
Generalized Cross Validation (GCV) criterion or an Un-Biased Risk Estimator criterion 
(UBRE) which is in practice an approximation to AIC. Smoothing parameters are chosen to 
minimize the GCV or UBRE score for the model, and the main computational challenge solved 
by the \code{mgcv} package is to do this efficiently and reliably. Two alternative 
numerical methods are provided, see \code{\link{mgcv}}, \code{\link{magic}} and 
\code{\link{gam.control}}. 

Broadly \code{gam} works by first constructing basis functions and a quadratic penalty 
coefficient matrix for each smooth term in the model formula, obtaining a model matrix for 
the strictly parametric part of the model formula, and combining these to obtain a 
complete model matrix (/design matrix) and a set of penalty matrices for the smooth terms. 
Some linear identifiability constraints are also obtained at this point. The model is 
fit using \code{\link{gam.fit}}, a modification of \code{\link{glm.fit}}. The GAM 
penalized likelihood maximization problem is solved by penalized Iteratively 
Reweighted  Least Squares (IRLS) (see e.g. Wood 2000). At each iteration a penalized 
weighted least squares problem is solved, and the smoothing parameters of that problem are 
estimated by GCV or UBRE. Eventually both model parameter estimates and smoothing 
parameter estimates converge. 

The fitting approach just described, in which the smoothing parameters are estimated for 
each approximating linear model of the IRLS process was suggested by Chong Gu (see, e.g. 
Gu 2002), and is very computationally efficient. However, because the approach neglects 
the dependence of the iterative weights on the smoothing parameters, it is usually 
possible to find smoothing parameters which actually yield slightly lower GCV/UBRE score 
estimates than those resulting from this `performance iteration'. \code{gam} therefore also 
allows the user to `improve' the smoothing parameter estimates, by using O'Sullivan's 
(1986) suggested method, in which for each trial set of smoothing parameters the IRLS is 
iterated to convergance before the UBRE/GCV score is evaluated. This requires much less 
efficient minimisation of the power iteration based on \code{\link{nlm}}, and is 
therefore quite slow. 





Two alternative bases are available for representing model
  terms. Univariate smooth terms can be represented using conventional
  cubic regression splines - which are very efficient
  computationally - or thin plate regression splines. Multivariate terms
  must be represented using thin plate regression splines. For either
  basis the user specifies the dimension of the basis for each smooth
  term. The dimension of the basis is one more than the maximum degrees of freedom that the 
  term can have, but usually the term will be fitted by penalized
  maximum likelihood estimation and the actual degrees of freedom will be
  chosen by GCV. However, the user can choose to fix the degrees of
  freedom of a term, in which case the actual degrees of freedom will be
  one less than the basis dimension.

  Thin plate regression splines are constructed by starting with the
  basis 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. The t.p.r.s. basis can become expensive to
  calculate for large datasets. In this case the user can supply a reduced 
  set of knots to use in basis construction (see knots, in the argument list).
  
  In the case of the cubic regression spline basis, knots  of the spline are placed evenly
  throughout the covariate values to which the term refers:  For
  example, if fitting 101 data with an 11 knot spline of \code{x} then
  there would be a knot at every 10th (ordered)  \code{x} value. The
  parameterization used represents the spline in terms of its
   values at the knots. The values at neighbouring knots
     are connected by sections of  cubic polynomial constrainted to be 
     continuous up to and including second derivative at the knots. The resulting curve
     is a natural cubic  spline through the values at the knots (given two extra conditions specifying 
     that the second derivative of the curve should be zero at the two end 
     knots). This parameterization gives the parameters a nice interpretability. 

     
     Details of \code{"mgcv"} GCV/UBRE minimization method are given in Wood (2000): 
the basis of the approach 
     is to alternate efficient global optimization with respect to one overall smoothing 
     parameter with Newton updates of a set of relative smoothing parameters for each smooth term. The Newton updates are
backed up by steepest descent, since the GCV/UBRE score functions are not positive definite everywhere.

 }

\arguments{ 
\item{formula}{ A GAM formula (see also \code{\link{gam.models}}). This is exactly like the formula for a
glm except that smooth terms can be added to the right hand side of the
formula (and a formula of the form \code{y ~ .} is not allowed).
Smooth terms are specified by expressions of the form: \cr
\code{s(var1,var2,...,k=12,fx=FALSE,bs="tp",by=a.var)} where \code{var1},
\code{var2}, etc. are the covariates which the smooth
is a function of and \code{k} is the dimension of the basis used to
represent the smooth term. If \code{k} is not
 specified then \code{k=10*3^(d-1)} is used where \code{d} is the number
 of covariates for this term. \code{fx} is used to indicate whether or
 not this term has a fixed muber of degrees of freedom (\code{fx=FALSE}
 to select d.f. by GCV/UBRE). \code{bs} indicates the basis to use, with
 \code{"cr"} indicating cubic regression spline, and \code{"tp"}
 indicating thin plate regression spline: \code{"cr"} can only be used
 with 1-d smooths. \code{by} can be used to specify a variable by which
 the smooth should be multiplied. For example \code{gam(y~z+s(x,by=z))}
   would specify a model \eqn{ E(y) = f(x)z}{E(y)=f(x)z} where
   \eqn{f(\cdot)}{f(.)} is a smooth function (the formula is
   \code{y~x+s(x,by=z)} rather than \code{y~s(x,by=z)} because in GAMs
   the smooths are always set up to sum to zero over the covariate
   values). The \code{by} option is particularly useful for models in
   which different functions of the same variable are required for
   each level of a factor: see \code{\link{s}}.

For backwards compatibility the formula may also include terms like 
\code{s(x,12|f)}, which specifies a regression spline which is not to be penalized
and has 12 knots, or \code{s(x,z,25)} indicating a rank 25 penalized
t.p.r.s. In such cases arguements \code{k}, \code{fx} and \code{bs} are
ignored if supplied and a one dimensional term will always use a cubic
regression spline basis. Note that a term of the form \code{s(x)} will
result in a term with a \code{"tp"} basis.

Formulae can involve nested or "overlapping" terms such as \cr
\code{y~s(x)+s(z)+s(x,z)} or \code{y~s(x,z)+s(z,v)}: see
\code{\link{gam.side.conditions}} for further details and examples. } 
 \item{family}{
This is a family object specifying the distribution and link to use in
fitting etc. See \code{\link{glm}} and \code{\link{family}} for more
details. The negative binomial families provided by the MASS library 
can be used, with or without known \eqn{\theta}{theta} parameter: see
\code{\link{gam.neg.bin}} for details.
} 

\item{data}{ A data frame containing the model response variable and covariates required by the
formula. If this is missing then the frame from which \code{gam} was called is
searched for the variables specified in the formula.} 

\item{weights}{ prior weights on the data.}

\item{subset}{ an optional vector specifying a subset of observations to be
          used in the fitting process.}

\item{na.action}{ a function which indicates what should happen when the data
          contain `NA's.  The default is set by the `na.action' setting
          of `options', and is `na.fail' if that is unset.  The
          ``factory-fresh'' default is `na.omit'.}

\item{control}{A list of fit control parameters returned by 
\code{\link{gam.control}}.}

\item{scale}{ If this is zero then GCV is used for all distributions
except Poisson and binomial where UBRE is used with scale parameter
assumed to be 1. If this is greater than 1 it is assumed to be the scale
parameter/variance and UBRE is used: to use the negative binomial in this case
\eqn{\theta}{theta} must be known. If \code{scale} is negative  GCV 
is always used, which means that the scale parameter will be estimated by GCV and the Pearson 
estimator, or in the case of the negative binomial \eqn{\theta}{theta} will be estimated 
in order to force the GCV/Pearson scale estimate to unity (if this is possible). For binomial models in 
particular, it is probably worth  comparing UBRE and GCV results; for ``over-dispersed Poisson'' GCV is
probably more appropriate than UBRE.} 

\item{knots}{this is an optional list containing user specified knot values to be used for basis construction. 
For the \code{cr} basis the user simply supplies the knots to be used, and there must be the same number as the basis
dimension, \code{k}, for the smooth concerned. For the \code{tp} basis \code{knots} has two uses. Firstly, for large datasets 
the calculation of the \code{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, \code{k}. 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 
basis dimension, \code{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). 
Different terms can use different 
numbers of knots, unless they share a covariate.
}

\item{sp}{A vector of smoothing parameters for each term can be provided here.
 Smoothing parameters must 
be supplied in the order that the smooth terms appear in the model 
formula. With fit method \code{"magic"} (see \code{\link{gam.control}} 
and \code{\link{magic}}) then negative elements indicate that the 
parameter should be estimated, and hence a mixture of fixed and estimated 
parameters is possible. With fit method \code{"mgcv"}, if \code{sp} is 
supplied then all its elements must be positive.  }

\item{min.sp}{for fit method \code{"magic"} only, lower bounds can be 
supplied for the smoothing parameters. Note that if this option is used then
the smoothing parameters, \code{sp}, in the returned object will need to be added to
what is supplied here to get the actual smoothing parameters. Lower bounds on the smoothing 
parameters can sometimes help stabilize otherwise divergent P-IRLS iterations.}

\item{H}{a user supplied fixed quadratic penalty on the parameters of the 
GAM can be supplied, with this as its coefficient matrix. A common use of this term is 
to add a ridge penalty to the parameters of the GAM in circumstances in which the model
is close to un-identifiable on the scale of the linear predictor, but perfectly well
defined on the response scale.}

\item{gamma}{It is sometimes useful to inflate the model degrees of 
freedom in the GCV or UBRE score by a constant multiplier. This allows 
such a multiplier to be supplied if fit method is \code{"magic"}.} 
\item{...}{further arguments for 
passing on e.g. to \code{gam.fit}} }

\value{ 
 The function returns an object of class \code{"gam"} which has the following elements: 

\item{boundary}{did parameters end up at boundary of parameter space?} 

\item{by}{a 2-d array of \code{by} variables (i.e. covariates that multiply a
smooth term) \code{by[i,j]} is the jth value for the ith \code{by}
variable. There are only as many rows of this array as there are
\code{by} variables in the model (often 0). The rownames of \code{by}
give the \code{by} variable names.}

\item{by.exists}{an array of logicals: \code{by.exists[i]} is
  \code{TRUE} if the ith smooth has a \code{by} variable associated with
it, \code{FALSE} otherwise.}

\item{call}{the matched call (allows \code{update} to be used with \code{gam} objects, for example). }

\item{coefficients}{the coefficients of the fitted model. Parametric
          coefficients are  first, followed  by coefficients for each
          spline term in turn.}

\item{converged}{indicates whether or not the iterative fitting method converged.} 

\item{covariate.shift}{covariates get shifted so that they are centred
  around zero - this is by how much.}

\item{deviance}{(unpenalized)}

\item{df}{The maximum degrees of freedom for each smooth term plus one. In the \code{"cr"} basis case 
this is the number of knots used. For the \code{"tp"} basis it is the rank of the smooth (one greater than the 
maximum degrees of freedom because of the centering constraint).}

\item{df.null}{null degrees of freedom.} 

\item{dim}{number of covariates of which term is a function.}

\item{edf}{estimated degrees of freedom for each smooth.}

\item{family}{family object specifying distribution and link used.}

\item{fit.method}{ The underlying multiple GCV/UBRE method used: \code{"magic"} 
for the new more stable method, \code{"mgcv"} for the Wood (2000) method.}

\item{fitted.values}{fitted model predictions of expected value for each
          datum.}

\item{formula}{the model formula.}

\item{full.formula}{the model formula with each smooth term fully
  expanded and with option arguments given explicitly (i.e. not with reference
  to other variables) - useful for later prediction from the model.}

\item{gcv.ubre}{The minimized GCV or UBRE score.}

\item{gcv.used}{TRUE if GCV used for smoothing parameter selection, FALSE if UBRE used.}

\item{hat}{array of elements from the leading diagonal of the `hat' (or `influence') matrix. Same length as response data vector.}

\item{iter}{number of iterations of P-IRLS taken to get convergence.}

\item{linear.predictor}{fitted model prediction of link function of
expected value for  each datum.}

\item{mgcv.conv}{ A list of covergence diagnostics relating to smoothing
parameter estimation. Differs for method \code{"magic"} and \code{"mgcv"}. Here is 
the \code{"mgcv"} version:

\item{score}{corresponding to edf, an array of GCV or UBRE scores for the model given the final 
     estimated relative smoothing parameters.}

\item{g}{the gradient of the GCV/UBRE score w.r.t. the relative smoothing parameters at termination.}

\item{h}{the second derivatives corresponding to \code{g} above - i.e. the leading diagonal of the Hessian.}

\item{e}{the eigen-values of the Hessian. All non-negative indicates a positive definite Hessian.}

\item{iter}{the number of iterations taken.}

\item{in.ok}{\code{TRUE} if the second smoothing parameter guess improved the GCV/UBRE score.}

\item{step.fail}{\code{TRUE} if the algorithm terminated by failing to improve the GCV/UBRE score rather than by "converging". 
Not necessarily a problem, but check the above derivative information quite carefully.}

In the case of \code{"magic"} the items are:

\item{full.rank}{The apparent rank of the problem given the model matrix and 
constraints.}

\item{rank}{The numerical rank of the problem.}

\item{fully.converged}{\code{TRUE} is multiple GCV/UBRE converged by meeting 
convergence criteria. \code{FALSE} if method stopped with a steepest descent step 
failure.}

\item{hess.pos.def}{Was the hessian of the gCV/UBRE score positive definite at 
smoothing parameter estimation convergence?}

\item{iter}{How many iterations were required to find the smoothing parameters?}

\item{score.calls}{and how many times did the GCV/UBRE score have to be evaluated?}

\item{rms.grad}{root mean square of the gradient of the GCV/UBRE score at 
convergence.}

} % end of mgcv.conv listing

\item{min.edf}{Minimum possible degrees of freedom for whole model.}

\item{model}{model frame containing all variables needed in original model fit.}

\item{nsdf}{number of parametric, non-smooth, model terms including the
          intercept.}

\item{null.deviance}{deviance for single parameter model.} 

\item{p.order}{the order of the penalty used for each term. 0 signals
  auto-selection.}

\item{prior.weights}{prior weights on observations.} 

\item{residuals}{the deviance residuals for the fitted model.}

\item{sig2}{estimated or supplied variance/scale parameter.}

\item{sp}{smoothing parameter for each smooth.}

\item{s.type}{type of spline basis used: 0 for conventional cubic regression spline, 1 for t.p.r.s.}

\item{UZ}{array storing the matrices for transforming from t.p.r.s. basis to
  equivalent t.p.s. basis - see \link{GAMsetup} for details of how the
  matrices are packed in this array.}

\item{Vp}{estimated covariance matrix for parameter estimators.}

\item{weights}{final weights used in IRLS iteration.}

\item{xp}{knot locations for each cubic regression spline based smooth.
  \code{xp[i,]+covariate.shift[i]} are the locations for the ith smooth.}

\item{Xu}{The set of unique covariate locations used to define t.p.s. from which
  t.p.r.s. basis was derived. Again see \link{GAMsetup} for details of
  the packing algorithm.}

\item{xu.length}{The number of unique covariate combinations in the
  data. }

\item{y}{response data.}

}

\references{

Key References:

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

Background References:

Green and Silverman (1994) Nonparametric Regression and Generalized  Linear Models. Chapman and Hall.
  
Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

Gu (2002) Smoothing Spline ANOVA Models, Springer.

Hastie and Tibshirani (1990) Generalized Additive Models. Chapman and Hall.

O'Sullivan, Yandall and Raynor (1986) Automatic smoothing of regression functions in generalized linear models.
J. Am. Statist.Ass. 81:96-103 

Wahba (1990) Spline Models of Observational Data. SIAM 

Wood (2001) mgcv:GAMs and Generalized Ridge Regression for R. R News 1(2):20-25
   
Wood and Augustin (2002) GAMs with integrated model selection using penalized regression splines and applications 
to environmental modelling. Ecological Modelling 157:157-177

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\section{WARNINGS }{

The \code{"mgcv"} code does not check for rank defficiency of the
model matrix that may result from lack of identifiability between the
parametric and smooth components of the model. 

You must have more unique combinations of covariates than the model has total
parameters. (Total parameters is sum of basis dimensions plus sum of non-spline 
terms less the number of spline terms). 

Automatic smoothing parameter selection is not likely to work well when 
fitting models to very few response data.

Relative scaling of covariates to a multi-dimensional smooth term has an
affect on the results: make sure that relative scalings are
sensible. For example, measuring one spatial co-ordinate in millimetres
and the other in lightyears will usually produce poor results.   

With large datasets (more than a few thousand data) the \code{"tp"}
basis gets very slow to use: use the \code{knots} argument as discussed above and 
shown in the examples. Alternatively, for 1-d smooths  you can use the \code{"cr"} basis.

For data with many  zeroes clustered together in the covariate space it is quite easy to set up 
GAMs which suffer from identifiability problems, particularly when using poisson or binomial
families. The problem is that with e.g. log or logit links, mean value zero corresponds to
an infinite range on the linear predictor scale. Some regularization is possible in such cases: see 
\code{\link{gam.control}} for details.
} 

\seealso{\code{\link{gam.models}}, \code{\link{s}}, \code{\link{predict.gam}},
  \code{\link{plot.gam}}, \code{\link{summary.gam}}, \code{\link{gam.side.conditions}},
\code{\link{gam.selection}},\code{\link{mgcv}}, \code{\link{gam.control}}
\code{\link{gam.check}}, \code{\link{gam.neg.bin}}, \code{\link{magic}},\code{\link{vis.gam}}
 }

\examples{
library(mgcv)
set.seed(0) 
n<-400
sig2<-4
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
pi <- asin(1) * 2
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, 0, sqrt(abs(sig2)))
y <- f + e
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))
summary(b)
plot(b,pages=1)
# an extra ridge penalty (useful with convergence problems) ....
bp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),H=diag(0.5,41)) 
print(b);print(bp);rm(bp)
# set the smoothing parameter for the first term, estimate rest ...
bp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),sp=c(0.01,-1,-1,-1))
plot(bp,pages=1);rm(bp)
# set lower bounds on smoothing parameters ....
bp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),min.sp=c(0.001,0.01,0,10)) 
print(b);print(bp);rm(bp)

# now a GAM with 3df regression spline term & 2 penalized terms
b0<-gam(y~s(x0,k=4,fx=TRUE,bs="tp")+s(x1,k=12)+s(x2,15))
plot(b0,pages=1)
# now fit a 2-d term to x0,x1
b1<-gam(y~s(x0,x1)+s(x2)+s(x3))
par(mfrow=c(2,2))
plot(b1)
par(mfrow=c(1,1))
# now simulate poisson data
g<-exp(f/5)
y<-rpois(rep(1,n),g)
b2<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson)
plot(b2,pages=1)
# and a pretty 2-d smoothing example....
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<-500
old.par<-par(mfrow=c(2,2))
x<-runif(n);z<-runif(n);
xs<-seq(0,1,length=30);zs<-seq(0,1,length=30)
pr<-data.frame(x=rep(xs,30),z=rep(zs,rep(30,30)))
truth<-matrix(test1(pr$x,pr$z),30,30)
contour(xs,zs,truth)
y<-test1(x,z)+rnorm(n)*0.1
b4<-gam(y~s(x,z))
fit1<-matrix(predict.gam(b4,pr,se=FALSE),30,30)
contour(xs,zs,fit1)
persp(xs,zs,truth)
vis.gam(b4)
par(old.par)
# very large dataset example using knots
n<-10000
x<-runif(n);z<-runif(n);
y<-test1(x,z)+rnorm(n)
ind<-sample(1:n,1000,replace=FALSE)
b5<-gam(y~s(x,z,k=50),knots=list(x=x[ind],z=z[ind]))
vis.gam(b5)
# and a pure "knot based" spline of the same data
b6<-gam(y~s(x,z,k=100),knots=list(x= rep((1:10-0.5)/10,10),
        z=rep((1:10-0.5)/10,rep(10,10))))
vis.gam(b6,color="heat")
}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..








\eof
\name{gam.check}
\alias{gam.check}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Some diagnostics for a fitted gam model.}
\description{ Takes a fitted \code{gam} object produced by \code{gam()} and produces some diagnostic information
about the fitting procedure and results. 
}
\usage{
gam.check(b)
}
%- maybe also `usage' for other objects documented here.
\arguments{ 
\item{b}{ a fitted \code{gam} object as produced by \code{gam()}.}
}
\details{ This function plots 4 standard diagnostic plots, and some other convergence diagnostics. Output differs depending on whether the 
underlying fitting method was \code{"mgcv"} or \code{"magic"} (see \code{\link{gam.control}}).

For fit method \code{"mgcv"}, the first plot shows the GCV or UBRE score against model 
degrees of freedom given the final estimates of the relative smoothing parameters for the model. This is a slice through the 
GCV/UBRE score function that passes through the minimum found during fitting. Although not conclusive (except in the single 
smoothing parameter case), a lack of multiple local minima on this plot is suggestive of a lack of multiple local minima in the 
GCV/UBRE function and is therefore a good thing. Multiple local minima on this plot indicates that the GCV/UBRE function 
may have multiple local minima, but in a multiple smoothing parameter case this is not conclusive - multiple local minima on one slice 
through a function do not necessarily imply that the function has multiple local minima. A "good" plot here is a smooth curve with 
only one local minimum (which is therefore its global minimum). In any case, note that possible problems with multiple local minima are 
substantially reduced by the global line search with respect to the overall model smoothing parameter that is performed at 
each iteration of the minimizer (and produces the data used in this plot). 

The location of the minimum used for the fitted model is also marked on the first plot. Sometimes this location may be a local minimum
that is not the global minimum on the plot. There is a legitimate reason for this to happen, and it does not always indicate problems.
Smoothing parameter selection is based on applying GCV/UBRE to the approximating linear model produced by the GLM IRLS fitting method
employed in \code{gam.fit()}. It is sometimes possible for these approximating models to develop "phantom" minima in their GCV/UBRE scores. These 
minima usually imply a big change in model parameters, and have the characteristic that the minimia will not be present in the GCV/UBRE score 
of the approximating model that would result from actually applying this parameter change. In other words, these are spurious minima in regions 
of parameter space well beyond those for which the weighted least squares problem can be expected to represent the real underlying likelihood well.
Such minima can lead to convergence problems. To help ensure convergence even in the presence of phantom minima, 
\code{gam.fit} switches to a cautious optimization mode after a user controlled number of iterations of the IRLS algorithm (see \link{gam.control}).
In the presence of local minima in the GCV/UBRE score, this method selects the minimum that leads to the smallest change in 
model estimated degrees of freedom. This approach is usually sufficient to deal with phantom minima. Setting \code{trace} to \code{TRUE} in 
\code{gam.control} will allow you to check exactly what is happening. 

If the location of the point indicating the minimum is not on the curve showing the GCV/UBRE function then there are numerical problems with the 
estimation of the effective degrees of freedom: this usually reflects problems with  the relative scaling of covariates that are arguments of a single smooth. 
In this circumstance reported estimated degrees of freedom can not be trusted, although the fitted model and term estimates are likely to be quite 
acceptable. 

If the fit method is \code{"magic"} then there is no global search and the problems with phantom local minima are much reduced. The first plot in this case will simply be a normal QQ plot of the standardixed residuals.

The other 3 plots are two residual plots and plot of fitted values against original data.

The function also prints out information about the convergence of the GCV minimization algorithm, indicating how 
many iterations  were required to minimise the GCV/UBRE score. A message is printed if the minimization terminated by 
failing to improve the score with a steepest descent step: otherwise minimization terminated by meeting convergence criteria.
The mean absolute gradient or RMS gradient of the GCV/UBRE function at the minimum is given. An indication of whether or not the Hessian of the GCV/UBRE function is positive definite is given. If some smoothing parameters 
are not well defined (effectively zero, or infinite) then it may not be, although this is not usually a problem. If the fit method is 
\code{"mgcv"}, a message is printed 
if the second guess smoothing parameters did not improve on the first guess - this is primarily there for the developers. 
For fit method \code{"magic"} the estimated rank of the model is printed.


The ideal results from this function have a smooth, single minima GCV/UBRE plot, good residual plots, and convergence to small 
gradients with a positive definite Hessian. However, failure to meet some of these criteria 
is often  acceptable,  and the information provided is primarily of use in diagnosing suspected problems. 
High gradients at convergence  are a clear indication of problems, however.

Fuller data can be extracted from \code{mgcv.conv} part of the \code{gam} object.  

}

\references{

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

\url{http://www.stats.gla.ac.uk/~simon/}

}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}


\seealso{  \code{\link{gam}} \code{\link{mgcv}}}

\examples{
library(mgcv)
set.seed(0)
n<-200
sig2<-4
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
pi <- asin(1) * 2
y <- 2 * sin(pi * x0)
y <- y + exp(2 * x1) - 3.75887
y <- y+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- y + e
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))
plot(b,pages=1)
gam.check(b)
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ...







\eof
\name{gam.control}
\alias{gam.control}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Setting GAM fitting defaults}
\description{ This is an internal function of package \code{mgcv} which allows control of the numerical options for fitting a GAM. 
Typically users will want to modify the defaults if model fitting fails to converge, or if the warnings are generated which suggest a 
loss of numerical stability during fitting.  
}
\usage{
gam.control(irls.reg=0.0,epsilon = 1e-04, maxit = 20,globit = 20,
            mgcv.tol=1e-6,mgcv.half=15,nb.theta.mult=10000, trace = FALSE,
            fit.method="magic",perf.iter=TRUE,rank.tol=.Machine$double.eps^0.5) 
}
\arguments{ 
\item{irls.reg}{For most models this should be 0. The iteratively re-weighted least squares method
by which GAMs are fitted  can fail to converge in some circumstances. For example data with many zeroes can cause 
problems in a model with a log link, because a mean of zero corresponds to an infinite range of linear predictor 
values. Such convergence problems are caused by a fundamental lack of identifiability, but do not show up as 
lack of identifiability in the penalized linear model problems that have to be solved at each stage of iteration.
In such circumstances it is possible to apply a ridge regression penalty to the model to impose identifiability, and 
\code{irls.reg} is the size of the penalty. The penalty can only be used if \code{fit.method=="magic"}.
}

\item{epsilon}{This is used for judging conversion of the GLM IRLS loop in \code{gam.fit}.}

\item{maxit}{Maximum number of IRLS iterations to perform using cautious GCV/UBRE opttimization, after \code{globit} 
IRLS iterations with normal GCV optimization have been performed. Note that fit method \code{"magic"} makes no distinction 
between cautious and global optimization.}

\item{globit}{Maximum number of IRLS iterations to perform with normal GCV/UBRE optimization. If convergence is not achieved after these 
iterations then a further \code{maxit} iterations will be performed using cautious GCV/UBRE optimization. } 

\item{mgcv.tol}{The convergence tolerance parameter to use in GCV/UBRE optimization.}

\item{mgcv.half}{If a step of  the GCV/UBRE optimization method leads to a worse GCV/UBRE score, then the step length is halved. This is
the number of halvings to try before giving up.}

\item{nb.theta.mult}{Controls the limits on theta when negative binomial parameter is to be estimated. Maximum theta is set to 
the initial value multiplied by \code{nb.theta.mult}, while the minimum value is set to the initial value divided by 
\code{nb.theta.mult}. }

\item{trace}{Set this to \code{TRUE} to turn on diagnostic output.}

\item{fit.method}{set to \code{"mgcv"} to use the method described in Wood (2000). Set to \code{"magic"} to use a 
newer numerically more stable method, which allows regularization and mixtures of fixed and estimated smoothing parameters. 
Set to \code{"fastest"} to use \code{"mgcv"} for one smooth models and \code{"magic"} otherwise.}

\item{perf.iter}{set to \code{TRUE} to use Gu's approach to finding smoothing 
parameters in which GCV or UBRE is
applied to the penalized linear modelling problem produced at each IRLS iteration. This method is very fast, but means that it 
is often possible to find smoothing parameters yielding a slightly lower GCV/UBRE score by perturbing some of the estimated 
smoothing parameters a little. Technically this occurs because the performance iteration effectively neglects the dependence of the 
iterative weights on the smoothing parameters. Setting \code{perf.iter} to \code{FALSE} uses O'Sullivan's approach in which
the IRLS is run to convergence for each trial set of smoothing parameters. This is much slower, since it uses \code{nlm} with 
finite differences and requires many more IRLS steps. } 

\item{rank.tol}{The tolerance used to estimate rank when using \code{fit.method="magic"}.}
}

\details{ 
With fit method \code{"mgcv"}, 
\code{maxit} and \code{globit} control the maximum iterations of the IRLS algorithm, as follows: 
the algorithm will first execute up to
\code{globit} steps in which the GCV/UBRE algorithm performs a global search for the best overall 
smoothing parameter at every iteration. If convergence is not achieved within \code{globit} iterations, then a further 
\code{maxit} steps are taken, in which the overall smoothing parameter estimate is taken as the 
one locally minimising the GCV/UBRE score and resulting in the lowest EDF change. The difference 
between the two phases is only significant if the GCV/UBRE function develops more than one minima. 
The reason for this approach is that the GCV/UBRE score for the IRLS problem can develop "phantom"
minimima for some models: these are minima which are not present in the GCV/UBRE score of the IRLS
problem resulting from moving the parameters to the minimum! Such minima can lead to convergence 
failures, which are usually fixed by the second phase. 
}


\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}


}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}


\seealso{   \code{\link{gam}}, \code{\link{gam.fit}}, \code{\link{glm.control}} }

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..





\eof
\name{gam.fit}
\alias{gam.fit}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Additive Models fitting using penalized regression splines and 
GCV}
\description{ This is an internal function of package \code{mgcv}. It is a modification
  of the function \code{glm.fit}, designed to be called from \code{gam}. The major
  modification is that rather than solving a weighted least squares problem at each IRLS step, 
  a weighted, penalized least squares problem
  is solved at each IRLS step with smoothing parameters associated with each penalty chosen by GCV or UBRE,
  using routine \code{mgcv}. For further information on usage see code for \code{gam}. Some regularization of the 
IRLS weights is also permitted as a way of addressing identifiability related problems (see 
\code{\link{gam.control}}). Negative binomial parameter estimation is supported.
}
}



\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\section{WARNING }{The code does not check for rank defficiency of the
model matrix -
      it will likely just fail instead!
 } 

\seealso{   \code{\link{gam}} \code{\link{mgcv}}}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..

\eof
\name{gam.models}
\alias{gam.models}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Specifying generalized Additive Models.}
\description{ This page is intended to provide some more information on
  how to specify GAMs. Assume that in general we are interested in
  modelling some response variable \eqn{y_i}{y} with an exponential
  family distribution (e.g. Normal, Poisson, binomial, etc.) using predictor variables
  \eqn{x_{1i}, x_{2i}, \ldots}{x1,x2,...}, and let \eqn{\mu_i \equiv
    E(y_i)}{m=E(y)}. A typical GAM might be:
  \deqn{g(\mu_i) = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} +
    f_1(x_{3i})+f_2(x_{4i},x_{4i})}{g(m)=b0 + b1x1 + b2x2 + f1(x3) +
    f2(x4,x5) }
  where \eqn{g} is a smooth monotonic "link" function, and \eqn{f_1}{f1}
  and \eqn{f_2}{f2} are smooth functions. The key idea here is that the
  dependence of the response on the predictors can be represented as a
  parametric sub-model plus the sum of some smooth fucntions of one or
  more of the predictor variables. Thus the model is quite flexible
  relative to strictly parametric linear or generalized linear models,
  but still has much more structure than the completely general model
  that says that the response is just some smooth function of all the
  covariates.

  Note one important point. In order for the model to be identifiable
  the smooth functions have to be constrained to have zero mean (usually
  taken over the set of covariate values). Such constraints are always
  applied by \code{gam}.

  Specification of the distribution and link function is done using the \code{\link{family}}
  argument to \code{\link{gam}} and works in the same way as for
  \code{\link{glm}}. This page hence assumes the default identity link
  normal error family, since the generalizations are easy.

  Starting with the model given above then, the \code{gam} formula would
  be \cr
  \code{y~x1+x2+s(x3)+s(x4,x5)}. \cr
  This would use the default basis for the smooths (a thin plate
  regression spline basis for each), with automatic selection of the
  effective degrees of freedom for both smooths. The dimension of the
  smoothing basis is given a default value as well (the dimension of the
  basis sets an upper limit on the maximum possible degrees of
  freedom for the basis - the limit is typically one less than basis
  dimension). Full details of how to control smooths are given in
  \code{\link{s}}. For the moment suppose that we would like to change
  the basis of the first smooth to a cubic regression spline basis with
  a dimension of 20, while fixing the second term at 25 degrees of
  freedom. The appropriate formula would be:\cr
  \code{y~x1+x2+s(x3,bs="cr",k=20)+s(x4,x5,k=26,fx=T)}.


  Now consider some more unusual models. Consider a model in which
  \eqn{y}{y} is a smooth function of \eqn{x}{x} except at a point
  \eqn{x^*}{x*} where the function jumps discontinuously. This model can
  be written as:
  \deqn{E(y_i)= \beta_0+\beta_1h(x*,x_i)+f_1(x_i)}{E(y) = b0 + b1
    h(x*,x) + f1(x)}
  where \eqn{h} is a step function jumping from 0 to 1 at
  \eqn{x^*}{x*}. The way to fit this model is to create a variable
  \code{h} which is zero for all \code{x} less than \eqn{x^*}{x*} and
  one otherwise. Then the model formula is:\cr
  \code{y~h+s(x)}. \cr

  Another situation that occurs quite often is the one in which we
  would like to find out if the model:
  \deqn{E(y_i) = f(x_i,z_i)}{E(y)=f(x,z)}
  is really necessary or whether:
  \deqn{E(y_i)=f_1(x_i)+f_2(z_i)}{E(y)=f1(x)+f2(z)}
  wouldn't do just as well. One way to do this is to look at the results
  of fitting:\cr
  \code{y~s(x)+s(z)+s(x,z)}.\cr
  \code{gam} automatically generates side conditions to make this model
  identifiable. You can also estimate "overlapping" models like:\cr
  \code{y~s(x,z)+s(z,v)}.
  
  Sometimes models of the form:
\deqn{E(y_i) = \beta_0+ f(x_i)z_i}{E(y)=b0+f(x)z}
  need to be estimated (where \eqn{f} is a smooth function, as usual.)
  The appropriate formula is:\cr
  \code{y~z+s(x,by=z)}\cr
  - the \code{by} argument ensures that the smooth function gets multiplied by
  covariate \code{z}, but GAM smooths are centred (average value zero),
  so the \code{z+} term is needed as well (\eqn{f} is being
  represented by a constant plus a centred smooth). If we'd wanted:
  \deqn{E(y_i) = f(x_i)z_i}{E(y)=f(x)z}
  then the appropriate formula would be:\cr
   \code{y~z+s(x,by=z)-1}.\cr
  
  The \code{by} mechanism also allows models to be estimated in which
  the form of a smooth depends on the level of a factor, but to do this
  the user must generate the dummy variables for each level of the
  factor. Suppose for example that \code{fac} is a factor with 3 levels
  \code{1}, \code{2}, \code{3},
  and at each level of this factor ther response depends smoothly on a
  variable \code{x} in a manner that is level dependent. Three dummy
  variables \code{fac.1}, \code{fac.2}, \code{fac.3}, can be generated for the factor   
  (e.g. \code{fac.1<-as.numeric(fac==1)}). Then the model formula would
  be:\cr
  \code{y~fac+s(x,by=fac.1)+s(x,by=fac.2)+s(x,by=fac.3)}.\cr

}

%- maybe also `usage' for other objects documented here.
\section{WARNING}{
There are no identifiability checks made between the smooth and
parametric parts of a gam formula.
}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\examples{
set.seed(10)
n<-400
sig2<-4
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
pi <- asin(1) * 2
f1 <- 2 * sin(pi * x2)
f2 <-  exp(2 * x2) - 3.75887
f3 <-  0.2 * x2^11 * (10 * (1 - x2))^6 + 
       10 * (10 * x2)^3 * (1 - x2)^10 - 1.396
e <- rnorm(n, 0, sqrt(abs(sig2)))
# A continuous `by' variable example.... 
y <- f3*x1 + e
b<-gam(y~x1-1+s(x2,by=x1))
plot(b,pages=1)
summary(b)
# A dummy `by' variable example (with a spurious covariate x0)
fac<-as.factor(c(rep(1,100),rep(2,100),rep(3,200)))
fac.1<-as.numeric(fac==1);fac.2<-as.numeric(fac==2);
fac.3<-as.numeric(fac==3)
y<-f1*fac.1+f2*fac.2+f3*fac.3+ e
b<-gam(y~fac+s(x2,by=fac.1)+s(x2,by=fac.2)+s(x2,by=fac.3)+s(x0))
plot(b,pages=1)
summary(b)
}
\keyword{models} \keyword{regression}%-- one or more ..



\eof
\name{gam.neg.bin}
\alias{gam.neg.bin}
\alias{mgcv.find.theta}
\alias{mgcv.get.scale}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{GAMs with the negative binomial distribution}
\description{The \code{gam} modelling function is designed to be able to use 
the \code{negative.binomial} and \code{neg.bin} families from the MASS library, 
with or without a known \eqn{\theta}{theta} parameter. A value for \code{theta} 
must always be passed to these families, but if \eqn{\theta}{theta} is to be
estimated then the passed value is treated as a starting value for estimation. 


If the \code{scale} argument passed to \code{gam} is positive, then it is used 
as the scale parameter \code{theta} is treated as a fixed known parameter and 
any smoothing parameters are chosen by UBRE. If \code{scale} is not positive then
\eqn{\theta}{theta} is estimated. The method of estimation is to choose \eqn{\hat \theta}{theta} 
so that the GCV (Pearson) estimate of the scale parameter is one (since the scale parameter 
is one for the negative binomial). 

\eqn{\theta}{theta} estimation is nested within the IRLS loop used for GAM fitting. After
each call to fit an iteratively weighted additive model to the IRLS pseudodata the \eqn{\theta}{theta}
estimate is updated. This is done by conditioning on all components of the current GCV/Pearson 
estimator of the scale parameter except \eqn{\theta}{theta} and then searching for the 
\eqn{\hat \theta}{theta} which equates this conditional  estimator to one. The search is 
a simple bisection search after an initial crude line search to bracket one. The search will 
terminate at the upper boundary of the search region is a Poisson fit would have yielded an estiamted 
scale parameter <1. Search limits can be set in \code{gam.control}.

Note that 
\code{neg.bin} only allows a log link, while \code{negative.binomial} also allows \code{"sqrt"} and
\code{"identity"}. In addition the \code{negative.binomial} family results in a more 
informative \code{gam} summary.  


}

%- maybe also `usage' for other objects documented here.

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\examples{
library(MASS) # required for negative binomial families
set.seed(3)
n<-400
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
pi <- asin(1) * 2
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
g<-exp(f/5)
# negative binomial data  
y<-rnbinom(g,size=3,mu=g)
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=negative.binomial(1)) # unknown theta
plot(b,pages=1)
print(b)
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=neg.bin(1)) # unknown theta
plot(b,pages=1)
print(b)
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=negative.binomial(3),scale=1) # known theta example
plot(b,pages=1)
print(b)
# Now use "sqrt" link available in negative.binomial (but not neg.bin)
f<-f-min(f);g<-f^2
y<-rnbinom(g,size=3,mu=g)
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=negative.binomial(1,link="sqrt")) 
plot(b,pages=1)
print(b)
}
\keyword{models} \keyword{regression}%-- one or more ..



\eof
\name{gam.parser}
\alias{gam.parser}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Additive Model fitting using penalized regression splines and 
GCV}
\description{ This is an internal function of package \code{mgcv}. It is a service routine for
 \code{gam} which splits of the strictly parametric part of the model formula and returns it as a formula and interprets the smooth parts 
of the model formula.  
 The processing of 
 the smooth terms results in the return of the basis dimension for each smooth term, whether 
 or not the smooth term is to be treated as a pure regression spline
 (i.e. one without a penalty) and the type of basis used to represent
 the spline. It also returns an array of covariate names which are arguments of smooths, and the name of the response variable.  
 For further information on usuage see code for \code{gam}.
}
}



\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}


\seealso{   \code{\link{gam}} \code{\link{mgcv}}}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..




\eof
\name{gam.selection}
\alias{gam.selection}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Additive Model Selection}
\description{ This page is intended to provide some more information on
  how to select GAMs. Given a model structure specified by a gam model formula,
\code{gam()} attempts to find the appropriate smoothness for each applicable model 
term using Generalized Cross Validation (GCV) or an Un-Biased Risk Estimator (UBRE),
the latter being used in cases in which the scale parameter is assumed known. GCV and 
UBRE are covered in Craven and Wahba (1979) and Wahba (1990). Fit method \code{"magic"}
uses Newton or failing that steepest descent updates of the smoothing parameters and is particularly numerically robust.
Fit method \code{"mgcv"} alternates 
grid searches for the correct overall level of smoothness for the whole model, given the 
relative smoothness of terms, with Newton/Steepest descent updates of the relative smoothness 
of terms, given the overall amount of smoothness. 

Automatic smoothness selection is unlikely to be successful with few data, particularly with
multiple terms to be selected. The methods can also fail to find the real minimum of the GCV/UBRE
score if the model contains many smooth terms that should really be completely smooth, or close to it (e.g. a straight line 
for a default 1-d smooth). The problem is that in this circumstance the optimal overall smoothness given the relative
smoothness of terms may make all terms completely smooth - but this will tend to move the smoothing parameters 
to a location where the GCV/UBRE score is nearly completely flat with respect to the smoothing parameters so that Newton and steepest  descent are both ineffective. These problems can usually be overcome by replacing some completely smooth terms with purely
parametric model terms. 

A good example of where smoothing parameter selection can ``fail'', but in an unimportant manner is provided by the 
\code{rock.gam} example in Venables and Ripley. In this case 3 smoothing parameters are to estimated from 48 data, which is 
probably over-ambitious. \code{gam} will estimate either 1.4 or 1 degrees of freedom for the smooth of \code{shape}, depending on 
the exact details of model specification (e.g. k value for each \code{s()} term). The lower GCV score is really at 1.4 (and if the other
2 terms are replaced by straight lines this estimate is always returned), but the \code{shape} term is in no way significant and the 
lowest GCV score is obtained by removing it altogether. The problem here is that the GCV score contains very little information on the optimal 
degrees of freedom to associate with a term that GCV would suggest should really be dropped. 


In general the most logically consistent method to use for deciding which terms to include in the model is to compare GCV/UBRE scores
for models with and without the term. More generally the score for the model with a smooth term can be compared to the score for the model with 
the smooth term replaced by appropriate parametric terms. Candidates for removal can be identified by reference to the approximate p-values provided by \code{summary.gam}. Candidates for replacement by parametric terms are smooth terms with estimated degrees of freedom close to their 
minimum possible.

}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\references{
Craven and Wahba (1979) Smoothing Noisy Data with Spline Functions. Numer. Math. 31:377-403

Venables and Ripley (1999) Modern Applied Statistics with S-PLUS

Wahba (1990) Spline Models of Observational Data. SIAM.

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

\url{http://www.stats.gla.ac.uk/~simon/}
}

\keyword{models} \keyword{regression}%-- one or more ..



\eof
\name{gam.setup}
\alias{gam.setup}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Additive Model set up.}
\description{ This is an internal function of package \code{mgcv}. It is called 
by \code{gam} to obtain the design matrix and penalty matrices for a GAM
set up using penalized regression splines. This is done by calling a mixture of
R routines and compiled C code.  For further information on usuage see code for \code{gam}.
}
}



\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}

}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}



\seealso{   \code{\link{gam}}, \code{\link{mgcv}}}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..









\eof
\name{gam.side.conditions}
\alias{gam.side.conditions}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Identifiability side conditions for a GAM.}
\description{ GAM formulae with repeated variables only correspond to
  identifiable models given some side conditions on the null space of
  the model penalties. This routine works out appropriate side conditions
  and returns the coefficient matrix that imposes them. It is called
  from \code{\link{gam}} and is not intended to be called by users.
}
\usage{
gam.side.conditions(G)
}
%- maybe also `usage' for other objects documented here.
\arguments{
  \item{G}{ the object obtained by calling \code{\link{gam.setup}} followed by
    \code{\link{GAMsetup}} }
}
\details{ Part of the basis of a thin plate regression spline (of the
  sort used to define GAMs in package \code{mgcv}) is a polynomial basis
  for the space of functions of zero wiggliness according to the spline
  wiggliness penalty. A GAM with repeated dependence on one or more
  variables will have multiple copies of these basis terms, and will
  hence not be identifiable without extra constraints. For example, the
  design matrix for the model 
  \code{y~s(x)+s(z)+s(x,z)} will feature two copies of the columns for
  \code{x} and \code{z}, but by imposing the constraint that the second
  parameter for each repeat variable is zero, the model becomes
  identifiable again.

  This routine automates the process of producing these constraints in
  general. The method used is to create labels for each null basis
  column that uniquely identify the term while allowing multiple copies
  of the column to be identified. The algorithm works upwards throught the
  dimensions - i.e. starts with 1-d terms, then does 2-d terms , then
  3-d and so on. A stack of null basis columns is maintained, and each
  new basis term tested against it - if the new term is a repeat then
  its parameter is constrained to be zero - otherwise it is added to
  the stack. The method is symbolic and not numerical - you
  can fool it by making two copies of the same covariate with different names.

  If the model can not be made identifiable then an error is generated.

  The code critically depends on equivalence between the R code
  null.space.basis.powers and the equivalent C code in tprs.c.

  Note that the routine does nothing to the constant terms in the basis
  - these are already dealt with by the constraints that centre all
  smooths. 

}
\value{
A matrix \eqn{ \bf C}{C} defining the constraints on the model. The constraints
are that \eqn{ { \bf Cp}={ \bf 0}}{Cp=0} where \eqn{ \bf p}{p} is the
parameter vector. \code{C} includes the original sum to zero constraints from
  \code{G$C} and is hence suitable for directly replacing that matrix. 
}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\references{
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

\url{http://www.stats.gla.ac.uk/~simon/}


}


\examples{
set.seed(0)
n<-400
sig2<-4
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
pi <- asin(1) * 2
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, 0, sqrt(abs(sig2)))
y <- f + e
b<-gam(y~s(x0)+s(x1)+s(x0,x1)+s(x2))
plot(b,pages=1)
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<-500
old.par<-par(mfrow=c(2,2))
x<-runif(n);z<-runif(n);
y<-test1(x,z)+rnorm(n)*0.1
b<-gam(y~s(x)+s(z)+s(x,z))
plot(b)
par(old.par)
rm(list=c("f","x0","x1","x2","x","z","y","b","test1","n","sig2","pi","e"))
}
\keyword{models} \keyword{regression}%-- one or more ..



\eof
\name{magic}
\alias{magic}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{ Stable Multiple Smoothing Parameter Estimation by GCV or UBRE, 
with optional fixed penalty}
\description{
Function to efficiently estimate smoothing parameters in generalized
ridge regression problems with multiple (quadratic) penalties, by GCV 
or UBRE. The function uses Newton's method in multi-dimensions, backed up by 
steepest descent to iteratively adjust the smoothing parameters for each penalty 
(one penalty may have a smoothing parameter fixed at unity ). 

For maximal numerical stability the method is based on orthogonal decomposition methods, 
and attempts to deal with numerical rank deficiency gracefully using a truncated singular 
value decomposition approach.
} %- end description

\usage{
magic(y,X,sp,S,off,rank=NULL,H=NULL,C=NULL,w=NULL,
      gamma=1,scale=1,gcv=TRUE,ridge.parameter=NULL,
      control=list(maxit=50,tol=1e-6,step.half=25,
              rank.tol=.Machine$double.eps^0.5))
}
%- maybe also `usage' for other objects documented here.
\arguments{

\item{y}{is the response data vector.}

\item{X}{is the model matrix.}

\item{sp}{is the array of smoothing parameters multiplying the penalty matrices stored in 
\code{S}. Any that are negative are autoinitialized, otherwise they are taken as supplying 
starting values. A supplied starting value will be reset to a default starting value if the
gradient of the GCV/UBRE score is too small at the supplied value.  }
 
\item{S}{ is a list of of penalty matrices. \code{S[[i]]} is the ith penalty matrix, but note
that it is not stored as a full matrix, but rather as the smallest square matrix including all 
the non-zero elements of the penalty matrix. Element 1,1 of \code{S[[i]]}  occupies 
element \code{off[i]}, \code{off[i]} of the ith penalty matrix. Each \code{S[[i]]} must be 
positive semi-definite.  }

\item{off}{is an array indicating the first parameter in the parameter vector that is 
penalized by the penalty involving \code{S[[i]]}.}

\item{rank}{ is an array specifying the ranks of the penalties. This is useful, but not 
essential, for forming square roots of the penalty matrices.}

\item{H}{ is the optional offset penalty - i.e. a penalty with a smoothing parameter fixed at 
1. This is useful for allowing regularization of the estimation process, fixed smoothing 
penalties etc.}

 \item{C}{ is the optional matrix specifying any linear equality constraints on the fitting 
problem. If \eqn{\bf b}{b} is the parameter vector then the parameters are forced to satisfy 
\eqn{ {\bf Cb} = {\bf 0} }{Cb=0}. }

\item{w}{ the regression weights. If this is a matrix then it is taken as being the 
square root of the inverse of the covariance matrix of \code{y}. If \code{w} is an array then 
it is taken as the diagonal of this matrix, or simply the weight for each element of 
\code{y}.}

\item{gamma}{is an inflation factor for the model degrees of freedom in the GCV or UBRE 
score.}

\item{scale}{ is the scale parameter for use with UBRE.}

\item{gcv}{ should be set to \code{TRUE} if GCV is to be used, \code{FALSE} for UBRE.}

\item{ridge.parameter}{It is sometimes useful to apply a ridge penalty to the fitting problem, 
penalizing the parameters in the constrained space directly. Setting this parameter to a value 
greater than zero will cause such a penalty to be used, with the magnitude given by the 
parameter value.}

\item{control}{ is a list of iteration control constants with the following elements:

\item{maxit}{The maximum number of iterations of the magic algorithm to allow.}

\item{tol}{The tolerance to use in judging convergence.}

\item{step.half}{If a trial step fails then the method tries halving it up to a maximum of 
\code{step.half} times.}

\item{rank.tol}{is a constant used to test for numerical rank deficiency of the problem. 
Basically any singular value less than \code{rank_tol} multiplied by the largest singular value of 
the  problem is set to zero.}

} %- end of control
}

\details{ 

The method is a computationally efficient means of applying GCV or UBRE (often approximately 
AIC) to the 
problem of smoothing parameter selection in generalized ridge regression problems 
of the form:
\deqn{ minimise~ \| { \bf W} ({ \bf Xb - y} ) \|^2 + {\bf b}^\prime {\bf Hb}  \sum_{i=1}^m
\theta_i {\bf b^\prime S}_i{\bf b} }{ minimise || W (Xb-y) ||^2 + b'Hb+
theta_1 b'S_1 b + theta_2 b'S_2 b + . . .}
possibly subject to constraints \eqn{ {\bf Cb}={\bf 0}}{Cb=0}. 
\eqn{ {\bf X}}{X} is a design matrix, \eqn{\bf b}{b} a parameter vector, 
\eqn{\bf y}{y} a data vector, \eqn{\bf W}{W} a weight matrix,
\eqn{ {\bf S}_i}{S_i} a positive semi-definite matrix  of coefficients
defining the ith penalty with associated smoothing parameter \eqn{\theta_i}{theta_i}, 
\eqn{\bf H}{H} is the positive semi-definite offset penalty matrix  and \eqn{\bf C}{C} a 
matrix of coefficients defining any linear equality constraints on the problem. 
\eqn{ {\bf X}}{X} need not be of full column rank.

The \eqn{\theta_i}{theta_i} are chosen to minize either the GCV score:

\deqn{V_g = \frac{n\|{\bf W}({\bf y} - {\bf Ay})\|^2}{[tr({\bf I} - \gamma {\bf A})]^2}}{V_g = n ||W(y-Ay)||^2/[tr(I - g A)]^2}

or the UBRE score:

\deqn{V_u=\|{\bf W}({\bf y}-{\bf Ay})\|^2/n-2 \phi tr({\bf I}-\gamma {\bf A})/n + \phi}{
V_u =||W(y-Ay||^2/n - 2 s tr(I - g A)/n + s }

where \eqn{\gamma}{g} is \code{gamma} the inflation factor for degrees of freedom (usually set to 1) and \eqn{\phi}{s} 
is \code{scale}, the scale parameter. \eqn{\bf A}{A} is the hat matrix (influence matrix) for the fitting problem (i.e
the matrix mapping data to fitted values). Dependence of the scores on the smoothing parameters is through \eqn{\bf A}{A}. 

The method operates by  Newton or steepest descent updates of the logs of the 
\eqn{\theta_i}{theta_i}. A key aspect of the method is stable and economical calculation of the 
first and second derivatives of the scores w.r.t. the log smoothing parameters. 
Because the GCV/UBRE scores are flat w.r.t. very large or very small \eqn{\theta_i}{theta_i}, 
it's important to get good starting parameters, and to be careful not to step into a flat region
of the smoothing parameter space. For this reason the algorithm rescales any Newton step that 
would result in a \eqn{log(\theta_i)}{log(theta_i)} change of more than 5. Newton steps are 
only used if the Hessian of the GCV/UBRE is postive definite, otherwise steepest descent is 
used. Similarly steepest descent is used if the Newton step has to be contracted too far 
(indicating that the quadratic model underlying Newton is poor). All initial steepest descent 
steps are scaled so that their largest component is 1. However a step is calculated, 
it is never expanded if it is successful (to avoid flat portions of the objective), 
but steps are successively halved if they do not decrease the GCV/UBRE score, until 
they do, or the direction is deemed to have failed. (Given the smoothing parameters the optimal 
\eqn{\bf b}{b} parameters are easily found.)


The method is coded in \code{C} with matrix factorizations performed using LINPACK and LAPACK routines.
}
\value{The function returns a list with the following items:

\item{b}{The best fit parameters given the estimated smoothing parameters.}

\item{scale}{the estimated (GCV) or supplied (UBRE) scale parameter.}

\item{score}{the minimized GCV or UBRE score.}

\item{sp}{an array of the estimated smoothing parameters.}

\item{rV}{a factored form of the parameter covariance matrix. The (Bayesian)  covariance
matrix of the parametes \code{b} is given by \code{rV\%*\%t(rV)*scale}. }

\item{gcv.info}{is a list of information about the performance of the method with the following elements:
\item{full.rank}{The apparent rank of the problem: number of parameters less number of equality constraints.}
\item{rank}{The estimated actual rank of the problem (at the final iteration of the method).}
\item{fully.converged}{is \code{TRUE} if the method converged by satisfying the convergence criteria, and \code{FALSE} if it coverged 
by failing to decrease the score along the search direction.}
\item{hess.pos.def}{is \code{TRUE} if the hessian of the UBRE or GCV score was positive definite at convergence.}
\item{iter}{is the number of Newton/Steepest descent iterations taken.}
\item{score.calls}{is the number of times that the GCV/UBRE score had to be evaluated.}
\item{rms.grad}{is the root mean square of the gradient of the UBRE/GCV score w.r.t. the smoothing parameters.}
}
Note that some further useful quantities can be obtained using \code{\link{magic.post.proc}}.

}
\references{
Wood, S.N. (2000)  Modelling and Smoothing Parameter Estimation
with Multiple  Quadratic Penalties. J.R.Statist.Soc.B 62(2):413-428

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}


\seealso{  
\code{\link{magic.post.proc}},
\code{\link{mgcv}},
\code{\link{gam}},
}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..






\eof
\name{magic.post.proc}
\alias{magic.post.proc}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Auxilliary information from magic fit}
\description{Obtains parameter covariance matrix, estimated degrees of 
freedom for each parameter and leading diagonal of influence/hat matrix, 
for a penalized regression estimated by \code{magic}.
}
\usage{
magic.post.proc(X,object,w)
}
%- maybe also `usage' for other objects documented here.
\arguments{
 \item{X}{ is the model matrix.}
\item{object}{is the list returned by \code{magic} after fitting the 
model with model matrix \code{X}.}
\item{w}{is the weight vector used in fitting, or the weight matrix used 
in fitting (i.e. supplied to 
\code{magic}, if one was.)}
}
\details{ \code{object} contains \code{rV} (\eqn{ {\bf V}}{V}, say), and 
\code{scale} (\eqn{ \phi}{s}, say) which can be 
used to obtain the require quantities as follows. The covariance matrix of 
the parameters is \eqn{ {\bf VV}^\prime \phi}{VV's}. The vector of 
estimated degrees of freedom for each parameter is the leading diagonal of 
\eqn{ {\bf VV}^\prime {\bf X}^\prime {\bf W}^\prime {\bf W}{\bf X}}{ VV'X'W'WX} 
where \eqn{\bf{W}}{W} is either the 
weight matrix \code{w} or the matrix \code{diag(w)}. The 
hat/influence  matrix is given by 
\eqn{ {\bf WX}{\bf VV}^\prime {\bf X}^\prime {\bf W}^\prime }{ WXVV'X'W'} 
.
}
\value{ A list with three items:
\item{Vb}{the covariance matrix of the model parameters.}
\item{hat}{the leading diagonal of the hat (influence) matrix.}
\item{edf}{the array giving the estimated degrees of freedom associated 
with each parameter.}
}
 
\seealso{\code{\link{magic}}}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..










\eof
\name{mgcv}
\alias{mgcv}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{ Multiple Smoothing Parameter Estimation by GCV or UBRE}
\description{
Function to efficiently estimate smoothing parameters in Generalized
Ridge Regression Problem with multiple (quadratic) penalties, by GCV 
or UBRE. The function uses Newton's method in multi-dimensions, backed up by steepest descent to iteratively 
adjust a set of relative smoothing parameters for each penalty. To ensure that the overall level of smoothing
is optimal, and to guard against trapping by local minima, a highly efficient global minimisation with respect to 
one overall smoothing parameter is also made at each iteration.

The other functions in the \code{mgcv} package are listed under `See Also'.
}
\usage{
mgcv(M)
}
%- maybe also `usage' for other objects documented here.
\arguments{
  \item{M}{ is the single argument to \code{mgcv} and should be a list
    with the elements listed below:
\itemize{ 
\item{M\$y}{The response data vector.}

\item{M\$w}{A vector of weights for the data (often proportional to the 
           reciprocal of the standard deviation of \code{M$y}). }
\item{M\$X}{The design matrix for the problem, note that \code{ncol(M$X)}
            must give the number of model parameters, while \code{nrow(M$X)} 
            should give the number of data.}
\item{M\$C}{Matrix containing any linear equality constraints 
            on the problem (i.e. \eqn{\bf C}{C} in \eqn{ {\bf Cp}={\bf 0} }{Cp=0}).}
\item{M\$S}{ A one dimensional array containing the non-zero elements of the
   penalty matrices. Let \code{start<-sum(M$df[1:(k-1)]^2)} or
   \code{start<-0} if \code{k==1}. Then penalty matrix \code{k} has
   \code{M$S[start+i+M$df[i]*(j-1)} on its ith row and jth column.
   The kth full penalty matrix is constructed by placing the matrix just
   derived into an appropriate matrix of zeroes with it's 1,1 element at
   \code{M$off[k],M$off[k]}
 }

\item{M\$m}{The number of penalty terms that are to be applied.}

\item{M\$off}{ Offset values locating the elements of \code{M$S[]} in the correct location within each penalty coefficient matrix. indexing starts at zero.}

\item{M\$fix}{ \code{M$fix[i]} is FALSE if the corresponding smoothing parameter 
is to be estimated, or TRUE if it is to be fixed at zero.}

\item{M\$df}{ \code{M$df[i]} gives the number of rows and columns of
               \code{M$S[i]} that are non-zero.}
\item{M\$sig2}{ This serves 2 purposes. If it is strictly positive then UBRE is 
used with sig2 taken as the known error variance. If it is 
non-positive then GCV is used (and on return it will contain an estimate of the error variance). 
Note that it is assumed that the variance of \eqn{y_i}{y_i} is 
given
                by \eqn{\sigma^2/w_i}{\code{sig2}/w_i}.}   

\item{M\$sp}{ An array of initial smoothing parameter estimates if any of 
these is not strictly positive then automatic generation of initial values
is used. If \code{M$fixed.sp} is non-zero then these parameters will used as known smoothing parameters.}

\item{M\$fixed.sp}{Non-zero to indicate that all smoothing parameters are supplied as known, and are not 
      to be estimated.}

\item{M\$conv.tol}{The convergence tolerence to use for the multiple
  smoothing parameter search. 1e-6 is usually ok.}

\item{M\$max.half}{Each step of the GCV minimisation is repeatedly halved if it fails to
  improve the GCV score. This is the maximum number of halvings to try
  at each step. 15 is usually enough.}

\item{M\$min.edf}{The minimum estimated degrees of freedom for the model: set to less than or equal to 0 to ignore.
                 Supplying this quantity helps to maintain numerical stability in cases in which the best model is 
                 actually very smooth. }

\item{M\$target.edf}{Set this to negative for it to be ignored. If it is positive then this is taken as the target edf
to use for the model - if the GCV/UBRE function has multiple minima w.r.t. the overall smoothing parameter then the 
minimum closest to the target edf will be selected. This feature is primarily useful when selecting smoothing parameters as part 
of iterative least squares methods - the GCV/UBRE score can develop spurious minima in these cases that are 
not representative of the problem being approximated.
}
} % end itemize
} % end M
}
\details{ 

This is documentation for the code implementing the method described in section 
4 of 
Wood (2000) . The method is a computationally efficient means of applying GCV to 
the 
problem of smoothing parameter selection in generalized ridge regression problems 
of 
the form:
\deqn{ minimise~ \| { \bf W} ({ \bf Xp - y} ) \|^2 \rho +  \sum_{i=1}^m
\lambda_i {\bf p^\prime S}_i{\bf p} }{ minimise || W (Xp-y) ||^2 rho + 
lambda_1 p'S_1 p + lambda_1 p'S_2 p + . . .}
possibly subject to constraints \eqn{ {\bf Cp}={\bf 0}}{Cp=0}. 
\eqn{ {\bf X}}{X} is a design matrix, \eqn{\bf p}{p} a parameter vector, 
\eqn{\bf y}{y} a data vector, \eqn{\bf W}{W} a diagonal weight matrix,
\eqn{ {\bf S}_i}{S_i} a positive semi-definite matrix  of coefficients
defining the ith penalty and \eqn{\bf C}{C} a matrix of coefficients 
defining any linear equality constraints on the problem. The smoothing
parameters are the \eqn{\lambda_i}{lambda_i} but there is an overall
smoothing parameter \eqn{\rho}{rho} as well. Note that \eqn{ {\bf X}}{X}
must be of full column rank, at least when projected  into the null space
of any equality constraints.  

The method operates by alternating very efficient direct searches for 
\eqn{\rho}{rho}
with Newton or steepest descent updates of the logs of the \eqn{\lambda_i}{lambda_i}. 
Because the GCV/UBRE scores are flat w.r.t. very large or very small \eqn{\lambda_i}{lambda_i}, 
it's important to get good starting parameters, and to be careful not to step into a flat region
of the smoothing parameter space. For this reason the algorithm rescales any Newton step that 
would result in a \eqn{log(\lambda_i)}{log(lambda_i)} change of more than 5. Newton steps are only used
if the Hessian of the GCV/UBRE is postive definite, otherwise steepest descent is used. Similarly steepest 
descent is used if the Newton step has to be contracted too far (indicating that the quadratic model 
underlying Newton is poor). All initial steepest descent steps are scaled so that their largest component is
1. However a step is calculated, it is never expanded if it is successful (to avoid flat portions of the objective), 
but steps are successively halved if they do not decrease the GCV/UBRE score, until they do, or the direction is deemed to have 
failed. \code{M$conv} provides some convergence diagnostics.

The method is coded in \code{C} and is intended to be portable. It should be 
noted that seriously ill conditioned problems (i.e. with close to column rank 
deficiency in the design matrix) may cause problems, especially if weights vary 
wildly between observations.  
}
\value{
  
  The function returns a list of the same form as the input list with the 
following
  elements added/modified:
  \item{M\$p}{The best fit parameters given the estimated smoothing parameters.}

  \item{M\$sp}{The estimated smoothing parameters 
(\eqn{\lambda_i/\rho}{lambda_i/rho})}

  \item{M\$sig2}{The estimated (GCV) or supplied (UBRE) error variance.}
  \item{M\$edf}{The estimated degrees of freedom associated with each penalized term. If \eqn{ {\bf p}_i}{p_i} is the vector of parameters penalized by the ith penalty then we can write \eqn{ {\bf p}_i = {\bf P}_i {\bf y}}{p_i = P_i y}, and the edf for the ith term is \eqn{ tr({\bf XP}_i)}{tr(XP_i)}... while this is sensible in the non-overlapping penalty case (e.g. GAMs) it makes alot less sense if penalties overlap!}

  \item{M\$Vp}{Estimated covariance matrix of model parameters.}

\item{M\$gcv.ubre}{The minimized GCV/UBRE score.}

\item{M\$conv}{A list of convergence diagnostics, with the following elements:
\item{edf}{Array of whole model estimated degrees of freedom.}
\item{hat}{Array of elements from leading diagonal of `hat' matrix (`influence' matrix). Same length as \code{M$y}.}
\item{score}{Array of ubre/gcv scores at the edfs for the final set of relatinve smoothing parameters.}
\item{g}{the gradient of the GCV/UBRE score w.r.t. the smoothing parameters at termination.}
\item{h}{the second derivatives corresponding to \code{g} above - i.e. the leading diagonal of the Hessian.}
\item{e}{the eigen-values of the Hessian. These should all be non-negative!}
\item{iter}{the number of iterations taken.}
\item{in.ok}{\code{TRUE} if the second smoothing parameter guess improved the GCV/UBRE score. (Please report examples 
where this is \code{FALSE})}
\item{step.fail}{\code{TRUE} if the algorithm terminated by failing to improve the GCV/UBRE score rather than by "converging". 
Not necessarily a problem, but check the above derivative information quite carefully.}
}
}
\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\section{WARNING }{ The method may not behave well with near column rank deficient \eqn{ {\bf
X}}{X}
especially in contexts where the weights vary wildly. } 

\seealso{  
\code{\link{gam}},
\code{\link{s}},
\code{\link{plot.gam}},
\code{\link{print.gam}},
\code{\link{summary.gam}},
\code{\link{predict.gam}},
\code{\link{vis.gam}},
\code{\link{residuals.gam}},
\code{\link{gam.models}},
\code{\link{gam.selection}},
\code{\link{gam.control}}
\code{\link{gam.check}}
\code{\link{pcls}},
\code{\link{mono.con}}
\code{\link{GAMsetup}},
\code{\link{gam.parser}},
\code{\link{gam.fit}},
\code{\link{gam.setup}},
\code{\link{QT}},
\code{\link{gam.neg.bin}},
\code{\link{null.space.dimension}},
\code{\link{uniquecombs}},
\code{\link{exclude.too.far}}
}

\examples{
set.seed(0)    
n<-100 # number of observations to simulate
x <- runif(4 * n, 0, 1) # simulate covariates
x <- array(x, dim = c(4, n)) # put into array for passing to GAMsetup
pi <- asin(1) * 2  # begin simulating some data
y <- 2 * sin(pi * x[1, ])
y <- y + exp(2 * x[2, ]) - 3.75887
y <- y + 0.2 * x[3, ]^11 * (10 * (1 - x[3, ]))^6 + 10 * (10 *
        x[3, ])^3 * (1 - x[3, ])^10 - 1.396
sig2<- -1    # set magnitude of variance
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- y + e          # simulated data
w <- matrix(1, n, 1) # weight matrix
par(mfrow = c(2, 2)) # scatter plots of simulated data   
plot(x[1, ], y);plot(x[2, ], y)
plot(x[3, ], y);plot(x[4, ], y)
# create list for passing to GAMsetup
G <- list(m = 4, n = n, nsdf = 0, df = c(15, 15, 15, 15),dim=c(1,1,1,1),
          s.type=c(0,0,0,0),by=0,by.exists=c(FALSE,FALSE,FALSE,FALSE),
          p.order=c(0,0,0,0), x = x,n.knots=rep(0,4),fit.method="mgcv") 
H <- GAMsetup(G)
H$y <- y    # add data to H
H$sig2 <- sig2  # add variance (signalling GCV use in this case) to H
H$w <- w       # add weights to H
H$sp<-array(-1,H$m)
H$fix<-array(FALSE,H$m)
H$conv.tol<-1e-6;H$max.half<-15
H$min.edf<-5;H$fixed.sp<-0
H <- mgcv(H)  # select smoothing parameters and fit model 
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..






\eof
\name{mono.con}
\alias{mono.con}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Monotonicity constraints for a cubic regression spline.}
\description{ Finds linear constraints sufficient for monotonicity (and
  optionally upper and/or lower boundedness) of a cubic regression
  spline. The basis representation assumed is that given by the
  \code{gam}, \code{"cr"} basis: that is the spline has a set of knots,
  which have fixed x values, but the y values of which constitute the
  parameters of the spline. }
\usage{
mono.con(x,up=TRUE,lower=NA,upper=NA)
}
%- maybe also `usage' for other objects documented here.
\arguments{
 \item{x}{The array of knot locations.}
 \item{up}{If \code{TRUE} then the constraints imply increase, if
   \code{FALSE} then decrease. }
 \item{lower}{This specifies the lower bound on the spline unless it is
   \code{NA} in which case no lower bound is imposed.}
 \item{upper}{This specifies the upper bound on the spline unless it is
   \code{NA} in which case no upper bound is imposed.}
}
\details{ Consider the natural cubic spline passing through the points:
  \eqn{ \{x_i,p_i:i=1 \ldots n \} }{ (x_i,p_i), i=1..n}. Then it is possible
  to find a relatively small set of linear constraints on \eqn{ \bf p}{p}
  sufficient to ensure monotonicity (and bounds if required):
  \eqn{ {\bf Ap}\ge{\bf b} }{Ap>=b}. Details are given in Wood (1994).
  This function returns a list containing \code{A} and \code{b}. 
}
\value{ The function returns a list containing constraint matrix
  \code{A} and constraint vector \code{b}.   
}
\references{

Gill, P.E., Murray, W. and Wright, M.H. (1981) Practical Optimization. Academic
Press, London. 

Wood, S.N. (1994) Monotonic smoothing splines fitted by cross validation SIAM
Journal on Scientific Computing 15(5):1126-1133


\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk} } 

\seealso{  \code{\link{mgcv} } \code{\link{pcls}}  }

\examples{
# Fit a monotonic penalized regression spline .....

# Generate data from a monotonic truth.
set.seed(10);x<-runif(100)*4-1;x<-sort(x);
f<-exp(4*x)/(1+exp(4*x));y<-f+rnorm(100)*0.1;plot(x,y)
dat<-data.frame(x=x,y=y)
# Show regular spline fit (and save fitted object)
f.ug<-gam(y~s(x,k=10,bs="cr"));lines(x,fitted(f.ug))
# Create Design matrix, constriants etc. for monotonic spline....
gam.setup(y~s(x,k=10,bs="cr")-1,dat,fit.method="mgcv")->G;
GAMsetup(G)->G;F<-mono.con(G$xp);
G$Ain<-F$A;G$bin<-F$b;G$C<-matrix(0,0,0);G$sp<-f.ug$sp;
G$p<-G$xp;G$y<-y;G$w<-y*0+1

p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

# now modify the gam object from unconstrained fit a little, to use it
# for predicting and plotting constrained fit. 
p<-c(0,p);f.ug$coefficients<-p; 
lines(x,predict.gam(f.ug,newdata=data.frame(x=x)),col=2)
}
\keyword{models} \keyword{smooth} \keyword{regression} %-- one or more ..







\eof
\name{mroot}
\alias{mroot}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Smallest square root of matrix}
\description{ Find a square root of a positive semi-definite matrix, 
having as few columns as possible. Uses either pivoted choleski 
decomposition or singular value decomposition to do this.  
}
\usage{
mroot(A,rank=NULL,method="chol")
}
%- maybe also `usage' for other objects documented here.
\arguments{
 \item{A}{ The positive semi-definite matrix, a square root of which is 
to be found.}
\item{rank}{if the rank of the matrix \code{A} is known then it should 
be supplied.}
\item{method}{ \code{"chol"} to use pivoted choloeski decompositon, 
which is fast but tends to over-estimate rank. \code{"svd"} to use 
singular value decomposition, which is slow, but is the most accurate way 
to estimate rank.}
}

\details{ The routine uses an LAPACK SVD routine, or the LINPACK pivoted 
Choleski routine. It is primarily of use for turning penalized regression 
problems into ordinary regression problems.} 
\value{ A matrix, \eqn{ {\bf B}}{B} with as many columns as the rank of 
\eqn{ {\bf A}}{A}, and such that \eqn{ {\bf A} = {\bf BB}^\prime}{A=BB'}.}
 

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..










\eof
\name{null.space.dimension}
\alias{null.space.dimension}
\alias{null.space.basis.powers}
\alias{null.space.basis.labels}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{The basis of the space of un-penalized functions for a t.p.r.s.}
\description{ The thin plate spline penalties give zero penalty to some
  functions. The space of these functions is spanned by a set of
  polynomial terms. \code{null.space.dimension} finds the dimension of this space, \eqn{M}{M}, given
  the number of covariates that the smoother is a function of, \eqn{d}{d},
  and the order of the smoothing penalty, \eqn{m}{m}. If \eqn{m}{m} does not
  satisfy \eqn{2m>d}{2m>d} then the smallest possible dimension
  for the null space is found given \eqn{d}{d} and the requirement that
  the smooth should be visually smooth. \code{null.space.basis.powers()}
  produces an \code{M} by \code{d} array, the ith row of which gives the
  powers of each of the d covariates featuring in the ith basis
  term. \code{null.space.basis.names()} gives the equivalent information
  as a list of basis labels - the key feature of which is that the same
  term will always have the same label, irrespective of the order in
  which variables are supplied to the routine. For some models a smooth
  may be multiplied by a \code{by} variable: in such cases the basis
  labels can include this name as well, but note that the constant label
  will be returned as if the \code{by} variable didn't exist (to meet
  requirements of function \code{gam.side.conditions}).
 
}
\usage{
null.space.dimension(d,m)
null.space.basis.powers(m,d)
null.space.basis.labels(names,m,by="NA")
}
%- maybe also `usage' for other objects documented here.
\arguments{
 \item{d}{ is a positive integer - the number of variables of which the
   t.p.s. is a function. }
 \item{m}{ a non-negative integer giving the order of the penalty
   functional, or signalling that the default order should be used.}
 \item{names}{is an array of \code{d} variable names}
 \item{by}{the name of any \code{by} variable, which mutliplies all the
   terms in the penalty null space basis.}
}
\details{ Thin plate splines are only visually smooth if the order of the
  wiggliness penalty, \eqn{m}{m}, satisfies \eqn{2m > d+1}{2m >
    d+1}. If \eqn{2m<d+1}{2m<d+1} then this routine finds the smallest
  \eqn{m}{m} giving visual smoothness
  for the given \eqn{d}{d}, otherwise the supplied \eqn{m}{m} is used. The null space dimension is given by:

  \eqn{M=(m+d+1)!/(d!(m-d)!)}{M=(m+d+1)!/(d!(m-d)!}

  which is the value returned.

  
}
\value{
\code{null.space.dimension} returns an integer (array), the null space dimension
\eqn{M}{M}: this is the only one of these functions that accepts array arguments. \cr  \code{null.space.basis.powers} produces a 2-d array,
\code{p}, say: each row corresponds to one of the \code{M} basis vectors, while each column
corresponds to one covariate. So if \code{x1}, \code{x2}, \code{x3}, etc
are the covariates, the ith basis vector is given by
\code{x1^p[i,1]*x2^p[i,2]*x3^p[i,3]...}. \cr
\code{null.space.basis.labels()} just produces uniques labels for the
basis terms.    
 
}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

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

\url{http://www.stats.gla.ac.uk/~simon/}
}

\section{WARNING}
{ This routine operates symbolically, not numerically: your model may
  still be numerically un-identifiable, and the routine can't tell if
  two variables with different names are really the same.
}  

\seealso{\code{\link{gam.side.conditions}}, \code{\link{gam}}}

\examples{
null.space.dimension(2,0)
null.space.basis.powers(2,2)
null.space.basis.labels(c("x","z"),m=2)
null.space.basis.labels(c("z","x"),m=2) # labels are invariant
}
\keyword{models} \keyword{regression}%-- one or more ..



\eof
\name{pcls}
\alias{pcls}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{ Penalized Constrained Least Squares Fitting}
\description{
Solves least squares problems with quadratic penalties subject to linear
equality and inequality constraints using quadratic programming.
}
\usage{
pcls(M)
}
%- maybe also `usage' for other objects documented here.
\arguments{
  \item{M}{is the single list argument to \code{pcls}. It should have  the 
  following elements (last 3 are not in argument for \link{mgcv}) :
\itemize{
 \item{y}{The response data vector.}
 \item{w}{A vector of weights for the data (often proportional to the 
           reciprocal of the variance). }
 \item{X}{The design matrix for the problem, note that \code{ncol(M$X)}
            must give the number of model parameters, while \code{nrow(M$X)} 
            should give the number of data.}
\item{C}{Matrix containing any linear equality constraints 
            on the problem (i.e. \eqn{ \bf C}{C} in \eqn{ {\bf Cp}={\bf
		0} }{Cp=0}). If you have no equality constraints
	    initialize this to a zero by zero matrix.}
 \item{S}{ A one dimensional array containing the non-zero elements of the
   penalty matrices. Let \code{start<-sum(M$df[1:(k-1)]^2)} or
   \code{start<-0} if \code{k==1}. Then penalty matrix \code{k} has
   \code{M$S[start+i+M$df[i]*(j-1)} on its ith row and jth column.
   The kth full penalty matrix is constructed by placing the matrix just
   derived into an appropriate matrix of zeroes with it's 1,1 element at
   \code{M$off[k],M$off[k]}
 }
 \item{off}{ Offset values locating the elements of \code{M$S[]} in
   the correct location within each penalty coefficient matrix. (Zero
   offset implies starting in first location)}

 \item{df}{ \code{M$df[i]} gives the number of rows and columns of
               \code{M$S[i]} that are non-zero.}
  
\item{sp}{ An array of  smoothing parameter estimates.}
\item{p}{An array of feasible initial parameter estimates - these must
satisfy the constraints, but should avoid satisfying the inequality
constraints as equality constraints.}
\item{Ain}{Matrix for the inequality constraints \eqn{ {\bf A}_{in}
    {\bf p} > {\bf b}_{in}}{A_in p > b}. }
\item{bin}{vector in the inequality constraints. }
} % end itemize
} % end M
}
\details{ 

  This solves the problem:
 
\deqn{ minimise~ \| { \bf W}^{1/2} ({ \bf Xp - y} ) \|^2  +  \sum_{i=1}^m
\lambda_i {\bf p^\prime S}_i{\bf p} }{ minimise || W^0.5 (Xp-y) ||^2 + 
lambda_1 p'S_1 p + lambda_1 p'S_2 p + . . .}
subject to constraints \eqn{ {\bf Cp}={\bf 0}}{Cp=0} and \eqn{ {\bf
    A}_{in}{\bf p}>{\bf b}_{in}}{A_in p > b_in}, w.r.t. \eqn{\bf p}{p} given the
smoothing parameters \eqn{\lambda_i}{lambda_i}.
\eqn{ {\bf X}}{X} is a design matrix, \eqn{\bf p}{p} a parameter vector, 
\eqn{\bf y}{y} a data vector, \eqn{\bf W}{W} a diagonal weight matrix,
\eqn{ {\bf S}_i}{S_i} a positive semi-definite matrix  of coefficients
defining the ith penalty and \eqn{\bf C}{C} a matrix of coefficients 
defining the linear equality constraints on the problem. The smoothing
parameters are the \eqn{\lambda_i}{lambda_i}. Note that \eqn{ {\bf X}}{X}
must be of full column rank, at least when projected  into the null space
of any equality constraints. \eqn{ {\bf A}_{in}}{A_in} is a matrix of
coefficients defining the inequality constraints, while \eqn{ {\bf
    b}_{in}}{b_in} is a vector involved in defining the inequality constraints.  

Quadratic programming is used to perform the solution. The method used
is designed for maximum stability with least squares problems:
i.e. \eqn{ {\bf X}^\prime {\bf X}}{X'X} is not formed explicitly. See
Gill et al. 1981.

}
\value{ The function returns an array containing the estimated parameter
  vector. 
   
}
\references{

Gill, P.E., Murray, W. and Wright, M.H. (1981) Practical Optimization. Academic
Press, London. 

Wood, S.N. (1994) Monotonic smoothing splines fitted by cross validation SIAM
Journal on Scientific Computing 15(5):1126-1133

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}} 

\seealso{  \code{\link{mgcv}} \code{\link{mono.con}}  }

\examples{
# first an un-penalized example - fit E(y)=a+bx subject to a>0
set.seed(0)
n<-100
x<-runif(n);y<-x-0.2+rnorm(n)*0.1
M<-list(X=matrix(0,n,2),p=c(0.1,0.5),off=array(0,0),df=array(0,0),S=0,
Ain=matrix(0,1,2),bin=0,C=matrix(0,0,0),sp<-0,y=y,w=y*0+1)
M$X[,1]<-1;M$X[,2]<-x;M$Ain[1,]<-c(1,0)
pcls(M)->M$p
plot(x,y);abline(M$p,col=2);abline(coef(lm(y~x)),col=3)

# and now a penalized example: a monotonic penalized regression spline .....

# Generate data from a monotonic truth.
x<-runif(100)*4-1;x<-sort(x);
f<-exp(4*x)/(1+exp(4*x));y<-f+rnorm(100)*0.1;plot(x,y)
dat<-data.frame(x=x,y=y)
# Show regular spline fit (and save fitted object)
f.ug<-gam(y~s(x,k=10,bs="cr"));lines(x,fitted(f.ug))
# Create Design matrix, constriants etc. for monotonic spline....
gam.setup(y~s(x,k=10,bs="cr")-1,dat,fit.method="mgcv")->G;
GAMsetup(G)->G;F<-mono.con(G$xp);
G$Ain<-F$A;G$bin<-F$b;G$C<-matrix(0,0,0);G$sp<-f.ug$sp;
G$p<-G$xp;G$y<-y;G$w<-y*0+1;

p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

# now modify the gam object from unconstrained fit a little, to use it
# for predicting and plotting constrained fit. 
p<-c(0,p);f.ug$coefficients<-p; 
lines(x,predict.gam(f.ug,newdata=data.frame(x=x)),col=2)

# now a tprs example of the same thing....

f.ug<-gam(y~s(x,k=10));lines(x,fitted(f.ug))
# Create Design matrix, constriants etc. for monotonic spline....
gam.setup(y~s(x,k=10),dat,fit.method="mgcv")->G;
GAMsetup(G)->G;
nc<-40         # number of constraints
xc<-0:nc/nc # points on [0,1]  
xc<-xc*4-1  # points at which to impose constraints
A0<-predict.gam(f.ug,data.frame(x=xc),type="lpmatrix") 
# ... A0%*%p evaluates spline at xc points
A1<-predict.gam(f.ug,data.frame(x=xc+1e-6),type="lpmatrix") 
A<-(A1-A0)/1e-6    
# ... approx. constraint matrix (A%*%p is -ve spline gradient at points xc)
G$Ain<-A;    # constraint matrix
G$bin<-rep(0,nc);  # constraint vector
G$sp<-f.ug$sp; # use smoothing parameters from un-constrained fit
G$p<-rep(0,11);G$p[11]<-0.1  
# ... monotonic start params, got by setting coefs of polynomial part
G$p[10]<- -mean(0.1*x) 
# ... must meet gam side conditions: sum of smooth over x's is zero 
G$y<-y;G$w<-y*0+1
p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

# now modify the gam object from unconstrained fit a little, to use it
# for predicting and plotting constrained fit. 
f.ug$coefficients<-p; 
lines(x,predict.gam(f.ug,newdata=data.frame(x=x)),col=3)

}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..





\eof
\name{plot.gam}
\alias{plot.gam}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Default GAM plotting}
\description{ Takes a fitted \code{gam} object produced by \code{gam()} and plots the 
   component smooth functions that make it up, on the scale of the linear predictor.}
}
\usage{
plot.gam(x,rug=TRUE,se=TRUE,pages=0,select=NULL,scale=-1,n=100,n2=40,
         pers=FALSE,theta=30,phi=30,jit=FALSE,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{ 
  
 \item{x}{ a fitted \code{gam} object as produced by \code{gam()}.
                }
 \item{rug}{ when TRUE (default) then the covariate to which the plot applies is displayed as a rug plot
 at the foot of each plot of a 1-d smooth, and the locations of the
 covariates are plotted as points on the contour plot representing a 2-d
 smooth.} 
  
\item{se}{ when TRUE (default) upper and lower lines are added to the
  1-d plots at 2 standard errors
        above and below the estimate of the smooth being plotted while for
	2-d plots, surfaces at +1 and -1 standard errors are contoured
	and overlayed on the contour plot for the estimate. If a
	positive number is supplied then this number is multiplied by
	the standard errors when calculating standard error curves or surfaces.}
 \item{pages}{ (default 0) the number of pages over which to spread the output. For example, 
if \code{pages=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 are. }
\item{select}{Allows the  plot for a single model term to be selected for printing. e.g. if you just want the plot for the second smooth term set \code{select=2}. }
\item{scale}{ set to -1 (default) to have the same y-axis scale for each plot, and to 0 for a 
          different y axis for each plot.  }
\item{n}{ number of points used for each 1-d plot - for a nice smooth plot this needs to be several times the estimated 
  degrees of freedom for the smooth. Default value 100.}
\item{n2}{Square root of number of points used to grid estimates of 2-d
  functions for contouring.}
\item{pers}{Set to \code{TRUE} if you want perspective plots for 2-d
  terms.}
\item{theta}{One of the perspective plot angles.}
\item{phi}{The other perspective plot angle.}
\item{jit}{Set to TRUE if you want rug plots for 1-d terms to be jittered.}
\item{...}{ other arguments.}
}
\details{ Produces default plot showing the smooth components of a
  fitted GAM.

  For plots of 1-d smooths, the x axis of each plot is labelled 
   with the covariate name, while the y axis is labelled \code{s(cov,edf) } where \code{cov}
   is the covariate name, and \code{edf} the estimated (or user defined for regression splines) degrees of freedom of the smooth.

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 \code{s(var1,var2,edf)}, indicating the
variables of which the term is a function, and the estimated degrees of
freedom for the term. When \code{se=TRUE}, estimator variability is shown by overlaying
contour plots at plus and minus 1 s.e. relative to the main
estimate. If \code{se} is a positive number then contour plots are at plus or minus \code{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.  
   
   Within the function, the data for the plots is obtained by direct
   calls to the compiled C code that \code{predict.gam} uses.   

   Smooths of more than 2 variables are not currently dealt with, but
   simply generate a warning, but see \code{\link{vis.gam}}.
   }

\value{ The function simply generates plots.
}

\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398


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

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\section{WARNING }{ Note that the behaviour of this function is not identical to 
  \code{plot.gam()} in Splus.

  The function can not deal with smooths of more than 2 variables!
} 

\seealso{  \code{\link{gam}}, \code{\link{predict.gam}}, \code{\link{vis.gam}}}

\examples{
library(mgcv)
set.seed(0)
n<-200
sig2<-4
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
pi <- asin(1) * 2
y <- 2 * sin(pi * x0)
y <- y + exp(2 * x1) - 3.75887
y <- y + 0.2 * x2^11 * (10 * (1 - x2))^6 + 10 * (10 * x2)^3 * (1 - x2)^10 - 1.396
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- y + e
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))
plot(b,pages=1)
# example with 2-d plots
b1<-gam(y~s(x0,x1)+s(x2)+s(x3))
op<-par(mfrow=c(2,2))
plot(b1)
par(op) 
}
\keyword{models} \keyword{smooth} \keyword{regression} \keyword{hplot}%-- one or more ...







\eof
\name{predict.gam}
\alias{predict.gam}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Prediction from fitted GAM model}
\description{ Takes a fitted \code{gam} object produced by \code{gam()} 
and produces predictions given a new set of values for the model covariates 
or the original values used for the model fit.}
}
\usage{
predict.gam(object,newdata,type="link",se.fit=FALSE,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{ 
  
 \item{object}{ a fitted \code{gam} object as produced by \code{gam()}.
                }
 \item{newdata}{ A data frame containing the values of the model covariates at which predictions
  are required. If this is not provided then predictions corresponding to the original data are returned. 
 
 } 
 
 
\item{type}{ When this has the value \code{"link"} (default) the linear predictor (possibly with
associated standard errors) is returned. When \code{type="terms"} each component of the 
linear predictor is returned seperately (possibly with standard errors): this includes 
parametric model components, followed by each smooth component, but excludes any offset. When \code{type="response"} predictions 
on the scale of the response are returned (possibly with approximate
standard errors). When \code{type="lpmatrix"} then a matrix is returned
which yields the values of the linear predictor (minus any offset) when applied to the
parameter vector (in this case \code{se.fit} is ignored). The latter
option is most useful for getting variance estimates for integrated quantities. }
 \item{se.fit}{ when this is TRUE (not default) standard error estimates are returned for each prediction.}
 \item{...}{ other arguments.}
}

\value{ If \code{type=="lpmatrix"} then a matrix is returned which will
  give a vector of linear predictor values (minus any offest) at the supplied covariate
  values, when applied to the model coefficient vector. 
Otherwise, if \code{se.fit} is \code{TRUE} then a 2 item list is returned with items (both arrays) \code{fit}
and \code{se.fit} containing predictions and associated standard error estimates, otherwise an 
array of predictions is returned. The dimensions of the returned arrays depends on whether 
\code{type} is \code{"terms"} or not: if it is then the array is 2 dimensional with each 
term in the linear predictor separate, otherwise the array is 1 dimensional and contains the 
linear predictor/predicted values (or corresponding s.e.s). The linear predictor returned termwise will 
not include the offset.
}

\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\section{WARNING }{ Note that the behaviour of this function is not identical to 
\code{predict.gam()} in Splus.} 

\seealso{  \code{\link{gam}} \code{\link{plot.gam}}}

\examples{
library(mgcv)
n<-200
sig2<-4
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
pi <- asin(1) * 2
y <- 2 * sin(pi * x0)
y <- y + exp(2 * x1) - 3.75887
y <- y + 0.2 * x2^11 * (10 * (1 - x2))^6 + 10 * (10 * x2)^3 * (1 - x2)^10 - 1.396
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- y + e
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))
newd<-data.frame(x0=(0:30)/30,x1=(0:30)/30,x2=(0:30)/30,x3=(0:30)/30)
pred<-predict.gam(b,newd)
 
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..


\eof
\name{print.gam}
\alias{print.gam}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Additive Model default print statement}
\description{ This is the default print statement for a GAM object. If you need a list of everything
that is part of a gam object see \code{\link{gam}}, or use \code{names()}. The family (including link), 
model formula, and estimated degrees of freedom for each model term (plus total) are printed, as 
well as the minimized GCV or UBRE score, depending on which was used.
}
}

\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@state.gla.ac.uk}}



\seealso{   \code{\link{gam}} \code{\link{mgcv}}}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..










\eof
\name{residuals.gam}
\alias{residuals.gam}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Additive Model residuals}
\description{Returns residuals for a fitted \code{gam} model
  object. Pearson, deviance, working and response residuals are
  available. 
}

\usage{
residuals.gam(object, type = c("deviance", "pearson","scaled.pearson", "working", "response"),...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
  \item{object}{ a \code{gam} fitted model object. }
  \item{type}{the type of residuals wanted. }
   \item{...}{other arguments.}
}
\details{Response residuals are the raw residuals (data minus fitted
  values). Scaled Pearson residuals are raw residuals divided by the standard
  deviation of the data according to the model mean variance
  relationship and estimated scale parameter. Pearson residuals are the same, but multiplied by the square root 
  of the scale parameter (so they are independent of the scale parameter):
  (\eqn{(y-\mu)/\sqrt{V(\mu)}}{(y-m)/V(m)^0.5}, where  \eqn{y}{y} is data \eqn{\mu}{m} 
  is model fitted value and \eqn{V}{V} is model mean-variance relationship.). Both are provided since not all texts 
  agree on the definition of Pearson residuals. Deviance residuals simply
  return the deviance residuals defined by the model family. Working
  residuals are the residuals returned from model fitting at convergence.

  There is  a special function for \code{gam} objects because of a bug
  in the calculation of Pearson residuals in some versions of
  residual.glm().
 
}
\value{ An array of residuals.
}


\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}



\seealso{   \code{\link{gam}} \code{\link{mgcv}}}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..













\eof
\name{s}
\alias{s}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Defining smooths in GAM formulae}
\description{ Function used in definition of smooth terms within
  \code{gam} model formulae. The function does not evaluate a (spline)
  smooth - it exists purely to help set up a model using spline based smooths.
}

\usage{s(..., k=-1,fx=FALSE,bs="tp",m=0,by=NA)}
}
%- maybe also `usage' for other objects documented here.
\arguments{
  \item{...}{ a list of variables that are the covariates that this
    smooth is a function of.}
  \item{k}{ the dimension of the basis used to represent the smooth term.
    The default depends on the number of variables that the smooth is a
    function of. \code{k} should not be less than the dimension of the
    null space of the penalty for the term (see
    \code{\link{null.space.dimension}}), but will be reset if
  it is.}
  \item{fx}{indicates whether the term is a fixed d.f. regression
  spline (\code{TRUE}) or a penalized regression spline (\code{FALSE}).}
  \item{bs}{this can be \code{"cr"} for a cubic regression spline or
    \code{"tp"} for a thin plate regression spline. Only thin plate
    regression splines can be used for multidimensional smooths, so this
    argument only has an effect for univariate smooths. Note that the
    \code{"cr"} basis is faster than the \code{"tp"} basis, particularly
  on large data sets.}
  \item{m}{The order of the penalty for this t.p.r.s. term (e.g. 2 for
    normal cubic spline penalty with 2nd derivatives). O signals
    autoinitialization, which sets the order to the lowest value
    satisfying 2m>d+1, where d is the number of covariates: this choise
    ensures visual smoothness. In addition, m must satisfy the technical
    restriction 2m>d, otherwise it will be autoinitialized.}
  \item{by}{specifies a covariate by which the whole smooth term is to
    be multiplied. This is particularly useful for creating models in
    which a smooth interacts with a factor: in this case the \code{by}
    variable would usually be the dummy variable coding one level of the
  factor. See the examples below.}
  }
}
\details{The function does not evaluate the variable arguments. It will
  correctly interpret calls like \code{s(x,14|f)} (results in pure
  regression spline smooth of one variable with a 14 knot cubic
  regression spline basis), \code{s(x,z,20)} (a penalized regression
  spline of 2 covariates using a 20 dimensional t.p.r.s. basis),
  etc. but this feature is purely for back compatibility reasons, and
  may not be maintained indefinitely. 

}
\value{ A list with the following items:
  \item{term}{An array of text strings giving the names of the covariates that 
the term is a function of.}
  \item{bs.dim}{The dimension of the basis used to represent the smooth.}
  \item{bs.type}{The type of basis. 0 is cubic regression spline. 1 is
    thin plate regression spline. 0 can only be used for 1-d smooths.}
  \item{fixed}{TRUE if the term is to be treated as a pure regression
    spline (with fixed degrees of freedom); FALSE if it is to be treated
  as a penalized regression spline}
  \item{dim}{The dimension of the smoother - i.e. the number of
    covariates that it is a function of.}
  \item{p.order}{The order of the t.p.r.s. penalty, or 0 for
    auto-selection of the penalty order.}
  \item{by}{is the name of any \code{by} variable as text ("NA" for none).}
  \item{full.call}{Text for pasting into a string to be converted to a
    gam formula, which has the values of function options given explicitly -
  this is useful for constructing a fully expanded gam formula which can
be used without needing access to any variables that may have been used
to define k, fx, bs or m in the original call. i.e. this is text which
when parsed and evaluated generates a call to \code{s()} with all the
options spelled out explicitly.}
}


\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}
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

\url{http://www.stats.gla.ac.uk/~simon/}

}

\seealso{   \code{\link{gam}}}

\examples{
# example utilising `by' variables
library(mgcv)
set.seed(0)
n<-200;sig2<-4
x1 <- runif(n, 0, 1);x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
fac<-c(rep(1,n/2),rep(2,n/2)) # create factor
fac.1<-rep(0,n)+(fac==1);fac.2<-1-fac.1 # and dummy variables
fac<-as.factor(fac)
f1 <-  exp(2 * x1) - 3.75887
f2 <-  0.2 * x1^11 * (10 * (1 - x1))^6 + 10 * (10 * x1)^3 * (1 - x1)^10 - 1.396
f<-f1*fac.1+f2*fac.2+x2
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- f + e
# NOTE: smooths will be centered, so need to include fac in model....
b<-gam(y~fac+s(x1,by=fac.1)+s(x1,by=fac.2)+x2) 
plot(b,pages=1)
}

\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..










\eof
\name{summary.gam}
\alias{summary.gam}
\alias{print.summary.gam}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Summary for a GAM fit}
\description{ Takes a fitted \code{gam} object produced by \code{gam()} and produces various useful
summaries from it. 
}
\usage{
summary.gam(object,...)
print.summary.gam(x,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{ 
\item{object}{ a fitted \code{gam} object as produced by \code{gam()}.}
\item{x}{a \code{summary.gam} object produced by \code{summary.gam()}.} 
\item{...}{ other arguments.}
}
\details{ Model degrees of freedom are taken as the trace of the influence (or hat) matrix \eqn{ {\bf A}}{A} for the model fit.
Residual degrees of freedom are taken as number of data minus model degrees of freedom. Let \eqn{ {\bf P}_i}{P_i} be the matrix 
giving the parameters of the ith smooth when applied to the data (or pseudodata in the generalized case) and let \eqn{ {\bf X}}{X} 
be the design matrix of the model. Then \eqn{ tr({\bf XP}_i )}{tr(XP_i)} is the edf for the ith term. Clearly this definition causes 
the edf's to add up properly! 

\code{print.summary.gam} tries to print various bits of summary information useful for term selection in a pretty way.
}

\value{\code{summary.gam} produces a list of summary information for a fitted \code{gam} object. 
\item{p.coeff}{is an array of estimates of the strictly parametric model coefficients.}
\item{p.t}{is an array of the \code{p.coeff}'s divided by their standard errors.}
\item{p.pv}{is an array of p-values for the null hypothesis that the corresponding parameter is zero. 
Calculated with reference to the t distribution with the estimated residual degrees of freedom for the model fit.}
\item{m}{The number of smooth terms in the model.}
\item{chi.sq}{An array of test statistics for assessing the significance of model smooth terms. If \eqn{ {\bf p}_i}{p_i} 
is the parameter vector for the ith smooth term, and this term has estimated covariance matrix \eqn{ {\bf V}_i}{V_i} then the 
statistic is \eqn{ {\bf p}_i^\prime {\bf V}_i^{k-} {\bf p}_i}{p_i'V_i^{k-}p_i}, where \eqn{ {\bf V}^{k-}_i}{V_i^{k-}} is the rank k-1 
pseudo-inverse of \eqn{ {\bf V_i}}{V_i}, and k is the basis dimension.}
\item{s.pv}{An array of approximate p-values for the null hypotheses that each smooth term is zero. Be warned, these are only 
approximate. In the case in which UBRE has been used, they are obtained by comparing the chi.sq statistic given above to the 
chi-squared distribution with degrees of  freedom given by the estimated degrees of freedom for the term. In the GCV case (in 
which the scale parameter will have been estimated) the statistic is compared to an F distribution with upper d.f. given by
the estimate degrees of freedom for the term, and lower d.f. given by the residual degrees of freedom for the model . 
Use at your own risk! Typically the p-values will be somewhat inaccurate
because they are conditional on the smoothing parameters and the distributional assumption doesn't have a firm theoretical basis. 
A pragmatic approach to the latter issue is to check p-values by refitting the model using regression splines
with each basis dimension set to one more than the rounded edf for
the term (see example, below). In this latter case the distributional
assumption would be fine if the smoothing parameters were
known, but of course they are not, and conditioning on the smoothing
parameters will always be problematic. The difficulty is that GCV has
used the data to select the most plausible smoothing parameters for each
term - not surprisingly this tends to mean that the terms appear
a little more `significant' than they should. Of course this problem
is no different to the standard difficulty in interpreting p-values for
model terms when the model has been selected by hypothesis testing
methods such as backwards elimination.   
}
\item{se}{array of standard error estimates for all parameter estimates.}
\item{r.sq}{The adjusted r-squared for the model. Defined as the proportion of variance explained, where original variance and 
residual variance are both estimated using unbiased estimators. This quantity can be negative if your model is worse than a one 
parameter constant model, and can be higher for the smaller of two nested models! Note that proportion null deviance 
explained is probably more appropriate for non-normal errors.}
\item{dev.expl}{The proportion of the null deviance explained by the model.}
\item{edf}{array of estimated degrees of freedom for the model terms.}
\item{residual.df}{estimated residual degrees of freedom.}
\item{n}{number of data.}
\item{gcv}{minimized GCV score for the model, if GCV used.}
\item{ubre}{minimized UBRE score for the model, if UBRE used.}
\item{scale}{estimated (or given) scale parameter.}
\item{family}{the family used.}
\item{formula}{the original GAM formula.}
}

\references{

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

\url{http://www.stats.gla.ac.uk/~simon/}
}
\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}

\section{WARNING }{ The supplied p-values are only approximate: they are
  conditional on the smoothing parameters, which will usually
  have been estimated.
} 

\seealso{  \code{\link{gam}}, \code{\link{predict.gam}}, \code{\link{gam.check}} }

\examples{
library(mgcv)
set.seed(0)
n<-200
sig2<-4
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
pi <- asin(1) * 2
y <- 2 * sin(pi * x0)
y <- y + exp(2 * x1) - 3.75887
y <- y + 0.2 * x2^11 * (10 * (1 - x2))^6 + 10 * (10 * x2)^3 * (1 - x2)^10 - 1.396
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- y + e
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))
plot(b,pages=1)
summary(b)
# now check the p-values by using a pure regression spline.....
b.d<-round(b$edf)+1 
b.d<-pmax(b.d,3) # can't have basis dimension less than this!
bc<-gam(y~s(x0,k=b.d[1],fx=TRUE)+s(x1,k=b.d[2],fx=TRUE)+
        s(x2,k=b.d[3],fx=TRUE)+s(x3,k=b.d[4],fx=TRUE))
plot(bc,pages=1)
summary(bc)
# p-value check - increase k to make this useful!
n<-200;p<-0;k<-20
for (i in 1:k)
{ b<-gam(y~s(x,z),data=data.frame(y=rnorm(n),x=runif(n),z=runif(n)))
  p[i]<-summary(b)$s.p[1]
}
plot(((1:k)-0.5)/k,sort(p))
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ...







\eof
\name{uniquecombs}
\alias{uniquecombs}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{find the unique rows in a matrix }
\description{
This routine returns a matrix containing all the unique rows of the
matrix supplied as its argument. That is, all the duplicate rows are
stripped out. Note that the ordering of the rows on exit is not the same
as on entry.
}
\usage{
uniquecombs(x)
}
%- maybe also `usage' for other objects documented here.
\arguments{
 \item{x}{ is an \R matrix }
}
\details{ Models with more parameters than unique combinations of
  covariates are not identifiable. This routine provides a means of
  evaluating the number of unique combinations of coavariates in a
  model. The routine calls compiled C code.
   
}
\value{
A matrix consisting of the unique rows of \code{x} (in arbitrary order).
 
}

\author{ Simon N. Wood \email{simon@stats.gla.ac.uk}}


\examples{
X<-matrix(c(1,2,3,1,2,3,4,5,6,1,3,2,4,5,6,1,1,1),6,3,byrow=TRUE)
print(X)
uniquecombs(X)
}
\keyword{models} \keyword{regression}%-- one or more ..



\eof
\name{vis.gam}
\alias{vis.gam}
\alias{persp.gam}
\title{Visualization of GAM objects}
\usage{
vis.gam(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,
        color="topo",se=-1,type="link",zlim=NULL,...)
}

\arguments{
  \item{x}{a \code{gam} object, produced by \code{gam()}}
  \item{view}{an array containing the names of the two predictor variables to be displayed on the 
   x and y dimensions of the plot. If omitted the first two suitable variables will be used}
  \item{cond}{a named list of the values to use for the other predictor variables (not in \code{view}). Variables omitted from 
this list will have their values set to their mean for continuous variables, or first level for factors.}
  \item{n.grid}{The number of grid nodes in each direction used for calculating the 
   plotted surface.} 
 \item{too.far}{ plot grid nodes that are too far from the points defined by the variables given in \code{view} 
can be excluded from the plot. \code{too.far} determines what is too far. The grid is scaled into the unit 
square along with the \code{view} variables and then grid nodes more than \code{too.far} from the predictor variables 
are excluded.}
\item{col}{The colours for the facets of the plot. If this is \code{NA} then if \code{se}>0 the facets are transparent, 
otherwise the colour scheme specified in \code{color} is used. If \code{col} is not \code{NA} then it is used as the facet 
colour.}
 \item{color}{ the colour scheme to use for plots when \code{se}<=0. One of \code{"topo"}, \code{"heat"}, \code{"cm"} or 
\code{"terrain"}.}
\item{se}{if less than or equal to zero then only the predicted surface is plotted, but if greater than zero, then 3 
surfaces are plotted, one at the predicted values minus \code{se} standard errors, one at the predicted values and one at
the predicted values plus \code{se} standard errors.}
\item{type}{\code{"link"} to plot on linear predictor scale and \code{"response"} to plot on the response scale.}
\item{zlim}{a two item array giving the lower and upper limits for the z-axis scale. \code{NULL} to choose automatically.}
\item{...}{other options to pass on to \code{\link{persp}}. In particular \code{ticktype="detailed"} will add proper axes 
labelling to the plots.}
}
\value{Simply produces a plot.}

\description{ Produces perspective plot views of \code{gam} model predictions, fixing all but the values in \code{view} to the 
values supplied in \code{cond}. Replaces \code{persp.gam}. 
}
\details{ The x and y limits are determined by the ranges of the variables supplied in \code{view}. If \code{se}<=0 then 
a single (height colour coded, by default) surface is produced, otherwise three (by default see-through) meshes are produced at 
mean and +/- \code{se} standard errors. Parts of the x-y plane too far from data can be excluded by setting \code{too.far}}

\author{Simon Wood \email{simon@stats.gla.ac.uk}, based on an original idea and design by Mike Lonergan.}

\seealso{
\code{\link{persp}} and \code{\link{gam}}.
}
\examples{
library(mgcv)
set.seed(0)
n<-200;sig2<-4
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
y<-x0^2+x1*x2 +runif(n,-0.3,0.3)
g<-gam(y~s(x0,x1,x2))
old.par<-par(mfrow=c(2,2))
# display the prediction surface in x0, x1 ....
vis.gam(g,ticktype="detailed",color="heat",theta=-35)  
vis.gam(g,se=2,theta=-35) # with twice standard error surfaces
vis.gam(g, view=c("x1","x2"),cond=list(x0=0.75)) # different view 
vis.gam(g, view=c("x1","x2"),cond=list(x0=0.75),theta=210,phi=40,too.far=0.07)
# ..... areas where there is no data are not plotted
par(old.par)

# Examples with factor and "by" variables

fac<-rep(1:4,20)
x<-runif(80)
y<-fac+2*x^2+rnorm(80)*0.1
fac<-factor(fac)
b<-gam(y~fac+s(x))

vis.gam(b,theta=-35,color="heat") # factor example

z<-rnorm(80)*0.4   
y<-as.numeric(fac)+3*x^2*z+rnorm(80)*0.1
b<-gam(y~fac+s(x,by=z))

vis.gam(b,theta=-35,color="heat",cond=list(z=1)) # by variable example

vis.gam(b,view=c("z","x"),theta= 35) # plot against by variable

}
\keyword{hplot} \keyword{models} \keyword{smooth} \keyword{regression}

\eof
