GAMs with GCV smoothness estimation and GAMMs by REML/PQL

GAMs with GCV smoothness estimation and GAMMs by REML/PQL [R logo]

Help pages for package `mgcv' version 1.3-22

A C D E F G I L M N P R S T U V

-- --

mgcv-package GAMs with GCV smoothness estimation and GAMMs by REML/PQL

-- A --

anova.gam Hypothesis tests related to GAM fits

-- C --

choose.k | Basis dimension choice for smooths
coef.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
coef.pdTens | Functions implementing a pdMat class for tensor product smooths
corMatrix.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix

-- D --

Dim.pdIdnot Overflow proof pdMat class for multiples of the identity matrix

-- E --

exclude.too.far | Exclude prediction grid points too far from data
extract.lme.cov | Extract the data covariance matrix from an lme object
extract.lme.cov2 | Extract the data covariance matrix from an lme object

-- F --

fix.family.link | Modify families for use in GAM fitting
fix.family.var | Modify families for use in GAM fitting
fixDependence | Detect linear dependencies of one matrix on another
formula.gam | Extract the formula from a gam object
formXtViX | Form component of GAMM covariance matrix
full.score | GCV/UBRE score for use within nlm

-- G --

gam | Generalized additive models with integrated smoothness estimation
gam.check | Some diagnostics for a fitted gam model
gam.control | Setting GAM fitting defaults
gam.convergence | GAM convergence and performance issues
gam.fit | GAM P-IRLS estimation with GCV/UBRE smoothness estimation
gam.fit2 | P-IRLS GAM estimation with GCV & UBRE derivative calculation
gam.fit3 | P-IRLS GAM estimation with GCV & UBRE derivative calculation
gam.method | Setting GAM fitting method
gam.models | Specifying generalized additive models
gam.neg.bin | GAMs with the negative binomial distribution
gam.outer | Minimize GCV or UBRE score of a GAM using `outer' iteration
gam.performance | GAM convergence and performance issues
gam.selection | Generalized Additive Model Selection
gam.setup | Generalized additive model set up
gam.side | Identifiability side conditions for a GAM
gam2derivative | Objective functions for GAM smoothing parameter estimation
gam2objective | Objective functions for GAM smoothing parameter estimation
gam3objective | Objective functions for GAM smoothing parameter estimation
gamm | Generalized Additive Mixed Models
gamm.setup | Generalized additive mixed model set up
gamObject | Fitted gam object
get.var | Get named variable or evaluate expression from list or data.frame

-- I --

influence.gam | Extract the diagonal of the influence/hat matrix for a GAM
initial.sp | Starting values for multiple smoothing parameter estimation
interpret.gam | Interpret a GAM formula

-- L --

logDet.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
logLik.gam | Extract the log likelihood for a fitted GAM

-- M --

magic | Stable Multiple Smoothing Parameter Estimation by GCV or UBRE, with optional fixed penalty
magic.post.proc | Auxilliary information from magic fit
mgcv | Multiple Smoothing Parameter Estimation by GCV or UBRE
mgcv.control | Setting mgcv defaults
mgcv.find.theta | GAMs with the negative binomial distribution
mgcv.get.scale | GAMs with the negative binomial distribution
mono.con | Monotonicity constraints for a cubic regression spline
mroot | Smallest square root of matrix

-- N --

new.name | Obtain a name for a new variable that is not already in use
notExp | Functions for better-than-log positive parameterization
notExp2 | Alternative to log parameterization for variance components
notLog | Functions for better-than-log positive parameterization
notLog2 | Alternative to log parameterization for variance components
null.space.dimension | The basis of the space of un-penalized functions for a TPRS

-- P --

p.spline | Constructor functions for smooth terms in a GAM
pcls | Penalized Constrained Least Squares Fitting
pdConstruct.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
pdConstruct.pdTens | Functions implementing a pdMat class for tensor product smooths
pdFactor.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
pdFactor.pdTens | Functions implementing a pdMat class for tensor product smooths
pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
pdMatrix.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
pdMatrix.pdTens | Functions implementing a pdMat class for tensor product smooths
pdTens | Functions implementing a pdMat class for tensor product smooths
persp.gam | Visualization of GAM objects
place.knots | Automatically place a set of knots evenly through covariate values
plot.gam | Default GAM plotting
predict.gam | Prediction from fitted GAM model
Predict.matrix | Prediction methods for smooth terms in a GAM
Predict.matrix.cr.smooth | Prediction methods for smooth terms in a GAM
Predict.matrix.cs.smooth | Prediction methods for smooth terms in a GAM
Predict.matrix.cyclic.smooth | Prediction methods for smooth terms in a GAM
Predict.matrix.tensor.smooth | Prediction methods for smooth terms in a GAM
Predict.matrix.tprs.smooth | Prediction methods for smooth terms in a GAM
Predict.matrix.ts.smooth | Prediction methods for smooth terms in a GAM
PredictMat | Prediction/Construction wrapper functions for GAM smooth terms
print.anova.gam | Hypothesis tests related to GAM fits
print.gam | Generalized Additive Model default print statement
print.summary.gam | Summary for a GAM fit

-- R --

residuals.gam Generalized Additive Model residuals

-- S --

s | Defining smooths in GAM formulae
smooth.construct | Constructor functions for smooth terms in a GAM
smooth.construct.cc.smooth.spec | Constructor functions for smooth terms in a GAM
smooth.construct.cr.smooth.spec | Constructor functions for smooth terms in a GAM
smooth.construct.cs.smooth.spec | Constructor functions for smooth terms in a GAM
smooth.construct.tensor.smooth.spec | Constructor functions for smooth terms in a GAM
smooth.construct.tp.smooth.spec | Constructor functions for smooth terms in a GAM
smooth.construct.ts.smooth.spec | Constructor functions for smooth terms in a GAM
smoothCon | Prediction/Construction wrapper functions for GAM smooth terms
solve.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
step.gam | Alternatives to step.gam
summary.gam | Summary for a GAM fit
summary.pdIdnot | Overflow proof pdMat class for multiples of the identity matrix
summary.pdTens | Functions implementing a pdMat class for tensor product smooths

-- T --

te | Define tensor product smooths in GAM formulae
tensor.prod.model.matrix | Utility functions for constructing tensor product smooths
tensor.prod.penalties | Utility functions for constructing tensor product smooths

-- U --

uniquecombs find the unique rows in a matrix

-- V --

vcov.gam | Extract parameter (estimator) covariance matrix from GAM fit
vis.gam | Visualization of GAM objects


Hypothesis tests related to GAM fits

Description

Performs hypothesis tests relating to one or more fitted gam objects. For a single fitted gam object, Wald tests of the significance of each parametric and smooth term are performed. Otherwise the fitted models are compared using an analysis of deviance table. The tests are usually approximate, unless the models are un-penalized.

Usage

## S3 method for class 'gam':
anova(object, ..., dispersion = NULL, test = NULL)
## S3 method for class 'anova.gam':
print(x, digits = max(3, getOption("digits") - 3),...)

Arguments

object,... | fitted model objects of class gam as produced by gam().
x | an anova.gam object produced by a single model call to anova.gam().
dispersion | a value for the dispersion parameter: not normally used.
test | what sort of test to perform for a multi-model call. One of "Chisq", "F" or "Cp".
digits | number of digits to use when printing output.

Details

If more than one fitted model is provided than anova.glm is used. If only one model is provided then the significance of each model term is assessed using Wald tests: see summary.gam for details of the actual computations. In the latter case print.anova.gam is used as the printing method.

P-values are usually reliable if the smoothing parameters are known, or the model is unpenalized. If smoothing parameters have been estimated then the p-values are typically somewhat too low under the null. This occurs because the uncertainty associated with the smoothing parameters is neglected in the calculations of the distributions under the null, which tends to lead to underdispersion in these distributions, and in turn to p-value estimates that are too low. (In simulations where the null is correct, I have seen p-values that are as low as half of what they should be.) Note however that tests can have low power if the estimated rank of the test statistic is much higher than the EDF, so that p-values can also be too low in some cases.

If it is important to have p-values that are as accurate as possible, then, at least in the single model case, it is probably advisable to perform tests using unpenalized smooths (i.e. s(...,fx=TRUE)) with the basis dimension, k, left at what would have been used with penalization. Such tests are not as powerful, of course, but the p-values are more accurate under the null. Whether or not extra accuracy is required will usually depend on whether or not hypothesis testing is a key objective of the analysis.

Value

In the multi-model case anova.gam produces output identical to anova.glm, which it in fact uses.
In the single model case an object of class anova.gam is produced, which is in fact an object returned from summary.gam.
print.anova.gam simply produces tabulated output.

WARNING

P-values are only approximate, particularly as a result of ignoring smoothing parameter uncertainty.

Author(s)

Simon N. Wood simon.wood@r-project.org with substantial improvements by Henric Nilsson.

See Also

gam, predict.gam, gam.check, summary.gam

Examples

library(mgcv)
set.seed(0)
n<-200
sig<-2
x0 <- rep(1:4,50)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
y <- 2 * x0
y <- y + exp(2 * x1)
y <- y + 0.2 * x2^11 * (10 * (1 - x2))^6 + 10 * (10 * x2)^3 * (1 - x2)^10
e <- rnorm(n, 0, sig)
y <- y + e
x0<-as.factor(x0)
b<-gam(y~x0+s(x1)+s(x2)+s(x3))
anova(b)
b1<-gam(y~x0+s(x1)+s(x2))
anova(b,b1,test="F")
[Package mgcv version 1.3-22 Index]


Basis dimension choice for smooths

Description

Choosing the basis dimension, and checking the choice, when using penalized regression smoothers.

Penalized regression smoothers gain computational efficiency by virtue of being defined using a basis of relatively modest size, k. When setting up models in the mgcv package, using s or te terms in a model formula, k must be chosen.

In practice k-1 sets the upper limit on the degrees of freedom associated with an s smooth (1 degree of freedom is lost to the identifiability constraint on the smooth). For te smooths the upper limit of the degrees of freedom is given by the product of the k values provided for each marginal smooth less one, for the constraint. However the actual effective degrees of freedom are controlled by the degree of penalization selected during fitting, by GCV, AIC or whatever is specified. The exception to this is if a smooth is specified using the fx=TRUE option, in which case it is unpenalized.

So, exact choice of k is not generally critical: it should be chosen to be large enough that you are reasonably sure of having enough degrees of freedom to represent the underlying `truth' reasonably well, but small enough to maintain reasonable computational efficiency. Clearly `large' and `small' are dependent on the particular problem being addressed.

As with all model assumptions, it is useful to be able to check the choice of k informally. If the effective degrees of freedom for a model term are estimated to be much less than k-1 then this is unlikely to be very worthwhile, but as the EDF approach k-1, checking can be important. A useful general purpose approach goes as follows: (i) fit your model and extract the deviance residuals; (ii) for each smooth term in your model, fit an equivalent, single, smooth to the residuals, using a substantially increased k to see if there is pattern in the residuals that could potentially be explained by increasing k. Examples are provided below.

More sophisticated approaches based on partial residuals are also possible.

One scenario that can cause confusion is this: a model is fitted with k=10 for a smooth term, and the EDF for the term is estimated as 7.6, some way below the maximum of 9. The model is then refitted with k=20 and the EDF increases to 8.7 - what is happening - how come the EDF was not 8.7 the first time around? The explanation is that the function space with k=20 contains a larger subspace of functions with EDF 8.7 than did the function space with k=10: one of the functions in this larger subspace fits the data a little better than did any function in the smaller subspace. These subtleties seldom have much impact on the statistical conclusions to be drawn from a model fit, however.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Wood, S.N. (2006) Generalized Additive Models: An Introduction with R. CRC.

http://www.maths.bath.ac.uk/~sw283/

Examples

## Simulate some data ....
library(mgcv)
set.seed(0)
n<-400;sig<-2
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sig)
y <- f + e
## fit a GAM with quite low `k'
b<-gam(y~s(x0,k=6)+s(x1,k=6)+s(x2,k=6)+s(x3,k=6))
plot(b,pages=1)

## check for residual pattern, removeable by increasing `k'
## typically `k', below, chould be substantially larger than
## the original, `k' but certainly less than n/2.
## Note use of cheap "cs" shrinkage smoothers, and gamma=1.4
## to reduce chance of overfitting...
rsd <- residuals(b)
gam(rsd~s(x0,k=40,bs="cs"),gamma=1.4) ## fine
gam(rsd~s(x1,k=40,bs="cs"),gamma=1.4) ## fine
gam(rsd~s(x2,k=40,bs="cs"),gamma=1.4) ## original `k' too low
gam(rsd~s(x3,k=40,bs="cs"),gamma=1.4) ## fine

## similar example with multi-dimensional smooth
b1 <- gam(y~s(x0)+s(x1,x2,k=15)+s(x3))
rsd <- residuals(b1)
gam(rsd~s(x0,k=40,bs="cs"),gamma=1.4) ## fine
gam(rsd~s(x1,x2,k=100,bs="ts"),gamma=1.4) ## original `k' too low
gam(rsd~s(x3,k=40,bs="cs"),gamma=1.4) ## fine

## and a `te' example
b2 <- gam(y~s(x0)+te(x1,x2,k=4)+s(x3))
rsd <- residuals(b2)
gam(rsd~s(x0,k=40,bs="cs"),gamma=1.4) ## fine
gam(rsd~te(x1,x2,k=10,bs="cs"),gamma=1.4) ## original `k' too low
gam(rsd~s(x3,k=40,bs="cs"),gamma=1.4) ## fine

## same approach works with other families in the original model
g<-exp(f/4)
y<-rpois(rep(1,n),g)
bp<-gam(y~s(x0,k=6)+s(x1,k=6)+s(x2,k=6)+s(x3,k=6),family=poisson)
rsd <- residuals(bp)
gam(rsd~s(x0,k=40,bs="cs"),gamma=1.4) ## fine
gam(rsd~s(x1,k=40,bs="cs"),gamma=1.4) ## fine
gam(rsd~s(x2,k=40,bs="cs"),gamma=1.4) ## original `k' too low
gam(rsd~s(x3,k=40,bs="cs"),gamma=1.4) ## fine


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 TRUE if the corresponding node is too far from data and FALSE otherwise. Basically a service routine for vis.gam and plot.gam.

Usage

exclude.too.far(g1,g2,d1,d2,dist)

Arguments

g1 | co-ordinates of grid relative to first axis.
g2 | co-ordinates of grid relative to second axis.
d1 | co-ordinates of data relative to first axis.
d2 | co-ordinates of data relative to second axis.
dist | how far away counts as too far. Grid and data are first scaled so that the grid lies exactly in the unit square, and 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 g1 and g2 lies exactly in the unit square (i.e. on [0,1] by [0,1]). These scalings are applied to g1, g2, d1 and d2. The minimum Euclidean distance from each node to a datum is then determined and if it is greater than dist the corresponding entry in the returned array is set to TRUE (otherwise to FALSE). The distance calculations are performed in compiled code for speed without storage overheads.

Value

A logical array with TRUE indicating a node in the grid defined by g1, g2 that is `too far' from any datum.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

vis.gam

Examples

library(mgcv)
x<-rnorm(100);y<-rnorm(100) # some "data"
n<-40 # 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)


Extract the data covariance matrix from an lme object

Description

This is a service routine for gamm. Extracts the estimated covariance matrix of the data from an lme object, allowing the user control about which levels of random effects to include in this calculation. extract.lme.cov forms the full matrix explicitly: extract.lme.cov2 tries to be more economical than this.

Usage

extract.lme.cov(b,data,start.level=1)
extract.lme.cov2(b,data,start.level=1)

Arguments

b | A fitted model object returned by a call to lme
data | The data frame/ model frame that was supplied to lme.
start.level | The level of nesting at which to start including random effects in the calculation. This is used to allow smooth terms to be estimated as random effects, but treated like fixed effects for variance calculations.

Details

The random effects, correlation structure and variance structure used for a linear mixed model combine to imply a covariance matrix for the response data being modelled. These routines extracts that covariance matrix. The process is slightly complicated, because different components of the fitted model object are stored in different orders (see function code for details!).

The extract.lme.cov calculation is not optimally efficient, since it forms the full matrix, which may in fact be sparse. extract.lme.cov2 is more efficient. If the covariance matrix is diagonal, then only the leading diagonal is returned; if it can be written as a block diagonal matrix (under some permutation of the original data) then a list of matrices defining the non-zero blocks is returned along with an index indicating which row of the original data each row/column of the block diagonal matrix relates to. The block sizes are defined by the coarsest level of grouping in the random effect structure.

gamm uses extract.lme.cov2.

extract.lme.cov does not currently deal with the situation in which the grouping factors for a correlation structure are finer than those for the random effects. extract.lme.cov2 does deal with this situation.

Value

For extract.lme.cov an estimated covariance matrix.
For extract.lme.cov2 a list containing the estimated covariance matrix and an indexing array. The covariance matrix is stored as the elements on the leading diagonal, a list of the matrices defining a block diagonal matrix, or a full matrix if the previous two options are not possible.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

For lme see:

Pinheiro J.C. and Bates, D.M. (2000) Mixed effects Models in S and S-PLUS. Springer

For details of how GAMMs are set up here for estimation using lme see:

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for Generalized Additive Mixed Models. Biometrics

http://www.maths.bath.ac.uk/~sw283/

See Also

gamm, formXtViX

Examples

library(nlme)
data(Rail)
b <- lme(travel~1,Rail,~1|Rail)
extract.lme.cov(b,Rail)


Modify families for use in GAM fitting

Description

Generalized Additive Model fitting by `outer' iteration, requires extra derivatives of the variance and link functions to be added to family objects. The functions add what is needed.

Usage

fix.family.link(fam)
fix.family.var(fam)

Arguments

fam A family.

Details

Outer iteration GAM estimation requires derivatives of the GCV/UBRE score, which are obtained by differentiating the P-IRLS GAM fitting iteration w.r.t. the model smoothing parameters. The expressions for the derivatives require the second derivative of the link w.r.t. the mean, and the first derivative of the variance function w.r.t. the mean. These functions add functions evaluating these quantities to a family.

If the family already has functions dvar and d2link then these functions simply return the family unmodified: this allows non-standard links to be used with gam when using outer iteration (performance iteration operates with unmodified families).

The dvar function is a function of a mean vector, mu, and returns a vector of corresponding first derivatives of the family variance function. The d2link function is also a function of a vector of mean values, mu: it returns a vector of second derivatives of the link, evaluated at mu.

If modifying your own family, note that you can often get away with supplying only a dvar function if your family only requires links that occur in one of the standard families.

Value

A family object with extra component functions dvar and d2link.

Author(s)

Simon N. Wood simon.wood@r-project.org


Detect linear dependencies of one matrix on another

Description

Identifies columns of a matrix X2 which are linearly dependent on columns of a matrix X1. Primarily of use in setting up identifiability constraints for nested GAMs.

Usage

fixDependence(X1,X2,tol=.Machine$double.eps^.5)

Arguments

X1 | A matrix.
X2 | A matrix, the columns of which may be partially linearly dependent on the columns of X1.
tol | The tolernce to use when assessing linear dependence.

Details

The algorithm uses a simple approach based on QR decomposition: see code for details.

Value

An array of the columns of X2 which are linearly dependent on columns of X1. NULL if the two matrices are independent.

Author(s)

Simon N. Wood simon.wood@r-project.org

Examples

n<-20;c1<-4;c2<-7
X1<-matrix(runif(n*c1),n,c1)
X2<-matrix(runif(n*c2),n,c2)
X2[,3]<-X1[,2]+X2[,4]*.1
X2[,5]<-X1[,1]*.2+X1[,2]*.04
fixDependence(X1,X2)


Extract the formula from a gam object

Description

Extracts the formula from a fitted gam object.

Usage

formula.gam(x,...)

Arguments

x | fitted model objects of class gam (see gamObject) as produced by gam().
... | un-used in this case

Details

Returns x$formula. Provided so that anova methods print an appropriate description of the model.

Value

A model formula.

Author(s)

Simon N. Wood simon.wood@r-project.org

See Also

gam


Form component of GAMM covariance matrix

Description

This is a service routine for gamm. Given, V, an estimated covariance matrix obtained using extract.lme.cov2 this routine forms X'inv(V)X as efficiently as possible, given the structure of V (usually sparse).

Usage

formXtViX(V,X)

Arguments

V | A data covariance matrix list returned from extract.lme.cov2
X | A model matrix.

Details

The covariance matrix returned by extract.lme.cov2 may be in a packed and re-ordered format, since it is usually sparse. Hence a special service routine is required to form the required products involving this matrix.

Value

A matrix.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

For lme see:

Pinheiro J.C. and Bates, D.M. (2000) Mixed effects Models in S and S-PLUS. Springer

For details of how GAMMs are set up here for estimation using lme see:

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for Generalized Additive Mixed Models. Biometrics

http://www.maths.bath.ac.uk/~sw283/

See Also

gamm, extract.lme.cov2


GCV/UBRE score for use within nlm

Description

Evaluates GCV/UBRE score for a GAM, given smoothing parameters. The routine calls gam.fit to fit the model, and is usually called by nlm to optimize the smoothing parameters.

This is basically a service routine for gam, and is not usually called directly by users. It is only used in this context for GAMs fitted by outer iteration (see gam.outer) when the the outer method is "nlm.fd" (see gam.method).

Usage

full.score(sp,G,family,control,gamma,pearson,...)

Arguments

sp | The logs of the smoothing parameters
G | a list returned by gam.setup
family | The family object for the GAM.
control | a list returned be gam.control
gamma | the degrees of freedom inflation factor (usually 1).
pearson | whether the GCV/UBRE score should be based on the Pearson statistic or the deviance. The deviance is usually preferable.
... | other arguments, typically for passing on to gam.fit.

Value

The value of the GCV/UBRE score, with attribute "full.gam.object" which is the full object returned by gam.fit.

Author(s)

Simon N. Wood simon.wood@r-project.org


Generalized additive models with integrated smoothness estimation

Description

Fits a generalized additive model (GAM) to data. The degree of smoothness of model terms is estimated as part of fitting; isotropic or scale invariant smooths of any number of variables are available as model terms; confidence/credible intervals are readily available for any quantity predicted using a fitted model; gam is extendable: i.e. users can add smooths.

Smooth terms are represented using penalized regression splines (or similar smoothers) 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 (isotropic) or tensor product splines (when an isotropic smooth is inappropriate). For more on specifying models see gam.models. For more on model selection see gam.selection. For faster fits use the "cr" bases for smooth terms, te smooths for smooths of several variables, and performance iteration for smoothing parameter estimation (see gam.method). For large datasets see warnings.

gam() is not a clone of what S-PLUS provides: the major differences are (i) that by default estimation of the degree of smoothness of model terms is part of model fitting, (ii) a Bayesian approach to variance estimation is employed that makes for easier confidence interval calculation (with good coverage probabilites) and (iii) the facilities for incorporating smooths of more than one variable are different: specifically there are no lo smooths, but instead (a) s terms can have more than one argument, implying an isotropic smooth and (b) te smooths are provided as an effective means for modelling smooth interactions of any number of variables via scale invariant tensor product smooths. If you want a clone of what S-PLUS provides use gam from package gam.

Usage

gam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
    na.action,offset=NULL,control=gam.control(),method=gam.method(),
    scale=0,knots=NULL,sp=NULL,min.sp=NULL,H=NULL,gamma=1,
    fit=TRUE,G=NULL,in.out,...)

Arguments

formula A GAM formula (see also 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 y ~ . is not allowed). Smooth terms are specified by expressions of the form:
s(var1,var2,...,k=12,fx=FALSE,bs="tp",by=a.var) where var1, var2, etc. are the covariates which the smooth is a function of and k is the dimension of the basis used to represent the smooth term. If k is not specified then k=10*3^(d-1) is used where d is the number of covariates for this term. fx is used to indicate whether or not this term has a fixed number of degrees of freedom (fx=FALSE to select d.f. by GCV/UBRE). bs indicates the basis to use for the smooth: for a full list see s, but note that the default "tp", while it possesses nice optimality properties is slow and memory hungry for very large datasets (but see examples for how to get around this). by can be used to specify a variable by which the smooth should be multiplied. For example gam(y~z+s(x,by=z)) would specify a model E(y)=f(x)z where f(.) is a smooth function (the formula is y~x+s(x,by=z) rather than y~s(x,by=z) because the smooths are always set up to sum to zero over the covariate values). The by option is particularly useful for models in which different functions of the same variable are required for each level of a factor and for `variable parameter models': see s.
An alternative for specifying smooths of more than one covariate is e.g.:
te(x,z,bs=c("tp","tp"),m=c(2,3),k=c(5,10)) which would specify a tensor product smooth of the two covariates x and z constructed from marginal t.p.r.s. bases of dimension 5 and 10 with marginal penalties of order 2 and 3. Any combination of basis types is possible, as is any number of covariates.
Formulae can involve nested or ``overlapping'' terms such as
y~s(x)+s(z)+s(x,z) or y~s(x,z)+s(z,v): see gam.side for further details and examples.
family | This is a family object specifying the distribution and link to use in fitting etc. See glm and family for more details. The negative binomial families provided by the MASS library can be used, with or without known theta parameter: see gam.neg.bin for details.
data | A data frame containing the model response variable and covariates required by the formula. By default the variables are taken from environment(formula): typically the environment from which gam is called.
weights | prior weights on the data.
subset | an optional vector specifying a subset of observations to be used in the fitting process.
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'.
offset | Can be used to supply a model offset for use in fitting. Note that this offset will always be completely ignored when predicting, unlike an offset included in formula: this conforms to the behaviour of lm and glm.
control | A list of fit control parameters returned by gam.control.
method | A list controlling the fitting methods used. This can make a big difference to computational speed, and, in some cases, reliability of convergence: see gam.method for details.
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 theta must be known. If 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 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.
knots | this is an optional list containing user specified knot values to be used for basis construction. For the cr and cc bases the user simply supplies the knots to be used, and there must be the same number as the basis dimension, k, for the smooth concerned. For the tp basis knots has two uses. Firstly, for large datasets the calculation of the tp basis can be time-consuming. The user can retain most of the advantages of the t.p.r.s. approach by supplying a reduced set of covariate values from which to obtain the basis - typically the number of covariate values used will be substantially smaller than the number of data, and substantially larger than the basis dimension, k. 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, 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.
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 "magic" (see gam.control and 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 "mgcv", if sp is supplied then all its elements must be positive. Note that fx=TRUE in a smooth term over-rides what is supplied here effectively setting the smoothing parameter to zero.
min.sp | for fit method "magic" only, lower bounds can be supplied for the smoothing parameters. Note that if this option is used then the smoothing parameters 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.
H | With fit method "magic" 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.
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 "magic".
fit | If this argument is TRUE then gam sets up the model and fits it, but if it is FALSE then the model is set up and an object G is returned which is the output from gam.setup plus some extra items required to complete the GAM fitting process.
G | Usually NULL, but may contain the object returned by a previous call to gam with fit=FALSE, in which case all other arguments are ignored except for gamma, in.out, control, method and fit.
in.out | optional list for initializing outer iteration. If supplied then this must contain two elements: sp should be an array of initialization values for all smoothing parameters (there must be a value for all smoothing parameters, whether fixed or to be estimated, but those for fixed s.p.s are not used); scale is the typical scale of the GCV/UBRE function, for passing to the outer optimizer.
... further arguments for passing on e.g. to gam.fit

(such as mustart).

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:

log(E(y_i))=f_1(x_1i)+f_2(x_2i)

where the (independent) response variables y_i~Poi, and f_1 and f_2 are smooth functions of covariates x_1 and 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 f_1 and 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 mgcv implementation of 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.

gam in mgcv solves the smoothing parameter estimation problem by using the Generalized Cross Validation (GCV) criterion

n D/(n - DoF)^2

or an Un-Biased Risk Estimator (UBRE )criterion

D/n + 2 s DoF / n -s

where D is the deviance, n the number of data, s the scale parameter and DoF the effective degrees of freedom of the model. Notice that UBRE is effectively just AIC rescaled, but is only used when s is known. It is also possible to replace D by the Pearson statistic (see gam.method), but this can lead to over smoothing. Smoothing parameters are chosen to minimize the GCV or UBRE score for the model, and the main computational challenge solved by the mgcv package is to do this efficiently and reliably. Various alternative numerical methods are provided: see gam.method.

Broadly gam works by first constructing basis functions and one or more quadratic penalty coefficient matrices 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 gam.fit, a modification of 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. Alternatively the P-IRLS scheme is iterated to convergence for each trial set of smoothing parameters, and GCV or UBRE scores are only evaluated on convergence - optimization is then `outer' to the P-IRLS loop: in this case extra derivatives have to be carried along with the P-IRLS iteration, to facilitate optimization, and gam.fit2 is used.

Five alternative basis-penalty types are built in for representing model smooths, but alternatives can easily be added (see smooth.construct which uses p-splines to illustrate how to add new smooths). The built in alternatives for univariate smooths terms are: a conventional penalized cubic regression spline basis, parameterized in terms of the function values at the knots; a cyclic cubic spline with a similar parameterization and thin plate regression splines. The cubic spline bases are computationally very efficient, but require `knot' locations to be chosen (automatically by default). The thin plate regression splines are optimal low rank smooths which do not have knots, but are more computationally costly to set up. Smooths of several variables can be represented using thin plate regression splines, or tensor products of any available basis including user defined bases (tensor product penalties are obtained automatically form the marginal basis penalties). The t.p.r.s. basis is isotropic, so if this is not appropriate tensor product terms should be used. Tensor product smooths have one penalty and smoothing parameter per marginal basis, which means that the relative scaling of covariates is essentially determined automatically by GCV/UBRE. The t.p.r.s. basis and cubic regression spline bases are both available with either conventional `wiggliness penalties' or penalties augmented with a shrinkage component: the conventional penalties treat some space of functions as `completely smooth' and do not penalize such functions at all; the penalties with extra shrinkage will zero a term altogether for high enough smoothing parameters: gam.selection has an example of the use of such terms.

For any 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), or use tensor products of cheaper bases.

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 x then there would be a knot at every 10th (ordered) 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 constrained 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 the underlying fitting methods are given in Wood (2000, 2004a, 2006b).

Value

If fit == FALSE the function returns a list G of items needed to fit a GAM, but doesn't actually fit it.
Otherwise the function returns an object of class "gam" as described in gamObject.

WARNINGS

If fit method "mgcv" is selected, the code does not check for rank deficiency 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.

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

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 gam.control for details.

Author(s)

Simon N. Wood simon.wood@r-project.org

Front end design inspired by the S function of the same name based on the work of Hastie and Tibshirani (1990). Underlying methods owe much to the work of Wahba (e.g. 1990) and Gu (e.g. 2002).

References

Key References on this implementation:

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

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

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

Wood, S.N. (2004b) On confidence intervals for GAMs based on penalized regression splines. Technical Report 04-12 Department of Statistics, University of Glasgow.

Wood, S.N. (2006a) Low rank scale invariant tensor product smooths for generalized additive mixed models. Biometrics

Wood S.N. (2006b) Generalized Additive Models: An Introduction with R. CRC Press.

Key Reference on GAMs and related models:

Hastie (1993) in Chambers and Hastie (1993) Statistical Models in S. Chapman and Hall.

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

Wahba (1990) Spline Models of Observational Data. SIAM

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.

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

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

http://www.maths.bath.ac.uk/~sw283/

See Also

mgcv-package, gamObject, gam.models, s, predict.gam,plot.gam, summary.gam, gam.side, gam.selection,mgcv, gam.control gam.check, gam.neg.bin, magic,vis.gam

Examples

library(mgcv)
set.seed(0)
n<-400
sig<-2
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
f0 <- function(x) 2 * sin(pi * x)
f1 <- function(x) exp(2 * x)
f2 <- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
f3 <- function(x) 0*x
f <- f0(x0) + f1(x1) + f2(x2)
e <- rnorm(n, 0, sig)
y <- f + e
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))
summary(b)
plot(b,pages=1,residuals=TRUE)
# same fit in two parts .....
G<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),fit=FALSE)
b<-gam(G=G)
# an extra ridge penalty (useful with convergence problems) ....
bp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),H=diag(0.5,37))
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,k=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/4)
y<-rpois(rep(1,n),g)
b2<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson)
plot(b2,pages=1)
# repeat fit using performance iteration
gm <- gam.method(gam="perf.magic")
b3<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,method=gm)
plot(b3,pages=1)

# a binary example
g <- (f-5)/3
g <- binomial()$linkinv(g)
y <- rbinom(g,1,g)
lr.fit <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=binomial)
## plot model components with truth overlaid in red
op <- par(mfrow=c(2,2))
for (k in 1:4) {
  plot(lr.fit,residuals=TRUE,select=k)
  xx <- sort(eval(parse(text=paste("x",k-1,sep=""))))
  ff <- eval(parse(text=paste("f",k-1,"(xx)",sep="")))
  lines(xx,(ff-mean(ff))/3,col=2)
}
par(op)
anova(lr.fit)
lr.fit1 <- gam(y~s(x0)+s(x1)+s(x2),family=binomial)
lr.fit2 <- gam(y~s(x1)+s(x2),family=binomial)
AIC(lr.fit,lr.fit1,lr.fit2)

# 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")


Objective functions for GAM smoothing parameter estimation

Description

Estimation of GAM smoothing parameters is most stable if optimization of the UBRE or GCV score is outer to the penalized iteratively re-weighted least squares scheme used to estimate the model given smoothing parameters. These functions evaluate the GCV/UBRE score of a GAM model, given smoothing parameters, in a manner sutable for use by optim or nlm. Not normally called directly, but rather service routines for gam.outer.

Usage

gam2objective(lsp,args,printWarn=FALSE,...)
gam2derivative(lsp,args,...)
gam3objective(lsp,args,...)

Arguments

lsp | The log smoothing parameters.
args | List of arguments required to call gam.fit2.
printWarn | Should gam.fit2 print some warnings? Used to suppress warnings that are only of interest for the final fitted model, until that model is reached.
... | Other arguments for passing to gam.fit2.

Details

gam2objective and gam2derivative are functions suitable for calling by optim, to evaluate the GCV/UBRE score and it's derivatives w.r.t. log smoothing parameters.

gam3objective is an equaivalent to gam2objective, suitable for optimization by nlm - derivatives of the GCV/UBRE function are calculated and returned as attributes.

The basic idea of optimizing smoothing parameters `outer' to the P-IRLS loop was first proposed in O'Sullivan et al. (1986).

Author(s)

Simon N. Wood simon.wood@r-project.org

References

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam.fit2, gam, mgcv, magic


Some diagnostics for a fitted gam model

Description

Takes a fitted gam object produced by gam() and produces some diagnostic information about the fitting procedure and results.

Usage

gam.check(b)

Arguments

b a fitted gam object as produced by gam().

Details

This function plots 4 standard diagnostic plots, and some other convergence diagnostics. Output differs depending on whether the underlying fitting method was mgcv or another method (see gam.method).

For mgcv based fits , 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).

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 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, gam.fit switches to a cautious optimization mode after a user controlled number of iterations of the IRLS algorithm (see 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 trace to TRUE in 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 based on magic or gam.fit2 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 standardized 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 mgcv, a message is printed if the second guess smoothing parameters did not improve on the first guess - this is primarily there for the developer. For other fit methods the estimated rank of the model is printed.

The ideal results from this function have a smooth, single minima GCV/UBRE plot (where appropriate), 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 mgcv.conv part of the gam object.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

choose.k, gam, mgcv, magic

Examples

library(mgcv)
set.seed(0)
n<-200
sig<-2
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
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, sig)
y <- y + e
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))
plot(b,pages=1)
gam.check(b)


Setting GAM fitting defaults

Description

This is an internal function of package 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. To change the default choise of fitting method, see gam.method.

Usage

gam.control(irls.reg=0.0,epsilon = 1e-06, maxit = 100,globit = 20,
            mgcv.tol=1e-7,mgcv.half=15,nb.theta.mult=10000, trace = FALSE,
            rank.tol=.Machine$double.eps^0.5,absorb.cons=TRUE,
            max.tprs.knots=5000,nlm=list(),optim=list(),newton=list(),
            outerPIsteps=1)

Arguments

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 irls.reg is the size of the penalty. The penalty can not be used if the underlying fitting method is mgcv (not the default - see gam.method for details).
epsilon | This is used for judging conversion of the GLM IRLS loop in gam.fit or gam.fit2 .
maxit | Maximum number of IRLS iterations to perform using cautious GCV/UBRE optimization, after globit IRLS iterations with normal GCV optimization have been performed. Note that only fitting based on mgcv (not default) makes any distinction between cautious and global optimization.
globit | Maximum number of IRLS iterations to perform with normal GCV/UBRE optimization. If convergence is not achieved after these iterations then a further maxit iterations will be performed using cautious GCV/UBRE optimization.
mgcv.tol | The convergence tolerance parameter to use in GCV/UBRE optimization.
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.
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 nb.theta.mult, while the minimum value is set to the initial value divided by nb.theta.mult.
trace | Set this to TRUE to turn on diagnostic output.
rank.tol | The tolerance used to estimate the rank of the fitting problem, for methods which deal with rank deficient cases (basically all except those based on mgcv).
absorb.cons | If TRUE then the GAM is set up using a parameterization which requires no further constraint. Usually this means that all the smooths are automatically centered (i.e. they sum to zero over the covariate values). If FALSE then the ordinary parameterizations of the smooths are used, which require constraints to be imposed during fitting.
max.tprs.knots | This is the default initial maximum number of knots to allow when constructing a t.p.r.s bases (bs="tp"). The set up cost (and storage) for these smooths scales as the square of the number of initial knots, so if it's too high you can appear to freeze R. Usually one would want to use an alternative smoothing basis (or te terms), or the approach illustrated in the examples in gam, rather than simply increasing this default.
nlm | list of control parameters to pass to nlm if this is used for outer estimation of smoothing parameters. See details.
optim | list of control parameters to pass to optim if this is used for outer estimation of smoothing parameters. See details.
newton | list of control parameters to pass to default Newton optimizer used for outer estimation of log smoothing parameters. See details.
outerPIsteps | The number of performance interation steps used to initialize outer iteration. Less than 1 means that only one performance iteration step is taken to get the function scale, but the corresponding smoothing parameter estimates are discarded.

Details

When outer iteration is used for fitting then the control list nlm stores control arguments for calls to routine nlm. The list has the following named elements: (i) ndigit is the number of significant digits in the GCV/UBRE score - by default this is worked out from epsilon; (ii) gradtol is the tolerance used to judge convergence of the gradient of the GCV/UBRE score to zero - by default set to 100*epsilon; (iii) stepmax is the maximum allowable log smoothing parameter step - defaults to 2; (iv) steptol is the minimum allowable step length - defaults to 1e-4; (v) iterlim is the maximum number of optimization steps allowed - defaults to 200; (vi) check.analyticals indicates whether the built in exact derivative calculations should be checked numerically - defaults to FALSE. Any of these which are not supplied and named in the list are set to their default values.

Outer iteration using newton is controlled by the list newton with the following elements : conv.tol (default 1e-6) is the relative convergence tolerance; maxNstep is the maximum length allowed for an element of the Newton search direction (default 5); maxSstep is the maximum length allowed for an element of the steepest descent direction (only used if Newton fails - default 2); maxHalf is the maximum number of step halvings to permit before giving up (default 30); use.svd determines whether to use an SVD step or QR step for the part of fitting that determines problem rank, etc. (defaults to FALSE).

Outer iteration using optim is controlled using list optim, which currently has one element: factr which takes default value 1e7.

When fitting is been done by calls to routine mgcv, maxit and globit control the maximum iterations of the IRLS algorithm, as follows: the algorithm will first execute up to 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 globit iterations, then a further 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.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

Wood, S.N. (2004) Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Amer. Statist. Ass.

http://www.maths.bath.ac.uk/~sw283/

See Also

gam.method gam, gam.fit, glm.control


GAM convergence and performance issues

Description

When fitting GAMs there is a tradeoff between speed of fitting and probability of fit convergence. The default fitting options specified by gam.method (as the default for argument method of gam), always opt for certainty of convergence over speed of fit. In the additive modelling contexts this means using fitting routine magic rather than the slightly faster routine mgcv. In the Generalized Additive Model case it means using `outer' iteration in preference to `performance iteration': see gam.outer for details.

It is possible for the default `outer' iteration to fail when finding intial smoothing parameters using a few steps of outer iteration (if you get a convergence failure message from magic when outer iterating, then this is what has happened): lower outerPIsteps in gam.control to fix this.

There are two things that you can do to speed up GAM fitting. (i) Change the method argument to gam so that `performance iteration' is used in place of the default outer iteration. See the perf.magic option under gam.method, for example. Usually performance iteration converges well and is quick. (ii) For large datasets it may be worth changing the smoothing basis to use bs="cr" (see s for details) for 1-d smooths, and to use te smooths in place of s smooths for smooths of more than one variable. This is because the default thin plate regression spline basis "tp" is costly to set up for large datasets (much over 1000 data, say). Alternatively see the last few examples for gam.

If the GAM estimation process fails to converge when using performance iteration, then switch to outer iteration via the method argument of gam (see gam.method). If it still fails, try increasing the number of IRLS iterations (see gam.control) or perhaps experiment with the convergence tolerance.

If you still have problems, it's worth noting that a GAM is just a (penalized) GLM and the IRLS scheme used to estimate GLMs is not guaranteed to converge. Hence non convergence of a GAM may relate to a lack of stability in the basic IRLS scheme. Therefore it is worth trying to establish whether the IRLS iterations are capable of converging. To do this fit the problematic GAM with all smooth terms specified with fx=TRUE so that the smoothing parameters are all fixed at zero. If this `largest' model can converge then, then the maintainer would quite like to know about your problem! If it doesn't converge, then its likely that your model is just too flexible for the IRLS process itself. Having tried increasing maxit in gam.control, there are several other possibilities for stabilizing the iteration. It is possible to try (i) setting lower bounds on the smoothing parameters using the min.sp argument of gam: this may or may not change the model being fitted; (ii) reducing the flexibility of the model by reducing the basis dimensions k in the specification of s and te model terms: this obviously changes the model being fitted somewhat; (iii) introduce a small regularization term into the fitting via the irls.reg argument of gam.control: this option obviously changes the nature of the fit somewhat, since parameter estimates are pulled towards zero by doing this.

Usually, a major contributer to fitting difficulties is that the model is a very poor description of the data.

Author(s)

Simon N. Wood simon.wood@r-project.org


GAM P-IRLS estimation with GCV/UBRE smoothness estimation

Description

This is an internal function of package mgcv. It is a modification of the function glm.fit, designed to be called from 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 mgcv or magic. For further information on usage see code for gam. Some regularization of the IRLS weights is also permitted as a way of addressing identifiability related problems (see gam.control). Negative binomial parameter estimation is supported.

The basic idea of estimating smoothing parameters at each step of the P-IRLS is due to Gu (1992), and is termed `performance iteration' or `performance oriented iteration'.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Gu (1992) Cross-validating non-Gaussian data. J. Comput. Graph. Statist. 1:169-179

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. (2004) Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Amer. Statist. Ass. 99:637-686

See Also

gam.fit2, gam, mgcv, magic


P-IRLS GAM estimation with GCV & UBRE derivative calculation

Description

Estimation of GAM smoothing parameters is most stable if optimization of the UBRE or GCV score is outer to the penalized iteratively re-weighted least squares scheme used to estimate the model given smoothing parameters.

These routines estimates a GAM given log smoothing paramaters, and evaluate derivatives of the GCV and UBRE scores of the model with respect to the log smoothing parameters. Calculation of exact derivatives is generally faster than approximating them by finite differencing, as well as generally improving the reliability of GCV/UBRE score minimization.

gam.fit2 evaluates first derivatives, by accumulating them as P-IRLS progresses. gam.fit3 uses a more efficient approach in which the P-IRLS is first run to convergence, and only then are the derivatives evaluated by a separate iteration. gam.fit3 can evaluate second as well as first derivatives.

Not normally called directly, but rather service routines for gam.

Usage

gam.fit2(x, y, sp, S=list(),rS=list(),off, H=NULL,
         weights = rep(1, nobs), start = NULL, etastart = NULL,
         mustart = NULL, offset = rep(0, nobs), family = gaussian(),
         control = gam.control(), intercept = TRUE,deriv=TRUE,
         gamma=1,scale=1,pearson=FALSE,printWarn=TRUE,...)
gam.fit3(x, y, sp, S=list(),rS=list(),off, H=NULL,
         weights = rep(1, nobs), start = NULL, etastart = NULL,
         mustart = NULL, offset = rep(0, nobs), family = gaussian(),
         control = gam.control(), intercept = TRUE,deriv=2,use.svd=TRUE,
         gamma=1,scale=1,printWarn=TRUE,...)

Arguments

x | The model matrix for the GAM.
y | The response variable.
sp | The log smoothing parameters.
S | A list of penalty matrices. Typically penalty matrices contain only a smallish square sub-matrix which is non-zero: this is what is actually stored. off[i] indicates which parameter is the first one penalized by S[[i]].
rS | List of square roots of penalty matrices, each having as few columns as possible, but as many rows as there are parameters.
off | off[i] indicates which parameter S[[i]][1,1] relates to.
H | The fixed penalty matrix for the model.
weights | prior weights for fitting.
start | optional starting parameter guesses.
etastart | optional starting values for the linear predictor.
mustart | optional starting values for the mean.
offset | the model offset
family | the family - actually this routine would never be called with gaussian()
control | control list as returned from glm.control
intercept | does the model have and intercept, TRUE or FALSE
deriv | Should derivatives of the GCV and UBRE scores be calculated? TRUE or FALSE for gam.fit2, but 0, 1 or 2, indicating the maximum order of differentiation to apply, for gam.fit3.
use.svd | Should the algorithm use SVD (TRUE) or the cheaper QR (FALSE) as the second matrix decomposition of the final derivative/coefficient evaluation method? Only used by gam.fit3.
gamma | The weight given to each degree of freedom in the GCV and UBRE scores can be varied (usually increased) using this parameter.
scale | The scale parameter - needed for the UBRE score.
pearson | The GCV/UBRE score can be based either on the Pearson statistic or the deviance. The latter is generally to be preferred, as it is less prone to severe undersmoothing. Only used by gam.fit2.
printWarn | Set to FALSE to suppress some warnings. Useful in order to ensure that some warnings are only printed if they apply to the final fitted model, rather than an intermediate used in optimization.
... | Other arguments: ignored.

Details

This routine is basically glm.fit with some modifications to allow (i) for quadratic penalties on the log likelihood; (ii) derivatives of the model coefficients with respect to log smoothing parameters to be obtained (by updating alongside the P-IRLS iteration) and (iii) derivatives of the GAM GCV and UBRE scores to be evaluated at convergence.

In addition the routine applies step halving to any step that increases the penalized deviance substantially.

The most costly parts of the calculation are performed by calls to compiled C code (which in turn calls LAPACK routines) in place of the compiled code that would usually perform least squares estimation on the working model in the IRLS iteration.

Estimation of smoothing parameters by optimizing GCV scores obtained at convergence of the P-IRLS iteration was proposed by O'Sullivan et al. (1986), and is here termed `outer' iteration.

Note that use of non-standard families with this routine requires modification of the families as described in fix.family.link.

Author(s)

Simon N. Wood simon.wood@r-project.org

The routine has been modified from glm.fit in R 2.0.1, written by the R core (see glm.fit for further credits).

References

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam.fit, gam, mgcv, magic


Setting GAM fitting method

Description

This is a function of package mgcv which allows selection of the numerical method used to optimize the smoothing parameter estimation criterion for a gam.

It is used to set argument method of gam.

Usage

gam.method(am="magic",gam="outer",outer="newton",gcv="deviance",
           family=NULL)

Arguments

am | Which method to use for a pure additive model (i.e. identity link, gaussian errors). Either "magic" if the Wood (2004) method (magic) is to be used, or "mgcv" if the faster, but less stable and general Wood (2000) method (mgcv) should be used.
gam | Which method to use in the generalized case (i.e. all case other than gaussian with identity link). "perf.magic" for the performance iteration (see details) with magic as the basic estimation engine. "perf.mgcv" for performance iteration with mgcv as the underlying estimation engine. "perf.outer" for magic based performance iteration followed by outer iteration (see details). "outer" for pure outer iteration.
outer | The optimization approach to use to optimize log smoothing parameters by outer iteration."newton" for modified Newton method backed up by steepest descent, based on exact first and second derivatives. "nlm" to use nlm with exact first derivatives to optimize the smoothness selection criterion. "nlm.fd" to use nlm with finite differenced first derivatives (slower and less reliable). "optim" to use the "L-BFGS-B" quasi-Newton method option of routine optim, with exact first derivatives.
gcv | One of "deviance", "GACV" or "pearson", specifying the flavour of GCV to use with outer iteration. "deviance" simply replaces the residual sum of squares term in a GCV score with the deviance, following Hastie and Tibshirani (1990, section 6.9). "GACV" (only available with outer method "newton") uses a varient of Gu and Xiang's (2001) generalized approximate cross validation, modified to deal with arbitrary link-error combinations. "pearson" uses a suggestion from the literature of replacing the RSS in a GCV score with the Pearson statistic: this is a very poor choise (see Wood, 2006, section 4.5.4). This latter option is not available with outer methods "newton" and "nlm".
family | The routine is called by gam to check the supplied method argument. In this circumstance the family argument is passed, to check that it works with the specified method. Negative binomial families only work with performance iteration, and the method is reset to this if necessary.

Details

The performance iteration was suggested by Gu (and is rather similar to the PQL method in generalized linear mixed modelling). At each step of the P-IRLS (penalized iteratively reweighted least squares) iteration, by which a gam is fitted, the smoothing parameters are estimated by GCV or UBRE applied to the working penalized linear modelling problem. In most cases, this process converges and gives smoothness estimates that perform well. It is usually very fast, since smoothing parameters are estimated alongside other model coefficients in a single P-IRLS process.

The performance iteration has two disadvantages. (i) in the presence of co-linearity or concurvity (a frequent problem when spatial smoothers are included in a model with other covariates) then the process can fail to converge. Suppose we start with some coefficient and smoothing parameter estimates, implying a working penalized linear model: the optimal smoothing parameters and coefficients for this working model may in turn imply a working model for which the original estimates are better than the most recent estimates. This sort of effect can prevent convergence.

Secondly it is often possible to find a set of smoothing parameters that result in a lower GCV or UBRE score, for the final working model, than the final score that results from the performance iterations. This is because the performance iteration is only approximately optimizing this score (since optimization is only performed on the working model). The disadvantage here is not that the model with lower score would perform better (it usually doesn't), but rather that it makes model comparison on the basis of GCV/UBRE score rather difficult.

Both disadvantages of performance iteration are surmountable by using what is basically O'Sullivan's (1986) suggestion. Here the P-IRLS scheme is iterated to convergence for a fixed set of smoothing parameters, with an appropriate GCV/UBRE score evaluated at convergence. This score at convergence is optimized in some way. This is termed "outer" optimization, since the optimization is outer to the P-IRLS loop. Outer iteration is slower than performance iteration.

The `appropriate GCV/UBRE' score in the previous paragraph can be defined in one of two ways either (i) the deviance, or (ii) the Pearson statistic can be used in place of the residual sum of squares in the GCV/UBRE score. (ii) makes the GCV/UBRE score correspond to the score for the working linear model at convergence of the P-IRLS, but in practice tends to result in oversmoothing, particularly with low n binomial data, or low mean counts. Hence the default is to use (i).

Several alternative optimisation methods can be used for outer optimization.Usually the fastest and most reliable approach is to use a modified Newton optimizer with exact first and second derivatives, and this is the default. nlm can be used with finite differenced first derivatives. This is not ideal theoretically, since it is possible for the finite difference estimates of derivatives to be very badly in error on rare occasions when the P-IRLS convergence tolerance is close to being matched exactly, so that two components of a finite differenced derivative require different numbers of iterations of P-IRLS in their evaluation. An alternative is provided in which nlm uses numerically exact first derivatives, this is faster and less problematic than the other scheme. An further alternative is to use a quasi-Newton scheme with exact derivtives, based on optim. In practice this usually seems to be slower than the nlm method.

In summary: performance iteration is fast, but can fail to converge. Outer iteration is a little slower, but more reliable. At present only performance iteration is available for negative binomial families.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

Wood, S.N. (2004) Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Amer. Statist. Ass.

http://www.maths.bath.ac.uk/~sw283/

See Also

gam.control gam, gam.fit, glm.control


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 y with an exponential family distribution (e.g. Normal, Poisson, binomial, etc.) using predictor variables x1,x2,..., and let m=E(y). A typical GAM might be:

g(m)=b0 + b1x1 + b2x2 + f1(x3) + f2(x4,x5)

where g is a smooth monotonic `link' function, and f1 and 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 gam.

Specification of the distribution and link function is done using the family argument to gam and works in the same way as for 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 gam formula would be
y~x1+x2+s(x3)+s(x4,x5).
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 s and te, and further discussion of basis dimension choice can be found in choose.k. 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:
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 y is a smooth function of x except at a point x* where the function jumps discontinuously. This model can be written as:

E(y) = b0 + b1 h(x*,x) + f1(x)

where h is a step function jumping from 0 to 1 at x*. The way to fit this model is to create a variable h which is zero for all x less than x* and one otherwise. Then the model formula is:
y~h+s(x).

Another situation that occurs quite often is the one in which we would like to find out if the model:

E(y)=f(x,z)

is really necessary or whether:

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:
y~s(x)+s(z)+s(x,z).
gam automatically generates side conditions to make this model identifiable. You can also estimate `overlapping' models like:
y~s(x,z)+s(z,v).

Sometimes models of the form:

E(y)=b0+f(x)z

need to be estimated (where f is a smooth function, as usual.) The appropriate formula is:
y~z+s(x,by=z)
- the by argument ensures that the smooth function gets multiplied by covariate z, but GAM smooths are centred (average value zero), so the z+ term is needed as well (f is being represented by a constant plus a centred smooth). If we'd wanted:

E(y)=f(x)z

then the appropriate formula would be:
y~z+s(x,by=z)-1.

The 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 fac is a factor with 3 levels 1, 2, 3, and at each level of this factor ther response depends smoothly on a variable x in a manner that is level dependent. Three dummy variables fac.1, fac.2, fac.3, can be generated for the factor (e.g. fac.1<-as.numeric(fac==1)). Then the model formula would be:
y~fac+s(x,by=fac.1)+s(x,by=fac.2)+s(x,by=fac.3).

In the above examples the smooths of more than one covariate have all employed single penalty thin plate regression splines. These isotropic smooths are not alway appropriate: if variables are not naturally `well scaled' relative to each other then it is often preferable to use tensor product smooths, with a wiggliness penalty for each covariate of the term. See te for examples.

WARNING

There are no identifiability checks made between the smooth and parametric parts of a gam formula, and any such lack of identifiability will cause problems if the underlying fitting routine is mgcv (not the default).

Author(s)

Simon N. Wood simon.wood@r-project.org

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)


GAMs with the negative binomial distribution

Description

The gam modelling function is designed to be able to use the negative.binomial and neg.bin families from the MASS library, with or without a known theta parameter. A value for theta must always be passed to these families, but if theta is to be estimated then the passed value is treated as a starting value for estimation.

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

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 theta estimate is updated. This is done by conditioning on all components of the current GCV/Pearson estimator of the scale parameter except theta and then searching for the 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 estimated scale parameter <1. Search limits can be set in gam.control.

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

The negative binomial families can not yet be used with `outer' estimation of smoothing parameters (see gam.method).

Author(s)

Simon N. Wood simon.wood@r-project.org

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)
# unknown theta ...
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=negative.binomial(1))
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)
# known theta example ...
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=negative.binomial(3),scale=1)
plot(b,pages=1)
print(b)
# Now use "sqrt" link available in negative.binomial (but not neg.bin)
set.seed(1)
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)


Minimize GCV or UBRE score of a GAM using 'outer' iteration

Description

Estimation of GAM smoothing parameters is most stable if optimization of the UBRE or GCV score is outer to the penalized iteratively re-weighted least squares scheme used to estimate the model given smoothing parameters.

This routine optimizes a GCV or UBRE score in this way. Basically the GCV or UBRE score is evaluated for each trial set of smoothing parameters by estimating the GAM for those smoothing parameters. The score is minimized w.r.t. the parameters numerically, using optim or nlm. Exact derivatives of the score can be used by fitting with gam.fit2, which improves efficiency and reliability relative to relying solely on finite difference derivatives.

Note that there is a choise between basing GCV/UBRE scores on the deviance or the Pearson statistic: see gam.method.

Not normally called directly, but rather a service routine for gam.

Usage

gam.outer(lsp,fscale,family,control,method,gamma,G,...)

Arguments

lsp | The log smoothing parameters.
fscale | Typical scale of the GCV or UBRE score.
family | the model family.
control | control argument to pass to gam.fit if pure finite differencing is being used.
method | method list reurned from gam.method. This defines the optimization method to use.
gamma | The degree of freedom inflation factor for the GCV/UBRE score.
G | List produced by gam.setup, containing most of what's needed to actually fit GAM.
... | other arguments, typically for passing on to gam.fit2 (ultimately).

Details

Estimation of smoothing parameters by optimizing GCV scores obtained at convergence of the P-IRLS iteration was proposed by O'Sullivan et al. (1986), and is here termed `outer' iteration.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam.fit2, gam, mgcv, magic


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, 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), see gam.method for more detail about the numerical optimization approaches available.

Automatic smoothness selection is unlikely to be successful with few data, particularly with multiple terms to be selected. In addition GCV and UBRE score can occasionally display local minima that can trap the minimisation algorithms. GCV/UBRE scores become constant with changing smoothing parameters at very low or very high smoothng parameters, and on occasion these `flat' regions can be separated from regions of lower score by a small `lip'. This seems to be the most common form of local minimum, but is usually avoidable by avoiding extreme smoothing parameters as starting values in optimization, and by avoiding big jumps in smoothing parameters while optimizing. Never the less, if you are suspicious of smoothing parameter estimates, try changing fit method (see gam.method) and see if the estimates change, or try changing some or all of the smoothing parameters `manually' (argument sp of gam).

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 summary.gam. Candidates for replacement by parametric terms are smooth terms with estimated degrees of freedom close to their minimum possible.

One appealing approach to model selection is via shrinkage. Smooth classes cs.smooth and tprs.smooth (specified by "cs" and "ts" respectively) have smoothness penalties which include a small shrinkage component, so that for large enough smoothing parameters the smooth becomes identically zero. This allows automatic smoothing parameter selection methods to effectively remove the term from the model altogether. The shrinkage component of the penalty is set at a level that usually makes negligable contribution to the penalization of the model, only becoming effective when the term is effectively `completely smooth' according to the conventional penalty.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

Examples

## an example of GCV based model selection
library(mgcv)
set.seed(0)
n<-400;sig<-2
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
x4 <- runif(n, 0, 1);x5 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sig)
y <- f + e
## Note the increased gamma parameter below to favour
## slightly smoother models (See Kim & Gu, 2004, JRSSB)...
b<-gam(y~s(x0,bs="ts")+s(x1,bs="ts")+s(x2,bs="ts")+
   s(x3,bs="ts")+s(x4,bs="ts")+s(x5,bs="ts"),gamma=1.4)
summary(b)
plot(b,pages=1)


Generalized additive model set up

Description

This is an internal function of package mgcv. It is called by 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 gam.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam, gamm, magic, mgcv


Identifiability side conditions for a GAM

Description

GAM formulae with repeated variables only correspond to identifiable models given some side conditions. This routine works out appropriate side conditions, based on zeroing redundant parameters. It is called from gam.setup and is not intended to be called by users.

The method identifies nested and repeated variables by their names, but evaluates what constraints to impose, numerically. Constraints are always applied to smooths of more variables in preference to smooths of fewer variables. The numerical approach allows appropriate constraints to be applied to models constructed using any smooths, including user defined smooths.

Usage

gam.side(sm,tol=.Machine$double.eps^.5)

Arguments

sm | A list of smooth objects as returned by smooth.construct.
tol | The tolernce to use when assessing linear dependence of smooths.

Details

Models such as y~s(x)+s(z)+s(x,z) can be estimated by gam, but require identifiability constraints to be applied, to make them identifiable. This routine does this, effectively setting redundant parameters to zero. When the redundancy is between smooths of lower and higher numbers of variables, the constraint is always applied to the smooth of the higher number of variables.

Dependent smooths are identified symbolically, but which constraints are needed to ensure identifiability of these smooths is determined numerically, using fixDependence. This makes the routine rather general, and not dependent on any particular basis.

Value

A list of smooths, with model matrices and penalty matrices adjusted to automatically impose the required constraints. Any smooth that has been modified will have an attribute "del.index", listing the columns of its model matrix that were deleted. This index is used in the creation of prediction matrices for the term.

Author(s)

Simon N. Wood simon.wood@r-project.org

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
## a fully nested tensor product example
b<-gam(y~s(x,bs="cr",k=6)+s(z,bs="cr",k=6)+te(x,z,k=6))
plot(b)
par(old.par)
rm(list=c("f","x0","x1","x2","x","z","y","b","test1","n","sig2","pi","e"))


Generalized Additive Mixed Models

Description

Fits the specified generalized additive mixed model (GAMM) to data, by a call to lme in the normal errors identity link case, or by a call to glmmPQL from the MASS library otherwise. In the latter case estimates are only approximately MLEs. The routine is typically slower than gam, and not quite as numerically robust.

Smooths are specified as in a call to gam as part of the fixed effects model formula, but the wiggly components of the smooth are treated as random effects. The random effects structures and correlation structures availabale for lme are used to specify other random effects and correlations.

It is assumed that the random effects and correlation structures are employed primarily to model residual correlation in the data and that the prime interest is in inference about the terms in the fixed effects model formula including the smooths. For this reason the routine calculates a posterior covariance matrix for the coefficients of all the terms in the fixed effects formula, including the smooths.

To use this function effectively it helps to be quite familiar with the use of gam and lme.

Usage

gamm(formula,random=NULL,correlation=NULL,family=gaussian(),
data=list(),weights=NULL,subset=NULL,na.action,knots=NULL,
control=nlme::lmeControl(niterEM=0,optimMethod="L-BFGS-B"),
niterPQL=20,verbosePQL=TRUE,method="ML",...)

Arguments

formula A GAM formula (see also 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 y ~ . is not allowed). Smooth terms are specified by expressions of the form:
s(var1,var2,...,k=12,fx=FALSE,bs="tp",by=a.var) where var1, var2, etc. are the covariates which the smooth is a function of and k is the dimension of the basis used to represent the smooth term. If k is not specified then k=10*3^(d-1) is used where d is the number of covariates for this term. fx is used to indicate whether or not this term has a fixed muber of degrees of freedom (fx=FALSE to select d.f. by GCV/UBRE). bs indicates the basis to use: see s for full details, but note that the default "tp" can be slow with large data sets. Tensor product smooths are specified using te terms.
random | The (optional) random effects structure as specified in a call to lme: only the list form is allowed, to facilitate manipulation of the random effects structure within gamm in order to deal with smooth terms. See example below.
correlation | An optional corStruct object (see corClasses) as used to define correlation structures in lme. Any grouping factors in the formula for this object are assumed to be nested within any random effect grouping factors, without the need to make this explicit in the formula (this is slightly different to the behaviour of lme). See examples below.
family | A family as used in a call to glm or gam. The default gaussian with identity link causes gamm to fit by a direct call to lme procided there is no offset term, otherwise glmmPQL from the MASS library is used.
data | A data frame containing the model response variable and covariates required by the formula. By default the variables are taken from environment(formula), typically the environment from which gamm is called.
weights | In the generalized case, weights with the same meaning as glm weights. An lme type weights argument may only be used in the identity link gaussian case, with no offset (see documentation for lme for details of how to use such an argument).
subset | an optional vector specifying a subset of observations to be used in the fitting process.
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'.
knots | this is an optional list containing user specified knot values to be used for basis construction. For the cr basis the user simply supplies the knots to be used, and there must be the same number as the basis dimension, k, for the smooth concerned. For the tp basis knots has two uses. Firstly, for large datasets the calculation of the tp basis can be time-consuming. The user can retain most of the advantages of the t.p.r.s. approach by supplying a reduced set of covariate values from which to obtain the basis - typically the number of covariate values used will be substantially smaller than the number of data, and substantially larger than the basis dimension, k. 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, 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.
control | A list of fit control parameters for lme returned by lmeControl. Note the default setting for the number of EM iterations used by lme: smooths are set up using custom pdMat classes, which are currently not supported by the EM iteration code. Only increase this number if you want to perturb the starting values used in model fitting (usually to worse values!). The optimMethod option is only used if your version of R does not have the nlminb optimizer function.
niterPQL | Maximum number of PQL iterations (if any).
verbosePQL | Should PQL report its progress as it goes along?
method | Which of "ML" or "REML" to use in the Gaussian additive mixed model case when lme is called directly. Ignored in the generalized case (or if the model has an offset), in which case glmmPQL is used.
... further arguments for passing on e.g. to lme

Details

The Bayesian model of spline smoothing introduced by Wahba (1983) and Silverman (1985) opens up the possibility of estimating the degree of smoothness of terms in a generalized additive model as variances of the wiggly components of the smooth terms treated as random effects. Several authors have recognised this (see Wang 1998; Ruppert, Wand and Carroll, 2003) and in the normal errors, identity link case estimation can be performed using general linear mixed effects modelling software such as lme. In the generalized case only approximate inference is so far available, for example using the Penalized Quasi-Likelihood approach of Breslow and Clayton (1993) as implemented in glmmPQL by Venables and Ripley (2002). One advantage of this approach is that it allows correlated errors to be dealt with via random effects or the correlation structures available in the nlme library.

Some brief details of how GAMs are represented as mixed models and estimated using lme or glmmPQL in gamm can be found in Wood (2004a,b). In addition gamm obtains a posterior covariance matrix for the parameters of all the fixed effects and the smooth terms. The approach is similar to that described in (Lin & Zhang, 1999) - the covariance matrix of the data (or pseudodata in the generalized case) implied by the weights, correlation and random effects structure is obtained, based on the estimates of the parameters of these terms and this is used to obtain the posterior covariance matrix of the fixed and smooth effects.

The bases used to represent smooth terms are the same as those used in gam.

In the event of lme convergence failures, consider modifying option(mgcv.vc.logrange): reducing it helps to remove indefiniteness in the likelihood, if that is the problem, but too large a reduction can force over or undersmoothing. See notExp2 for more information on this option. Failing that, you can try increasing the niterEM option in control: this will perturb the starting values used in fitting, but usually to values with lower likelihood! Note that this version of gamm works best with R 2.2.0 or above and nlme, 3.1-62 and above, since these use an improved optimizer.

Value

Returns a list with two items:

gam | an object of class gam, less information relating to GCV/UBRE model selection. At present this contains enough information to use predict, summary and print methods and vis.gam, but not to use e.g. the anova method function to compare models.
lme | the fitted model object returned by lme or glmmPQL. Note that the model formulae and grouping structures may appear to be rather bizarre, because of the manner in which the GAMM is split up and the calls to lme and glmmPQL are constructed.

WARNINGS

This version of gamm works best with nlme 3.1-62 and above and R 2.2.0 and above. This is because it is designed to work with the optimizers used from these versions onwards, rather than the earlier default optimizers.

The routine will be very slow and memory intensive if correlation structures are used for the very large groups of data. e.g. attempting to run the spatial example in the examples section with many 1000's of data is definitely not recommended: often the correlations should only apply within clusters that can be defined by a grouping factor, and provided these clusters do not get too huge then fitting is usually possible.

Models must contain at least one random effect: either a smooth with non-zero smoothing parameter, or a random effect specified in argument random.

Models like s(z)+s(x)+s(x,z) are not currently supported.

gamm is not as numerically stable as gam: an lme call will occasionally fail. See details section for suggestions.

gamm is usually much slower than gam, and on some platforms you may need to increase the memory available to R in order to use it with large data sets (see mem.limits).

Note that the weights returned in the fitted GAM object are dummy, and not those used by the PQL iteration: this makes partial residual plots look odd.

Note that the gam object part of the returned object is not complete in the sense of having all the elements defined in gamObject and does not inherit from glm: hence e.g. multi-model anova calls will not work.

The parameterization used for the smoothing parameters in gamm, bounds them above and below by an effective infinity and effective zero. See notExp2 for details of how to change this.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Breslow, N. E. and Clayton, D. G. (1993) Approximate inference in generalized linear mixed models. Journal of the American Statistical Association 88, 9-25.

Lin, X and Zhang, D. (1999) Inference in generalized additive mixed models by using smoothing splines. JRSSB. 55(2):381-400

Pinheiro J.C. and Bates, D.M. (2000) Mixed effects Models in S and S-PLUS. Springer

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

Silverman, B.W. (1985) Some aspects of the spline smoothing approach to nonparametric regression. JRSSB 47:1-52

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth edition. Springer.

Wahba, G. (1983) Bayesian confidence intervals for the cross validated smoothing spline. JRSSB 45:133-150

Wood, S.N. (2004a) Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of the American Statistical Association. 99:673-686

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

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for Generalized Additive Mixed Models. Biometrics.

Wang, Y. (1998) Mixed effects smoothing spline analysis of variance. J.R. Statist. Soc. B 60, 159-174

http://www.maths.bath.ac.uk/~sw283/

See Also

magic for an alternative for correlated data, te, s, predict.gam, plot.gam, summary.gam, gam.neg.bin, vis.gam,pdTens,gamm.setup

Examples

library(mgcv)
## simple examples using gamm as alternative to gam
set.seed(0)
n <- 400
sig <- 2
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sig)
y <- f + e
b <- gamm(y~s(x0)+s(x1)+s(x2)+s(x3))
plot(b$gam,pages=1)
summary(b$lme) # details of underlying lme fit
summary(b$gam) # gam style summary of fitted model
anova(b$gam)

b <- gamm(y~te(x0,x1)+s(x2)+s(x3))
op <- par(mfrow=c(2,2))
plot(b$gam)
par(op)

## Add a factor to the linear predictor, to be modelled as random
fac <- rep(1:4,n/4)
f <- f + fac*3
fac<-as.factor(fac)

g<-exp(f/5)
y<-rpois(rep(1,n),g)
b2<-gamm(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,random=list(fac=~1))
plot(b2$gam,pages=1)

## now an example with autocorrelated errors....
x <- 0:(n-1)/(n-1)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10-1.396
e <- rnorm(n,0,sig)
for (i in 2:n) e[i] <- 0.6*e[i-1] + e[i]
y <- f + e
op <- par(mfrow=c(2,2))
b <- gamm(y~s(x,k=20),correlation=corAR1())
plot(b$gam);lines(x,f-mean(f),col=2)
b <- gamm(y~s(x,k=20))
plot(b$gam);lines(x,f-mean(f),col=2)
b <- gam(y~s(x,k=20))
plot(b);lines(x,f-mean(f),col=2)

## more complicated autocorrelation example - AR errors
## only within groups defined by `fac'
e <- rnorm(n,0,sig)
for (i in 2:n) e[i] <- 0.6*e[i-1]*(fac[i-1]==fac[i]) + e[i]
y <- f + e
b <- gamm(y~s(x,k=20),correlation=corAR1(form=~1|fac))
plot(b$gam);lines(x,f-mean(f),col=2)
par(op)

## more complex situation with nested random effects and within
## group correlation

set.seed(0)
n.g <- 10
n<-n.g*10*4
sig <- 2
## simulate smooth part
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
## simulate nested random effects....
fa <- as.factor(rep(1:10,rep(4*n.g,10)))
ra <- rep(rnorm(10),rep(4*n.g,10))
fb <- as.factor(rep(rep(1:4,rep(n.g,4)),10))
rb <- rep(rnorm(4),rep(n.g,4))
for (i in 1:9) rb <- c(rb,rep(rnorm(4),rep(n.g,4)))
## simulate auto-correlated errors within groups
e<-array(0,0)
for (i in 1:40) {
eg <- rnorm(n.g, 0, sig)
for (j in 2:n.g) eg[j] <- eg[j-1]*0.6+ eg[j]
e<-c(e,eg)
}
y <- f + ra + rb + e
dat<-data.frame(y=y,x0=x0,x1=x1,x2=x2,x3=x3,fa=fa,fb=fb)
## fit model ....
b <- gamm(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+
  s(x3,bs="cr"),data=dat,random=list(fa=~1,fb=~1),
  correlation=corAR1())
plot(b$gam,pages=1)

## and a "spatial" example
library(nlme);set.seed(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<-200
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)  # true function
f <- test1(x,z)  # true expectation of response
## Now simulate correlated errors...
cstr <- corGaus(.1,form = ~x+z)
cstr <- Initialize(cstr,data.frame(x=x,z=z))
V <- corMatrix(cstr) # correlation matrix for data
Cv <- chol(V)
e <- t(Cv) %*% rnorm(n)*0.05 # correlated errors
## next add correlated simulated errors to expected values
y <- f + e ## ... to produce response
b<- gamm(y~s(x,z,k=50),correlation=corGaus(.1,form=~x+z))
plot(b$gam) # gamm fit accounting for correlation
# overfits when correlation ignored.....
b1 <- gamm(y~s(x,z,k=50));plot(b1$gam)
b2 <- gam(y~s(x,z,k=50));plot(b2)
par(old.par)


Generalized additive mixed model set up

Description

This is an internal function of package mgcv. It is called by gamm to set up a generalized additive mixed model in a form suitable for fitting by calls to glmmPQL from the MASS library or lme from the nlme library. The main task is the representation of the smooth terms as random effects.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Wood, S.N. (2004) Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of the American Statistical Association. 99:673-686

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for Generalized Additive Mixed Models. Biometrics

http://www.maths.bath.ac.uk/~sw283/

See Also

gamm


Fitted gam object

Description

A fitted GAM object returned by function gam and of class "gam" inheriting from classes "glm" and "lm". Method functions anova, logLik, influence, plot, predict, print, residuals and summary exist for this class.

All compulsory elements of "glm" and "lm" objects are present, but the fitting method for a GAM is different to a linear model or GLM, so that the elements relating to the QR decomposition of the model matrix are absent.

Value

A gam object has the following elements:

aic | AIC of the fitted model: bear in mind that the degrees of freedom used to calculate this are the effective degrees of freedom of the model, and the likelihood is evaluated at the maximum of the penalized likelihood in most cases, not at the MLE.
assign | Array whose elements indicate which model term (listed in pterms) each parameter relates to: applies only to non-smooth terms.
boundary | did parameters end up at boundary of parameter space?
call | the matched call (allows update to be used with gam objects, for example).
coefficients | the coefficients of the fitted model. Parametric coefficients are first, followed by coefficients for each spline term in turn.
control | the gam control list used in the fit.
converged | indicates whether or not the iterative fitting method converged.
data | the original supplied data argument (for class "glm" compatibility).
deviance | model deviance (not penalized deviance).
df.null | null degrees of freedom.
df.residual | effective residual degrees of freedom of the model.
edf | estimated degrees of freedom for each model parameter. Penalization means that many of these are less than 1.
family | family object specifying distribution and link used.
fit.method | Character string describing the multiple GCV/UBRE smoothing parameter estimation method used.
fitted.values | fitted model predictions of expected value for each datum.
formula | the model formula.
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.
gcv.ubre | The minimized GCV or UBRE score.
hat | array of elements from the leading diagonal of the `hat' (or `influence') matrix. Same length as response data vector.
iter | number of iterations of P-IRLS taken to get convergence.
linear.predictors | fitted model prediction of link function of expected value for each datum.
method | One of "GCV" or "UBRE", depending on the fitting criterion used.
mgcv.conv A list of convergence diagnostics relating to the "mgcv" or "magic" parts of smoothing parameter estimation - this will not be very meaningful for pure "outer" estimation of smoothing parameters. mgcv.conv differs for method "magic" and "mgcv". Here is the "mgcv" version:
score
corresponding to edf, an array of GCV or UBRE scores for the model given the final estimated relative smoothing parameters.

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

h
the second derivatives corresponding to g above - i.e. the leading diagonal of the Hessian.

e
the eigen-values of the Hessian. All non-negative indicates a positive definite Hessian.

iter
the number of iterations taken.

in.ok
TRUE if the second smoothing parameter guess improved the GCV/UBRE score.

step.fail
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 "magic" the items are:
full.rank
The apparent rank of the problem given the model matrix and constraints.

rank
The numerical rank of the problem.

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

hess.pos.def
Was the hessian of the GCV/UBRE score positive definite at smoothing parameter estimation convergence?

iter
How many iterations were required to find the smoothing parameters?

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

rms.grad
root mean square of the gradient of the GCV/UBRE score at convergence.

min.edf | Minimum possible degrees of freedom for whole model.
model | model frame containing all variables needed in original model fit.
na.action | The na.action used in fitting.
nsdf | number of parametric, non-smooth, model terms including the intercept.
null.deviance | deviance for single parameter model.
offset | model offset.
outer.info | If `outer' iteration has been used to fit the model (see gam.method) then this is present and contains whatever was returned by the optimization routine used (currently nlm or optim).
prior.weights | prior weights on observations.
pterms | terms object for strictly parametric part of model.
rank | apparent rank of fitted model.
residuals | the working residuals for the fitted model.
sig2 | estimated or supplied variance/scale parameter.
smooth | list of smooth objects, containing the basis information for each term in the model formula in the order in which they appear. These smooth objects are what gets returned by the smooth.construct objects.
sp | smoothing parameter for each smooth.
terms | terms object of model model frame.
Vp | estimated covariance matrix for the parameters. This is a Bayesian posterior covariance matrix that results from adopting a particular Bayesian model of the smoothing process. Paricularly useful for creating credible/confidence intervals.
Ve | frequentist estimated covariance matrix for the parameter estimators. Particularly useful for testing whether terms are zero. Not so useful for CI's as smooths are usually biased.
weights | final weights used in IRLS iteration.
y | response data.

WARNINGS

This model object is different to that described in Chambers and Hastie (1993) in order to allow smoothing parameter estimation etc.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Key References on this implementation:

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

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

Wood, S.N. (in press) Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Amer. Statist. Ass.

Wood, S.N. (2004) On confidence intervals for GAMs based on penalized regression splines. Technical Report 04-12 Department of Statistics, University of Glasgow.

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for generalized additive mixed models. Biometrics.

Key Reference on GAMs and related models:

Hastie (1993) in Chambers and Hastie (1993) Statistical Models in S. Chapman and Hall.

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

Wahba (1990) Spline Models of Observational Data. SIAM

See Also

gam


Get named variable or evaluate expression from list or data.frame

Description

This routine takes a text string and a data frame or list. It first sees if the string is the name of a variable in the data frame/ list. If it is then the value of this variable is returned. Otherwise the routine tries to evaluate the expression within the data.frame/list (but nowhere else) and if successful returns the result. If neither step works then NULL is returned. The routine is useful for processing gam formulae.

Usage

 get.var(txt,data)

Arguments

txt | a text string which is either the name of a variable in data or when parsed is an expression that can be evaluated in data. It can also be neither in which case the function returns NULL.
data | A data frame or list.

Value

The evaluated variable, or NULL

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

gam

Examples

y <- 1:4;dat<-data.frame(x=5:10)
get.var("x",dat)
get.var("y",dat)
get.var("x==6",dat)


Extract the diagonal of the influence/hat matrix for a GAM

Description

Extracts the leading diagonal of the influence matrix (hat matrix) of a fitted gam object.

Usage

## S3 method for class 'gam':
influence(model,...)

Arguments

model | fitted model objects of class gam as produced by gam().
... | un-used in this case

Details

Simply extracts hat array from fitted model. (More may follow!)

Value

An array (see above).

Author(s)

Simon N. Wood simon.wood@r-project.org

See Also

gam


Starting values for multiple smoothing parameter estimation

Description

Finds initial smoothing parameter guesses for multiple smoothing parameter estimation. The idea is to find values such that the estimated degrees of freedom per penalized parameter should be well away from 0 and 1 for each penalized parameter, thus ensuring that the values are in a region of parameter space where the smoothing parameter estimation criterion is varying substantially with smoothing parameter value.

Usage

initial.sp(X,S,off,expensive=FALSE)

Arguments

X | is the model matrix.
S | is a list of of penalty matrices. 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 S[[i]] occupies element off[i], off[i] of the ith penalty matrix. Each S[[i]] must be positive semi-definite.
off | is an array indicating the first parameter in the parameter vector that is penalized by the penalty involving S[[i]].
expensive | if TRUE then the overall amount of smoothing is adjusted so that the average degrees of freedom per penalized parameter is exactly 0.5: this is numerically costly.

Details

Basically uses a crude approximation to the estimated degrees of freedom per model coefficient, to try and find smoothing parameters which bound these e.d.f.'s away from 0 and 1.

Usually only called by magic and gam.

Value

An array of initial smoothing parameter estimates.

Author(s)

Simon N. Wood simon.wood@r-project.org

See Also

magic, gam.outer, gam,


Interpret a GAM formula

Description

This is an internal function of package mgcv. It is a service routine for gam which splits off the strictly parametric part of the model formula, returning it as a formula, and interprets the smooth parts of the model formula.

Not normally called directly.

Usage

interpret.gam(gf)

Arguments

gf A GAM formula as supplied to gam or gamm.

Value

An object of class split.gam.formula with the following items:

pf | A model formula for the strictly parametric part of the model.
pfok | TRUE if there is a pf formula.
smooth.spec | A list of class xx.smooth.spec objects where xx depends on the basis specified for the term. (These can be passed to smooth constructor method functions to actually set up penalties and bases.)
full.formula | An expanded version of the model formula in which the options are fully expanded, and the options do not depend on variables which might not be available later.
fake.formula | A formula suitable for use in evaluating a model frame.
response | Name of the response variable.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam gamm


Extract the log likelihood for a fitted GAM

Description

Function to extract the log-likelihood for a fitted gam model (note that the models are usually fitted by penalized likelihood maximization).

Usage

logLik.gam(object,...)

Arguments

object | fitted model objects of class gam as produced by gam().
... | un-used in this case

Details

Modification of logLik.glm which corrects the degrees of freedom for use with gam objects.

The function is provided so that AIC functions correctly with gam objects, and uses the appropriate degrees of freedom (accounting for penalization). Note, when using AIC for penalized models, that the degrees of freedom are the effective degrees of freedom and not the number of parameters, and the model maximizes the penalized likelihood, not the actual likelihood! This seems to be reasonably well founded for the known scale parameter case (see Hastie and Tibshirani, 1990, section 6.8.3), and in fact in this case the default smoothing parameter estimation criterion is effectively this modified AIC.

Value

Standard logLik object: see logLik.

Author(s)

Simon N. Wood simon.wood@r-project.org based directly on logLik.glm

References

Hastie and Tibshirani, 1990, Generalized Additive Models.

See Also

AIC


Stable Multiple Smoothing Parameter Estimation by GCV or UBRE, with optional fixed penalty

magic(mgcv) R Documentation

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.

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),extra.rss=0,n.score=length(y))

Arguments

y | is the response data vector.
X | is the model matrix.
sp | is the array of smoothing parameters multiplying the penalty matrices stored in 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.
S | is a list of of penalty matrices. 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 S[[i]] occupies element off[i], off[i] of the ith penalty matrix. Each S[[i]] must be positive semi-definite.
off | is an array indicating the first parameter in the parameter vector that is penalized by the penalty involving S[[i]].
rank | is an array specifying the ranks of the penalties. This is useful, but not essential, for forming square roots of the penalty matrices.
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.
C | is the optional matrix specifying any linear equality constraints on the fitting problem. If b is the parameter vector then the parameters are forced to satisfy Cb=0.
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 y, specifically V_y^{-1}=w'w. If w is an array then it is taken as the diagonal of this matrix, or simply the weight for each element of y. See below for an example using this.
gamma | is an inflation factor for the model degrees of freedom in the GCV or UBRE score.
scale | is the scale parameter for use with UBRE.
gcv | should be set to TRUE if GCV is to be used, FALSE for UBRE.
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.
control | is a list of iteration control constants with the following elements: maxit The maximum number of iterations of the magic algorithm to allow. tol The tolerance to use in judging convergence. step.half If a trial step fails then the method tries halving it up to a maximum of step.half times. rank.tol is a constant used to test for numerical rank deficiency of the problem. Basically any singular value less than rank_tol multiplied by the largest singular value of the problem is set to zero.
extra.rss | is a constant to be added to the residual sum of squares (squared norm) term in the calculation of the GCV, UBRE and scale parameter estimate. In conjuction with n.score, this is useful for certain methods for dealing with very large data sets.
n.score | number to use as the number of data in GCV/UBRE score calculation: usually the actual number of data, but there are methods for dealing with very large datasets that change this.

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:

minimise || W (Xb-y) ||^2 + b'Hb + theta_1 b'S_1 b + theta_2 b'S_2 b + . . .

possibly subject to constraints Cb=0. X is a design matrix, b a parameter vector, y a data vector, W a weight matrix, S_i a positive semi-definite matrix of coefficients defining the ith penalty with associated smoothing parameter theta_i, H is the positive semi-definite offset penalty matrix and C a matrix of coefficients defining any linear equality constraints on the problem. X need not be of full column rank.

The theta_i are chosen to minimize either the GCV score:

V_g = n ||W(y-Ay)||^2/[tr(I - g A)]^2

or the UBRE score:

V_u =||W(y-Ay||^2/n - 2 s tr(I - g A)/n + s

where g is gamma the inflation factor for degrees of freedom (usually set to 1) and s is scale, the scale parameter. 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 A.

The method operates by Newton or steepest descent updates of the logs of the 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 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 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 b parameters are easily found.)

The method is coded in C with matrix factorizations performed using LINPACK and LAPACK routines.

Value

The function returns a list with the following items:

b | The best fit parameters given the estimated smoothing parameters.
scale | the estimated (GCV) or supplied (UBRE) scale parameter.
score | the minimized GCV or UBRE score.
sp | an array of the estimated smoothing parameters.
rV | a factored form of the parameter covariance matrix. The (Bayesian) covariance matrix of the parametes b is given by rV%*%t(rV)*scale.
gcv.info | is a list of information about the performance of the method with the following elements: full.rank The apparent rank of the problem: number of parameters less number of equality constraints. rank The estimated actual rank of the problem (at the final iteration of the method). fully.converged is TRUE if the method converged by satisfying the convergence criteria, and FALSE if it coverged by failing to decrease the score along the search direction. hess.pos.def is TRUE if the hessian of the UBRE or GCV score was positive definite at convergence. iter is the number of Newton/Steepest descent iterations taken. score.calls is the number of times that the GCV/UBRE score had to be evaluated. 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 magic.post.proc.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

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

http://www.maths.bath.ac.uk/~sw283/

See Also

magic.post.proc, mgcv, gam,

Examples

## Use `magic' for a standard additive model fit ...
   library(mgcv)
   set.seed(1);n <- 400;sig2 <- 4
   x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
   x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
   f <- 2 * sin(pi * x0)
   f <- f + exp(2 * x1) - 3.75887
   f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10
   e <- rnorm(n, 0, sqrt(sig2))
   y <- f + e
## set up additive model
   G <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),fit=FALSE)
## fit using magic
   mgfit <- magic(G$y,G$X,G$sp,G$S,G$off,G$rank,C=G$C)
## and fit using gam as consistency check
   b <- gam(G=G)
   mgfit$sp;b$sp  # compare smoothing parameter estimates
   edf <- magic.post.proc(G$X,mgfit,G$w)$edf  # extract e.d.f. per parameter
## get termwise e.d.f.s
   twedf <- 0;for (i in 1:4) twedf[i] <- sum(edf[((i-1)*10+1):(i*10)])
   twedf;b$edf  # compare

## Now a correlated data example ...
    library(nlme)
## simulate truth
    set.seed(1);n<-400;sig<-2
    x <- 0:(n-1)/(n-1)
    f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10-1.396
## produce scaled covariance matrix for AR1 errors...
    V <- corMatrix(Initialize(corAR1(.6),data.frame(x=x)))
    Cv <- chol(V)  # t(Cv)
## Simulate AR1 errors ...
    e <- t(Cv)%*%rnorm(n,0,sig) # so cov(e) = V * sig^2
## Observe truth + AR1 errors
    y <- f + e
## GAM ignoring correlation
    par(mfrow=c(1,2))
    b <- gam(y~s(x,k=20))
    plot(b);lines(x,f-mean(f),col=2);title("Ignoring correlation")
## Fit smooth, taking account of *known* correlation...
    w <- solve(t(Cv)) # V^{-1} = w'w
    ## Use `gam' to set up model for fitting...
    G <- gam(y~s(x,k=20),fit=FALSE)
    ## fit using magic, with weight *matrix*
    mgfit <- magic(G$y,G$X,G$sp,G$S,G$off,G$rank,C=G$C,w=w)
## Modify previous gam object using new fit, for plotting...
    mg.stuff <- magic.post.proc(G$X,mgfit,w)
    b$edf <- mg.stuff$edf;b$Vp <- mg.stuff$Vb
    b$coefficients <- mgfit$b
    plot(b);lines(x,f-mean(f),col=2);title("Known correlation")


Auxilliary information from magic fit

Description

Obtains Bayesian parameter covariance matrix, frequentist parameter estimator covariance matrix, estimated degrees of freedom for each parameter and leading diagonal of influence/hat matrix, for a penalized regression estimated by magic.

Usage

magic.post.proc(X,object,w)

Arguments

X | is the model matrix.
object | is the list returned by magic after fitting the model with model matrix X.
w | is the weight vector used in fitting, or the weight matrix used in fitting (i.e. supplied to magic, if one was.) t(w)%*%w should typically give the inverse of the covariance matrix of the response data supplied to magic.

Details

object contains rV (V, say), and scale (s, say) which can be used to obtain the require quantities as follows. The Bayesian covariance matrix of the parameters is VV's. The vector of estimated degrees of freedom for each parameter is the leading diagonal of VV'X'W'WX where W is either the weight matrix w or the matrix diag(w). The hat/influence matrix is given by WXVV'X'W' .

The frequentist parameter estimator covariance matrix is VV'X'W'WXVV's: it is useful for testing terms for equality to zero.

Value

A list with three items:

Vb | the Bayesian covariance matrix of the model parameters.
Ve | the frequentist covariance matrix for the parameter estimators.
hat | the leading diagonal of the hat (influence) matrix.
edf | the array giving the estimated degrees of freedom associated with each parameter.

Author(s)

Simon N. Wood simon.wood@r-project.org

See Also

magic


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.

For a listing of all routines in the mgcv package type:
library(help="mgcv")

Usage

mgcv(y,X,sp,S,off,C=NULL,w=rep(1,length(y)),H=NULL,
     scale=1,gcv=TRUE,control=mgcv.control())

Arguments

y | The response data vector.
X | The design matrix for the problem, note that ncol(X) must give the number of model parameters, while nrow(X) should give the number of data.
sp | An array of smoothing parameters. If control$fixed==TRUE then these are taken as being the smoothing parameters. Otherwise any positive values are assumed to be initial estimates and negative values to signal auto-initialization.
S | A list of penalty matrices. Only the smallest square block containing all non-zero matrix elements is actually stored, and off[i] indicates the element of the parameter vector that S[[i]][1,1] relates to.
off | Offset values indicating where in the overall parameter a particular stored penalty starts operating. For example if p is the model parameter vector and k=nrow(S[[i]])-1, then the ith penalty is given by t(p[off[i]:(off[i]+k)])%*%S[[i]]%*%p[off[i]:(off[i]+k)].
C | Matrix containing any linear equality constraints on the problem (i.e. C in Cp=0).
w | A vector of weights for the data (often proportional to the reciprocal of the standard deviation of y).
H | A single fixed penalty matrix to be used in place of the multiple penalty matrices in S. mgcv cannot mix fixed and estimated penalties.
scale | This is the known scale parameter/error variance to use with UBRE. Note that it is assumed that the variance of y_i is given by scale/w_i.
gcv | If gcv is TRUE then smoothing parameters are estimated by GCV, otherwise UBRE is used.
control | A list of control options returned by mgcv.control.

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:

minimise || W (Xp-y) ||^2 rho + lambda_1 p'S_1 p + lambda_1 p'S_2 p + . . .

possibly subject to constraints Cp=0. X is a design matrix, p a parameter vector, y a data vector, W a diagonal weight matrix, S_i a positive semi-definite matrix of coefficients defining the ith penalty and C a matrix of coefficients defining any linear equality constraints on the problem. The smoothing parameters are the lambda_i but there is an overall smoothing parameter rho as well. Note that 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 rho with Newton or steepest descent updates of the logs of the lambda_i. Because the GCV/UBRE scores are flat w.r.t. very large or very small 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 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. M$conv provides some convergence diagnostics.

The method is coded in 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

An object is returned with the following elements:

b | The best fit parameters given the estimated smoothing parameters.
scale | The estimated or supplied scale parameter/error variance.
score | The UBRE or GCV score.
sp | The estimated (or supplied) smoothing parameters (lambda_i/rho)
Vb | Estimated covariance matrix of model parameters.
hat | diagonal of the hat/influence matrix.
edf | array of estimated degrees of freedom for each parameter.
info A list of convergence diagnostics, with the following elements:
    edf
    Array of whole model estimated degrees of freedom.
    score
    Array of ubre/gcv scores at the edfs for the final set of relative smoothing parameters.
    g
    the gradient of the GCV/UBRE score w.r.t. the smoothing parameters at termination.
    h
    the second derivatives corresponding to g above - i.e. the leading diagonal of the Hessian.
    e
    the eigenvalues of the Hessian. These should all be non-negative!
    iter
    the number of iterations taken.
    in.ok
    TRUE if the second smoothing parameter guess improved the GCV/UBRE score. (Please report examples where this is FALSE)
    step.fail
    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.

WARNING

The method may not behave well with near column rank deficient X especially in contexts where the weights vary wildly.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam, magic

Examples

library(help="mgcv") # listing of all routines

set.seed(1);n<-400;sig2<-4
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sqrt(sig2))
y <- f + e
# set up additive model
G<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),fit=FALSE)
# fit using mgcv
mgfit<-mgcv(G$y,G$X,G$sp,G$S,G$off,C=G$C)


GAMs with GCV smoothness estimation and GAMMs by REML/PQL

Description

mgcv provides functions for generalized additive modelling and generalized additive mixed modelling. Particular features of the package are facilities for automatic smoothness selection, and the provision of a variety of smooths of more than one variable. User defined smooths are also supported. A Bayesian approach to confidence/credible interval calculation is provided. Lower level routines for generalized ridge regression and penalized linearly constrained least squares are also provided.

Details

mgcv provides generalized additive modelling functions gam, predict.gam and plot.gam, which are very similar in use to the S functions of the same name designed by Trevor Hastie. However the underlying representation and estimation of the models is based on a penalized regression spline approach, with automatic smoothness selection. A number of other functions such as summary.gam and anova.gam are also provided, for extracting information from a fitted gamObject.

Use of gam is much like use of glm, except that within a gam model formula, isotropic smooths of any number of predictors can be specified using s terms, while scale invariant smooths of any number of predictors can be specified using te terms. Estimation is by penalized likelihood or quasi-likelihood maximization, with smoothness selection by GCV or gAIC/ UBRE. See gam, gam.models and gam.selection for some discussion of model specification and selection. For detailed control of fitting see gam.convergence, gam.method and gam.control. For checking and visualization see gam.check, choose.k, vis.gam and plot.gam. While a number of types of smoother are built into the package, it is also extendable with user defined smooths, see p.spline, for example.

A Bayesian approach to smooth modelling is used to derive standard errors on predictions, and hence credible intervals. The Bayesian covariance matrix for the model coefficients is returned in Vp of the gamObject. See predict.gam for examples of how this can be used to obtain credible regions for any quantity derived from the fitted model, either directly, or by direct simulation from the posterior distribution of the model coefficients. Frequentist approximations can be used for hypothesis testing: see anova.gam and summary.gam, but note that the underlying approximations are not always good in this case.

The package also provides a generalized additive mixed modelling function, gamm, based on glmmPQL from the MASS library and lme from the nlme library. gamm is particularly useful for modelling correlated data (i.e. where a simple independence model for the residual variation is inappropriate). In addition, low level routine magic can fit models to data with a known correlation structure.

Some of the underlying GAM fitting methods are available as low level fitting functions: see magic and mgcv. Penalized weighted least squares with linear equality and inequality constraints is provided by pcls.

For a complete list of functions type library(help=mgcv).

Author(s)

Simon Wood <simon.wood@r-project.org>

with contributions and/or help from Kurt Hornik, Mike Lonergan, Henric Nilsson and Brian Ripley.

Maintainer: Simon Wood <simon.wood@r-project.org>

References

Wood, S.N. (2006) Generalized Additive Models: an introduction with R, CRC

Examples

## see examples for gam and gamm


Setting mgcv defaults

Description

This is an internal function of package mgcv which allows control of the numerical options for fitting a generalized ridge regression problem using routine mgcv.

Usage

mgcv.control(conv.tol=1e-7,max.half=20,target.edf=NULL,min.edf=-1)

Arguments

conv.tol | The convergence tolerance.
max.half | successive step halvings are employed if the Newton method and then the steepest descent backup fail to improve the UBRE/GCV score. This is how many to use before giving up.
target.edf | If this is non-null it indicates that cautious optimization should be used, which opts for the local minimum closest to the target model edf if there are multiple local minima in the GCV/UBRE score.
min.edf | Lower bound on the model edf. Useful for avoiding numerical problems at high smoothing parameter values. Negative for none.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

mgcv


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 gam, "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)

Arguments

x | The array of knot locations.
up | If TRUE then the constraints imply increase, if FALSE then decrease.
lower | This specifies the lower bound on the spline unless it is NA in which case no lower bound is imposed.
upper | This specifies the upper bound on the spline unless it is NA in which case no upper bound is imposed.

Details

Consider the natural cubic spline passing through the points: (x_i,p_i), i=1..n. Then it is possible to find a relatively small set of linear constraints on p sufficient to ensure monotonicity (and bounds if required): Ap>=b. Details are given in Wood (1994). This function returns a list containing A and b.

Value

The function returns a list containing constraint matrix A and constraint vector b.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

mgcv pcls

Examples

## see ?pcls


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")

Arguments

A | The positive semi-definite matrix, a square root of which is to be found.
rank | if the rank of the matrix A is known then it should be supplied.
method | "chol" to use pivoted choloeski decompositon, which is fast but tends to over-estimate rank. "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, B with as many columns as the rank of A, and such that A=BB'.

Author(s)

Simon N. Wood simon.wood@r-project.org

Examples

  set.seed(0)
  a <- matrix(runif(24),6,4)
  A <- a%*%t(a) ## A is +ve semi-definite, rank 4
  B <- mroot(A) ## default pivoted choleski method
  tol <- 100*.Machine$double.eps
  chol.err <- max(abs(A-B%*%t(B)));chol.err
  if (chol.err>tol) warning("mroot (chol) suspect")
  B <- mroot(A,method="svd") ## svd method
  svd.err <- max(abs(A-B%*%t(B)));svd.err
  if (svd.err>tol) warning("mroot (svd) suspect")


Obtain a name for a new variable that is not already in use

Description

gamm works by transforming a GAMM into something that can be estimated by lme, but this involves creating new variables, the names of which should not clash with the names of other variables on which the model depends. This simple service routine checks a suggested name against a list of those in use, and if neccesary modifies it so that there is no clash.

Usage

 new.name(proposed,old.names)

Arguments

proposed | a suggested name
old.names | An array of names that must not be duplicated

Value

A name that is not in old.names.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

gamm

Examples

old <- c("a","tuba","is","tubby")
new.name("tubby",old)


Functions for better-than-log positive parameterization

Description

It is common practice in statistical optimization to use log-parameterizations when a parameter ought to be positive. i.e. if an optimization parameter a should be non-negative then we use a=exp(b) and optimize with respect to the unconstrained parameter b. This often works well, but it does imply a rather limited working range for b: using 8 byte doubles, for example, if b's magnitude gets much above 700 then a overflows or underflows. This can cause problems for numerical optimization methods.

notExp is a monotonic function for mapping the real line into the positive real line with much less extreme underflow and overflow behaviour than exp. It is a piece-wise function, but is continuous to second derivative: see the source code for the exact definition, and the example below to see what it looks like.

notLog is the inverse function of notExp.

The major use of these functions was originally to provide more robust pdMat classes for lme for use by gamm. Currently the notExp2 and notLog2 functions are used in their place, as a result of changes to the nlme optimization routines.

Usage

notExp(x)

notLog(x)

Arguments

x Argument array of real numbers (notExp) or positive real numbers (notLog).

Value

An array of function values evaluated at the supplied argument values.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

pdTens, pdIdnot, gamm

Examples

## Illustrate the notExp function:
## less steep than exp, but still monotonic.
x <- -100:100/10
op <- par(mfrow=c(2,2))
plot(x,notExp(x),type="l")
lines(x,exp(x),col=2)
plot(x,log(notExp(x)),type="l")
lines(x,log(exp(x)),col=2) # redundancy intended
x <- x/4
plot(x,notExp(x),type="l")
lines(x,exp(x),col=2)
plot(x,log(notExp(x)),type="l")
lines(x,log(exp(x)),col=2) # redundancy intended
par(op)
range(notLog(notExp(x))-x) # show that inverse works!


Alternative to log parameterization for variance components

Description

notLog2 and notExp2 are alternatives to log and exp or notLog and notExp for re-parameterization of variance parameters. They are used by the pdTens and pdIdnot classes which in turn implement smooths for gamm.

The functions are typically used to ensure that smoothing parameters are positive, but the notExp2 is not monotonic: rather it cycles between `effective zero' and `effective infinity' as its argument changes. The notLog2 is the inverse function of the notExp2 only over an interval centered on zero.

Parameterizations using these functions ensure that estimated smoothing parameters remain positive, but also help to ensure that the likelihood is never indefinite: once a working parameter pushes a smoothing parameter below `effetive zero' or above `effective infinity' the cyclic nature of the notExp2 causes the likelihood to decrease, where otherwise it might simply have flattened.

This parameterization is really just a numerical trick, in order to get lme to fit gamm models, without failing due to indefiniteness. Note in particular that asymptoticresults on the likelihood/REML criterion are not invalidated by the trick, unless parameter estimates end up close to the effective zero or effective infinity: but if this is the case then the asymptotics would also have been invalid for a conventional monotonic parameterization.

This reparameterization was made necessary by some modifications to the underlying optimization method in lme introduced in nlme 3.1-62. It is possible that future releases will return to the notExp parameterization.

Note that you can reset `effective zero' and `effective infinity': see below.

Usage

notExp2(x,d=.Options$mgcv.vc.logrange,b=1/d)

notLog2(x,d=.Options$mgcv.vc.logrange,b=1/d)

Arguments

x | Argument array of real numbers (notExp) or positive real numbers (notLog).
d | the range of notExp2 runs from exp(-d) to exp(d). To change the range used by gamm reset mgcv.vc.logrange using options.
b | determines the period of the cycle of notExp2.

Value

An array of function values evaluated at the supplied argument values.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

pdTens, pdIdnot, gamm

Examples

## Illustrate the notExp2 function:
x <- seq(-50,50,length=1000)
op <- par(mfrow=c(2,2))
plot(x,notExp2(x),type="l")
lines(x,exp(x),col=2)
plot(x,log(notExp2(x)),type="l")
lines(x,log(exp(x)),col=2) # redundancy intended
x <- x/4
plot(x,notExp2(x),type="l")
lines(x,exp(x),col=2)
plot(x,log(notExp2(x)),type="l")
lines(x,log(exp(x)),col=2) # redundancy intended
par(op)


The basis of the space of un-penalized functions for a TPRS

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. null.space.dimension finds the dimension of this space, M, given the number of covariates that the smoother is a function of, d, and the order of the smoothing penalty, m. If m does not satisfy 2m>d then the smallest possible dimension for the null space is found given d and the requirement that the smooth should be visually smooth.

Usage

null.space.dimension(d,m)

Arguments

d | is a positive integer - the number of variables of which the t.p.s. is a function.
m | a non-negative integer giving the order of the penalty functional, or signalling that the default order should be used.

Details

Thin plate splines are only visually smooth if the order of the wiggliness penalty, m, satisfies 2m > d+1. If 2m<d+1 then this routine finds the smallest m giving visual smoothness for the given d, otherwise the supplied m is used. The null space dimension is given by:

M=(m+d+1)!/(d!(m-d)!

which is the value returned.

Value

An integer (array), the null space dimension M.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam

Examples

null.space.dimension(2,0)


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)

Arguments

M is the single list argument to pcls. It should have the following elements:
    y
    The response data vector.
    w
    A vector of weights for the data (often proportional to the reciprocal of the variance).
    X
    The design matrix for the problem, note that ncol(M$X) must give the number of model parameters, while nrow(M$X) should give the number of data.
    C
    Matrix containing any linear equality constraints on the problem (e.g. C in Cp=c). If you have no equality constraints initialize this to a zero by zero matrix. Note that there is no need to supply the vector c, it is defined implicitly by the initial parameter estimates p.
    S
    A list of penalty matrices. S[[i]] is the smallest contiguous matrix including all the non-zero elements of the ith penalty matrix. The first parameter it penalizes is given by off[i]+1 (starting counting at 1).
    off
    Offset values locating the elements of M$S in the correct location within each penalty coefficient matrix. (Zero offset implies starting in first location)
    sp
    An array of smoothing parameter estimates.
    p
    An array of feasible initial parameter estimates - these must satisfy the constraints, but should avoid satisfying the inequality constraints as equality constraints.
    Ain
    Matrix for the inequality constraints A_in p > b.
    bin
    vector in the inequality constraints.

Details

This solves the problem:

minimise || W^0.5 (Xp-y) ||^2 + lambda_1 p'S_1 p + lambda_1 p'S_2 p + . . .

subject to constraints Cp=c and A_in p > b_in, w.r.t. p given the smoothing parameters lambda_i. X is a design matrix, p a parameter vector, y a data vector, W a diagonal weight matrix, S_i a positive semi-definite matrix of coefficients defining the ith penalty and C a matrix of coefficients defining the linear equality constraints on the problem. The smoothing parameters are the lambda_i. Note that X must be of full column rank, at least when projected into the null space of any equality constraints. A_in is a matrix of coefficients defining the inequality constraints, while 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. X'X is not formed explicitly. See Gill et al. 1981.

Value

The function returns an array containing the estimated parameter vector.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

mgcv 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),S=list(),
Ain=matrix(0,1,2),bin=0,C=matrix(0,0,0),sp=array(0,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)

# Penalized example: 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, constraints etc. for monotonic spline....
sm<-smoothCon(s(x,k=10,bs="cr"),dat,knots=NULL)
F<-mono.con(sm$xp);   # get constraints
G<-list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,p=sm$xp,y=y,w=y*0+1)
G$Ain<-F$A;G$bin<-F$b;G$S<-sm$S;G$off<-0

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

fv<-Predict.matrix(sm,data.frame(x=x))%*%p
lines(x,fv,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....
sm<-smoothCon(s(x,k=10,bs="tp"),dat,knots=NULL)
xc<-0:39/39 # points on [0,1]
nc<-length(xc)  # number of constraints
xc<-xc*4-1  # points at which to impose constraints
A0<-Predict.matrix(sm,data.frame(x=xc))
# ... A0
A1<-Predict.matrix(sm,data.frame(x=xc+1e-6))
A<-(A1-A0)/1e-6
# ... approx. constraint matrix (A%*%p is -ve spline gradient at points xc)
G<-list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,y=y,w=y*0+1,S=sm$S,off=0)
G$Ain<-A;    # constraint matrix
G$bin<-rep(0,nc);  # constraint vector
G$p<-rep(0,10);G$p[10]<-0.1
# ... monotonic start params, got by setting coefs of polynomial part
p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

fv2<-Predict.matrix(sm,data.frame(x=x))%*%p
lines(x,fv2,col=3)


Overflow proof pdMat class for multiples of the identity matrix

Description

This set of functions is a modification of the pdMat class pdIdent from library nlme. The modification is to replace the log parameterization used in pdMat with a notLog2 parameterization, since the latter avoids indefiniteness in the likelihood and associated convergence problems: the parameters also relate to variances rather than standard deviations, for consistency with the pdTens class. The functions are particularly useful for working with Generalized Additive Mixed Models where variance parameters/smoothing parameters can be very large or very small, so that overflow or underflow can be a problem.

These functions would not normally be called directly, although unlike the pdTens class it is easy to do so.

Usage

pdIdnot(value = numeric(0), form = NULL,
       nam = NULL, data = sys.frame(sys.parent()))

Arguments

value | Initialization values for parameters. Not normally used.
form | A one sided formula specifying the random effects structure.
nam | a names argument, not normally used with this class.
data | data frame in which to evaluate formula.

Details

Note that while the pdFactor and pdMatrix functions return the inverse of the scaled random effect covariance matrix or its factor, the pdConstruct function is initialised with estimates of the scaled covariance matrix itself.

Value

A class pdIdnot object, or related quantities. See the nlme documentation for further details.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Pinheiro J.C. and Bates, D.M. (2000) Mixed effects Models in S and S-PLUS. Springer

The nlme source code.

http://www.maths.bath.ac.uk/~sw283/

See Also

te, pdTens, notLog2, gamm

Examples

# see gamm


Functions implementing a pdMat class for tensor product smooths

Description

This set of functions implements an nlme library pdMat class to allow tensor product smooths to be estimated by lme as called by gamm. Tensor product smooths have a penalty matrix made up of a weighted sum of penalty matrices, where the weights are the smoothing parameters. In the mixed model formulation the penalty matrix is the inverse of the covariance matrix for the random effects of a term, and the smoothing parameters (times a half) are variance parameters to be estimated. It's not possible to transform the problem to make the required random effects covariance matrix look like one of the standard pdMat classes: hence the need for the pdTens class. A notLog2 parameterization ensures that the parameters are positive.

These functions (pdTens, pdConstruct.pdTens, pdFactor.pdTens, pdMatrix.pdTens, coef.pdTens and summary.pdTens) would not normally be called directly.

Usage

pdTens(value = numeric(0), form = NULL,
       nam = NULL, data = sys.frame(sys.parent()))

Arguments

value | Initialization values for parameters. Not normally used.
form | A one sided formula specifying the random effects structure. The formula should have an attribute S which is a list of the penalty matrices the weighted sum of which gives the inverse of the covariance matrix for these random effects.
nam | a names argument, not normally used with this class.
data | data frame in which to evaluate formula.

Details

If using this class directly note that it is worthwhile scaling the S matrices to be of `moderate size', for example by dividing each matrix by its largest singular value: this avoids problems with lme defaults (\ref{smooth.construct.tensor.smooth.spec} does this automatically).

This appears to be the minimum set of functions required to implement a new pdMat class.

Note that while the pdFactor and pdMatrix functions return the inverse of the scaled random effect covariance matrix or its factor, the pdConstruct function is sometimes initialised with estimates of the scaled covariance matrix, and sometimes intialized with its inverse.

Value

A class pdTens object, or its coefficients or the matrix it represents or the factor of that matrix. pdFactor returns the factor as a vector (packed column-wise) (pdMatrix always returns a matrix).

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Pinheiro J.C. and Bates, D.M. (2000) Mixed effects Models in S and S-PLUS. Springer

The nlme source code.

http://www.maths.bath.ac.uk/~sw283/

See Also

te gamm

Examples

# see gamm


Automatically place a set of knots evenly through covariate values

Description

Given a univariate array of covariate values, places a set of knots for a regression spline evenly through the covariate values.

Usage

 place.knots(x,nk)

Arguments

x | array of covariate values (need not be sorted).
nk | integer indicating the required number of knots.

Details

Places knots evenly throughout a set of covariates. For example, if you had 11 covariate values and wanted 6 knots then a knot would be placed at the first (sorted) covariate value and every second (sorted) value thereafter. With less convenient numbers of data and knots the knots are placed within intervals between data in order to achieve even coverage, where even means having approximately the same number of data between each pair of knots.

Value

An array of knot locations.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

smooth.construct.cc.smooth.spec

Examples

x<-runif(30)
place.knots(x,7)
rm(x)


Default GAM plotting

Description

Takes a fitted gam object produced by gam() and plots the component smooth functions that make it up, on the scale of the linear predictor. Optionally produces term plots for parametric model components as well.

Usage

## S3 method for class 'gam':
plot(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scale=-1,
         n=100,n2=40,pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,
         ylab=NULL,main=NULL,ylim=NULL,xlim=NULL,too.far=0.1,
         all.terms=FALSE,shade=FALSE,shade.col="gray80",
         shift=0,trans=I,...)

Arguments

x | a fitted gam object as produced by gam().
residuals | If TRUE then partial residuals are added to plots of 1-D smooths. If FALSE then no residuals are added. If this is an array of the correct length then it is used as the array of residuals to be used for producing partial residuals. If TRUE then the residuals are the working residuals from the IRLS iteration weighted by the IRLS weights. Partial residuals for a smooth term are the residuals that would be obtained by dropping the term concerned from the model, while leaving all other estimates fixed (i.e. the estimates for the term plus the residuals).
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.
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. See also shade, below.
pages | (default 0) the number of pages over which to spread the output. For example, if 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.
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 select=2.
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. Ignored if ylim supplied.
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.
n2 | Square root of number of points used to grid estimates of 2-d functions for contouring.
pers | Set to TRUE if you want perspective plots for 2-d terms.
theta | One of the perspective plot angles.
phi | The other perspective plot angle.
jit | Set to TRUE if you want rug plots for 1-d terms to be jittered.
xlab | If supplied then this will be used as the x label for all plots.
ylab | If supplied then this will be used as the y label for all plots.
main | Used as title (or z axis label) for plots if supplied.
ylim | If supplied then this pair of numbers are used as the y limits for each plot.
xlim | If supplied then this pair of numbers are used as the x limits for each plot.
too.far | If greater than 0 then this is used to determine when a location is too far from data to be plotted when plotting 2-D smooths. This is useful since smooths tend to go wild away from data. The data are scaled into the unit square before deciding what to exclude, and too.far is a distance within the unit square.
all.terms | if set to TRUE then the partial effects of parametric model components are also plotted, via a call to termplot. Only terms of order 1 can be plotted in this way.
shade | Set to TRUE to produce shaded regions as confidence bands for smooths (not avaliable for parametric terms, which are plotted using termplot).
shade.col | define the color used for shading confidence bands.
shift | constant to add to each smooth (on the scale of the linear predictor) before plotting. Can be useful for some diagnostics, or with trans.
trans | function to apply to each smooth (after any shift), before plotting. shift and trans are occasionally useful as a means for getting plots on the response scale, when the model consists only of a single smooth.
... | other graphics parameters to pass on to plotting commands.

Details

Produces default plot showing the smooth components of a fitted GAM, and optionally parametric terms as well, when these can be handled by termplot.

For plots of 1-d smooths, the x axis of each plot is labelled with the covariate name, while the y axis is labelled s(cov,edf) where cov is the covariate name, and 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 s(var1,var2,edf), indicating the variables of which the term is a function, and the estimated degrees of freedom for the term. When se=TRUE, estimator variability is shown by overlaying contour plots at plus and minus 1 s.e. relative to the main estimate. If se is a positive number then contour plots are at plus or minus 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. Note that these plots can not be modified to the same extent as the other plot.

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

Fine control of plots for parametric terms can be obtained by calling termplot directly, taking care to use its terms argument.

Value

The function simply generates plots.

WARNING

Note that the behaviour of this function is not identical to plot.gam() in S-PLUS.

Plots of 2-D smooths with standard error contours shown can not easily be customized.

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

Author(s)

Simon N. Wood simon.wood@r-project.org

Henric Nilsson henric.nilsson@statisticon.se donated the code for the shade option.

The design is inspired by the S function of the same name described in Chambers and Hastie (1993) (but is not a clone).

References

Chambers and Hastie (1993) Statistical Models in S. Chapman & 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

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

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

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam, predict.gam, vis.gam

Examples

library(mgcv)
set.seed(0)
## fake some data...

f1 <- function(x) {exp(2 * x)}
f2 <- function(x) {
  0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
}
f3 <- function(x) {x*0}

n<-200
sig2<-4
x0 <- rep(1:4,50)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
e <- rnorm(n, 0, sqrt(sig2))
y <- 2*x0 + f1(x1) + f2(x2) + f3(x3) + e
x0 <- factor(x0)

## fit and plot...
b<-gam(y~x0+s(x1)+s(x2)+s(x3))
plot(b,pages=1,residuals=TRUE,all.terms=TRUE,shade=TRUE,shade.col=2)

## just parametric term alone...
termplot(b,terms="x0",se=TRUE)

## more use of color...
op <- par(mfrow=c(2,2),bg="blue")
x <- 0:1000/1000
for (i in 1:3) {
  plot(b,select=i,rug=FALSE,col="green",
    col.axis="white",col.lab="white",all.terms=TRUE)
  for (j in 1:2) axis(j,col="white",labels=FALSE)
  box(col="white")
  eval(parse(text=paste("fx <- f",i,"(x)",sep="")))
  fx <- fx-mean(fx)
  lines(x,fx,col=2) ## overlay `truth' in red
}
par(op)

## example with 2-d plots...
b1<-gam(y~x0+s(x1,x2)+s(x3))
op<-par(mfrow=c(2,2))
plot(b1,all.terms=TRUE)
par(op)


Prediction from fitted GAM model

Description

Takes a fitted gam object produced by gam() and produces predictions given a new set of values for the model covariates or the original values used for the model fit. Predictions can be accompanied by standard errors, based on the posterior distribution of the model coefficients. The routine can optionally return the matrix by which the model coefficients must be pre-multiplied in order to yield the values of the linear predictor at the supplied covariate values: this is useful for obtaining credible regions for quantities derived from the model, and for lookup table prediction outside R (see example code below).

Usage

## S3 method for class 'gam':
predict(object,newdata,type="link",se.fit=FALSE,terms=NULL,
        block.size=1000,newdata.guaranteed=FALSE,na.action=na.pass,...)

Arguments

object | a fitted gam object as produced by gam().
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. If newdata is provided then it should contain all the variables needed for prediction: a warning is generated if not.
type | When this has the value "link" (default) the linear predictor (possibly with associated standard errors) is returned. When 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 and any intercept. When type="response" predictions on the scale of the response are returned (possibly with approximate standard errors). When type="lpmatrix" then a matrix is returned which yields the values of the linear predictor (minus any offset) when postmultiplied by the parameter vector (in this case se.fit is ignored). The latter option is most useful for getting variance estimates for quantities derived from the model: for example integrated quantities, or derivatives of smooths. A linear predictor matrix can also be used to implement approximate prediction outside R (see example code, below).
se.fit | when this is TRUE (not default) standard error estimates are returned for each prediction.
terms | if type=="terms" then only results for the terms given in this array will be returned.
block.size | maximum number of predictions to process per call to underlying code: larger is quicker, but more memory intensive. Set to < 1 to use total number of predictions as this.
newdata.guaranteed | Set to TRUE to turn off all checking of newdata except for sanity of factor levels: this can speed things up for large prediction tasks, but newdata must be complete, with no NA values for predictors required in the model.
na.action | what to do about NA values in newdata. With the default na.pass, any row of newdata containing NA values for required predictors, gives rise to NA predictions (even if the term concerned has no NA predictors). na.exclude or na.omit result in the dropping of newdata rows, if they contain any NA values for required predictors. If newdata is missing then NA handling is determined from object$na.action.
... | other arguments.

Details

The standard errors produced by predict.gam are based on the Bayesian posterior covariance matrix of the parameters Vp in the fitted gam object.

To facilitate plotting with termplot, if object possesses an attribute "para.only" and type=="terms" then only parametric terms of order 1 are returned (i.e. those that termplot can handle).

Note that, in common with other prediction functions, any offset supplied to gam as an argument is always ignored when predicting, unlike offsets specified in the gam model formula.

See the examples for how to use the lpmatrix for obtaining credible regions for quantities derived from the model.

Value

If 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 se.fit is TRUE then a 2 item list is returned with items (both arrays) fit and 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 type is "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 or the intercept.
newdata can be a data frame, list or model.frame: if it's a model frame then all variables must be supplied.

WARNING

Note that the behaviour of this function is not identical to predict.gam() in Splus.

type=="terms" does not exactly match what predict.lm does for parametric model components.

Author(s)

Simon N. Wood simon.wood@r-project.org

The design is inspired by the S function of the same name described in Chambers and Hastie (1993) (but is not a clone).

References

Chambers and Hastie (1993) Statistical Models in S. Chapman & 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

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

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

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam, gamm, plot.gam

Examples

library(mgcv)
n<-200
sig <- 2
x0 <- runif(n, 0, 1)
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
y <- 2 * sin(pi * x0)
y <- y + exp(2 * x1)
y <- y + 0.2 * x2^11 * (10 * (1 - x2))^6 + 10 * (10 * x2)^3 * (1 - x2)^10
y <- y + x3
e <- rnorm(n, 0, sig)
y <- y + e
b<-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3))
rm(y,x0,x1,x2,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)

## now get variance of sum of predictions using lpmatrix

Xp <- predict(b,newd,type="lpmatrix")

## Xp %*% coef(b) yields vector of predictions

a <- rep(1,31)
Xs <- t(a) %*% Xp ## Xs %*% coef(b) gives sum of predictions
var.sum <- Xs %*% b$Vp %*% t(Xs)

## Now get the variance of non-linear function of predictions
## by simulation from posterior distribution of the params

library(MASS)
br<-mvrnorm(1000,coef(b),b$Vp) ## 1000 replicate param. vectors
res <- rep(0,1000)
for (i in 1:1000)
{ pr <- Xp %*% br[i,] ## replicate predictions
  res[i] <- sum(log(abs(pr))) ## example non-linear function
}
mean(res);var(res)

## loop is replace-able by following ....

res <- colSums(log(abs(Xp %*% t(br))))

## The following shows how to use use an "lpmatrix" as a lookup
## table for approximate prediction. The idea is to create
## approximate prediction matrix rows by appropriate linear
## interpolation of an existing prediction matrix. The additivity
## of a GAM makes this possible.
## There is no reason to ever do this in R, but the following
## code provides a useful template for predicting from a fitted
## gam *outside* R: all that is needed is the coefficient vector
## and the prediction matrix. Use larger `Xp'/ smaller `dx' and/or
## higher order interpolation for higher accuracy.

xn <- c(.341,.122,.476,.981) ## want prediction at these values
x0 <- 1         ## intercept column
dx <- 1/30      ## covariate spacing in `newd'
for (j in 0:2) { ## loop through smooth terms
  cols <- 1+j*9 +1:9      ## relevant cols of Xp
  i <- floor(xn[j+1]*30)  ## find relevant rows of Xp
  w1 <- (xn[j+1]-i*dx)/dx ## interpolation weights
  ## find approx. predict matrix row portion, by interpolation
  x0 <- c(x0,Xp[i+2,cols]*w1 + Xp[i+1,cols]*(1-w1))
}
dim(x0)<-c(1,28)
fv <- x0%*%coef(b) + xn[4];fv    ## evaluate and add offset
se <- sqrt(x0%*%b$Vp%*%t(x0));se ## get standard error
## compare to normal prediction
predict(b,newdata=data.frame(x0=xn[1],x1=xn[2],
        x2=xn[3],x3=xn[4]),se=TRUE)


Prediction methods for smooth terms in a GAM

Description

Takes smooth objects produced by smooth.construct methods and obtains the matrix mapping the parameters associated with such a smooth to the predicted values of the smooth at a set of new covariate values.

In practice this method is often called via the wrapper function PredictMat.

Usage

Predict.matrix(object,data)

Arguments

object | is a smooth object produced by a smooth.construct method function. The object contains all the information required to specify the basis for a term of its class, and this information is used by the appropriate Predict.matrix function to produce a prediction matrix for new covariate values. Further details are given in smooth.construct.
data | A data frame containing the values of the (named) covariates at which the smooth term is to be evaluated.

Details

Smooth terms in a GAM formula are turned into smooth specification objects of class xx.smooth.spec during processing of the formula. Each of these objects is converted to a smooth object using an appropriate smooth.construct function. The Predict.matrix functions are used to obtain the matrix that will map the parameters associated with a smooth term to the predicted values for the term at new covariate values.

Note that new smooth classes can be added by writing a new smooth.construct method function and a corresponding Predict.matrix method function: see the example code provided for smooth.construct for details.

Value

A matrix which will map the parameters associated with the smooth to the vector of values of the smooth evaluated at the covariate values given in object.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

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

http://www.maths.bath.ac.uk/~sw283/

See Also

gam,gamm, smooth.construct, PredictMat

Examples

# See smooth.construct examples


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 gam, or use 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.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Wood, S.N. (2006) Generalized Additive Models: An Introduction with R. CRC/ Chapmand and Hall, Boca Raton, Florida.

http://www.maths.bath.ac.uk/~sw283/

See Also

gam,mgcv


Generalized Additive Model residuals

Description

Returns residuals for a fitted gam model object. Pearson, deviance, working and response residuals are available.

Usage

## S3 method for class 'gam':
residuals(object, type = c("deviance", "pearson","scaled.pearson", "working", "response"),...)

Arguments

object | a gam fitted model object.
type | the type of residuals wanted.
... | 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): ((y-m)/V(m)^0.5, where y is data m is model fitted value and 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 gam objects because of a bug in the calculation of Pearson residuals in some earlier versions of residual.glm.

Value

An array of residuals.

Author(s)

Simon N. Wood simon.wood@r-project.org

See Also

gam


Defining smooths in GAM formulae

Description

Function used in definition of smooth terms within 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)

Arguments

... | a list of variables that are the covariates that this smooth is a function of.
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. k should not be less than the dimension of the null space of the penalty for the term (see null.space.dimension), but will be reset if it is. See choose.k for further information.
fx | indicates whether the term is a fixed d.f. regression spline (TRUE) or a penalized regression spline (FALSE).
bs | this can be "cr" for a cubic regression spline, "cs" for a cubic regression spline with shrinkage, "cc" for a cyclic (periodic) spline, "tp" for a thin plate regression spline, "ts" for a thin plate regression spline with shrinkage or a user defined character string for other user defined smooth classes. Of the built in alternatives, only thin plate regression splines can be used for multidimensional smooths. Note that the "cr" and "cc" bases are faster to set up than the "tp" basis, particularly on large data sets (although the knots argument to gam can be used to get round this).
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.
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 by variable would usually be the dummy variable coding one level of the factor. See the examples below. This is the means by which `variable coefficient models' can be specified.

Details

The function does not evaluate the variable arguments. To use this function to specify use of your own smooths, note the relationships between the inputs and the output object and see the example in smooth.construct.

Value

A class xx.smooth.spec object, where xx is a basis identifying code given by the bs argument of s. These smooth.spec objects define smooths and are turned into bases and penalties by smooth.construct method functions.
The returned object contains the following items:

term | An array of text strings giving the names of the covariates that the term is a function of.
bs.dim | The dimension of the basis used to represent the smooth.
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
dim | The dimension of the smoother - i.e. the number of covariates that it is a function of.
p.order | The order of the t.p.r.s. penalty, or 0 for auto-selection of the penalty order.
by | is the name of any by variable as text ("NA" for none).
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 s() with all the options spelled out explicitly.
label | A suitable text label for this smooth term.

Author(s)

Simon N. Wood simon.wood@r-project.org

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

http://www.maths.bath.ac.uk/~sw283/

See Also

te, gam, gamm

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
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)


Constructor functions for smooth terms in a GAM

Description

Smooth terms in a GAM formula are turned into smooth specification objects of class xx.smooth.spec during processing of the formula. Each of these objects is converted to a smooth object using an appropriate smooth.construct function. New smooth classes can be added by writing a new smooth.construct method function and a corresponding Predict.matrix method function (see example code below).

In practice, smooth.construct is usually called via the wrapper function smoothCon.

Usage

smooth.construct(object,data,knots)

Arguments

object is a smooth specification object, generated by an s or te term in a GAM formula. Objects generated by s terms have class xx.smooth.spec where xx is given by the bs argument of s (this convention allows the user to add their own smoothers). If object is not class tensor.smooth.spec it will have the following elements:
    term
    The names of the covariates for this smooth, in an array.
    bs.dim
    Argument k of the s term generating the object. This is the dimension of the basis used to represent the term (or, arguably, 1 greater than the basis dimension for cc terms).
    fixed
    TRUE if the term is to be unpenalized, otherwise FALSE.
    dim
    the number covariates of which this smooth is a function.
    p.order
    the order of the smoothness penalty or 0 for autoselection of this. This is argument m of the s term that generated object.
    by
    the name of any by variable to multiply this term as supplied as an argument to s. "NA" if there is no such term.
    full.call
    The full version of the s term, with all defaults expanded explicitly.
    label
    A suitable label for use with this term.
    null.space.dim
    The dimension of the null space of the wiggliness penalty.
If object is of class tensor.smooth.spec then it was generated by a te term in the GAM formula, and specifies a smooth of several variables with a basis generated as a tensor product of lower dimensional bases. In this case the object will be different and will have the following elements:
    margin
    is a list of smooth specification objects of the type listed above, defining the bases which have their tensor product formed in order to construct this term.
    term
    is the array of names of the covariates that are arguments of the smooth.
    by
    is the name of any by variable, or "NA".
    fx
    is an array, the elements of which indicate whether (TRUE) any of the margins in the tensor product should be unpenalized.
    full.call
    The full version of the s term, with all defaults expanded explicitly.
    label
    A suitable label for use with this term.
    dim
    is the number of covariates of which this smooth is a function.
    mp
    TRUE if multiple penalties are to be used.
    np
    TRUE if 1-D marginal smooths are to be re-parameterized in terms of function values.
data | a data frame in which the covariates and any by variable can be found.
knots | an optional data frame specifying knot locations for each covariate. If it is null then the knot locations are generated automatically.

Details

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

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

tprs.smooth objects require several items to be stored in order to define the basis. These are:

    shift
    A record of the shift applied to each covariate in order to center it around zero and avoid any co-linearity problems that might otehrwise occur in the penalty null space basis of the term.
    Xu
    A matrix of the unique covariate combinations for this smooth (the basis is constructed by first stripping out duplicate locations).
    UZ
    The matrix mapping the t.p.r.s. parameters back to the parameters of a full thin plate spline.
    null.space.dimension
    The dimension of the space of functions that have zero wiggliness according to the wiggliness penalty for this term.

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

The classes cs.smooth and ts.smooth are modifications of the classes cr.smooth and tprs.smooth respsectively: they differ only in having an extra shrinkage component added to their penalties, so that automatic model selection can effectively remove such terms from a model altogether.

Value

The input argument object, assigned a new class to indicate what type of smooth it is and with at least the following items added:

X | The model matrix from this term.
C | The matrix defining any constraints on the term - usually a one row matrix giving the column sums of the model matrix, which defines the constraint that each term should sum to zero over the covariate values.
S | A list of positive semi-definite penalty matrices that apply to this term. The list will be empty if the term is to be left un-penalized.
rank | an array giving the ranks of the penalties.
df | the degrees of freedom associated with this term (at least when unpenalized).

Usually the returned object will also include extra information required to define the basis, and used by Predict.matrix methods to make predictions using the basis. See the Details section for the information included for the built in smooth classes.
tensor.smooth returned objects will additionally have each element of the margin list updated in the same way. tensor.smooths also have a list, XP, containing re-parameterization matrices for any 1-D marginal terms re-parameterized in terms of function values. This list will have NULL entries for marginal smooths that are not re-parameterized, and is only long enough to reach the last re-parameterized marginal in the list.

WARNING

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

Author(s)

Simon N. Wood simon.wood@r-project.org

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

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

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

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

http://www.maths.bath.ac.uk/~sw283/

See Also

get.var, gamm, gam, Predict.matrix, smoothCon, PredictMat

Examples

# adding "p-spline" classes and methods

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

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

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

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


Prediction/Construction wrapper functions for GAM smooth terms

Description

Wrapper functions for construction of and prediction from smooth terms in a GAM. The purpose of the wrappers is to allow user-transparant re-parameterization of smooth terms, in order to allow identifiability constraints to be absorbed into the parameterization of each term, if required.

Usage

smoothCon(object,data,knots,absorb.cons=FALSE,scale.penalty=TRUE)
PredictMat(object,data)

Arguments

object | is a smooth specification object or a smooth object.
data | A data frame containing the values of the (named) covariates at which the smooth term is to be evaluated.
knots | An optional data frame supplying any knot locations to be supplied for basis construction.
absorb.cons | Set to TRUE in order to have identifiability constraints absorbed into the basis.
scale.penalty | should the penalty coefficient matrix be scaled to have approximately the same `size' as the inner product of the terms model matrix with itself? This can improve the performance of gamm fitting.

Details

These wrapper functions exist to allow smooths specified using smooth.construct and Predict.matrix method functions to be re-parameterized so that identifiability constraints are no longer required in fitting. This is done in a user transparent manner, but is typically of no importance in use of GAMs.

The parameterization used by gam can be controlled via gam.control.

Value

From smoothCon a smooth object returned by the appropriate smooth.construct method function. If constraints are to be absorbed then the object will have attributes "qrc" and "nCons", the qr decomposition of the constraint matrix (returned by qr) and the number of constraints, respectively: these are used in the re-parameterization.
For predictMat a matrix which will map the parameters associated with the smooth to the vector of values of the smooth evaluated at the covariate values given in object.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

gam.control, smooth.construct, Predict.matrix


Alternatives to step.gam

Description

There is no step.gam in package mgcv. The mgcv default for model selection is to use MSE/KL-distance criteria such as GCV or UBRE/AIC. Since the smoothness estimation part of model selection is done in this way it is logically most consistent to perform model selection on the basis of such criteria: i.e. to decide which terms to include or omit by looking at changes in GCV/UBRE score.

To facilitate fully automatic model selection the package includes 2 classes of smoothers ("cs" and "ts": see s) which can be penalized to zero for sufficiently high smoothing parameter estimates: use of such smooths provides an effective alternative to step-wise model selection. The example below shows an example of the application of this approach, where selection is a fully integrated part of model estimation.

Author(s)

Simon N. Wood simon.wood@r-project.org

Examples

## an example of GCV based model selection as
## an alternative to stepwise selection
library(mgcv)
set.seed(0)
n<-400;sig<-2
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
x4 <- runif(n, 0, 1);x5 <- runif(n, 0, 1)
f <- 2 * sin(pi * x0)
f <- f + exp(2 * x1) - 3.75887
f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
e <- rnorm(n, 0, sig)
y <- f + e
## Note the increased gamma parameter below to favour
## slightly smoother models...
b<-gam(y~s(x0,bs="ts")+s(x1,bs="ts")+s(x2,bs="ts")+
   s(x3,bs="ts")+s(x4,bs="ts")+s(x5,bs="ts"),gamma=1.4)
summary(b)
plot(b,pages=1)


Summary for a GAM fit

Description

Takes a fitted gam object produced by gam() and produces various useful summaries from it.

Usage

## S3 method for class 'gam':
summary(object, dispersion=NULL, freq=TRUE, ...)

## S3 method for class 'summary.gam':
print(x,digits = max(3, getOption("digits") - 3),
                  signif.stars = getOption("show.signif.stars"),...)

Arguments

object | a fitted gam object as produced by gam().
x | a summary.gam object produced by summary.gam().
dispersion | A known dispersion parameter. NULL to use estimate or default (e.g. 1 for Poisson).
freq | By default p-values for individual terms are calculated using the frequentist estimated covariance matrix of the parameter estimators. If this is set to FALSE then the Bayesian covariance matrix of the parameters is used instead. See details.
digits | controls number of digits printed in output.
signif.stars | Should significance stars be printed alongside output.
... | other arguments.

Details

Model degrees of freedom are taken as the trace of the influence (or hat) matrix A for the model fit. Residual degrees of freedom are taken as number of data minus model degrees of freedom. Let 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 X be the design matrix of the model. Then tr(XP_i) is the edf for the ith term. Clearly this definition causes the edf's to add up properly!

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

If freq=FALSE then the Bayesian parameter covariance matrix, object$Vp, is used to calculate test statistics for terms, and the degrees of freedom for reference distributions is taken as the estimated degrees of freedom for the term concerned. This is not easy to justify theoretically, and the resulting `Bayesian p-values' are difficult to interpret and often have much worse frequentist performance than the default p-values.

Value

summary.gam produces a list of summary information for a fitted gam object.

p.coeff | is an array of estimates of the strictly parametric model coefficients.
p.t | is an array of the p.coeff's divided by their standard errors.
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 if the dispersion parameter has been estimated, and the standard normal if not.
m | The number of smooth terms in the model.
chi.sq | An array of test statistics for assessing the significance of model smooth terms. If p_i is the parameter vector for the ith smooth term, and this term has estimated covariance matrix V_i then the statistic is p_i'V_i^{k-}p_i, where V_i^{k-} is the rank k pseudo-inverse of V_i, and k is estimated rank of V_i.
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 of known dispersion parameter, they are obtained by comparing the chi.sq statistic given above to the chi-squared distribution with k degrees of freedom, where k is the estimated rank of V_i. If the dispersion parameter is unknown (in which case it will have been estimated) the statistic is compared to an F distribution with k upper d.f. and lower d.f. given by the residual degrees of freedom for the model . Typically the p-values will be somewhat too low, because they are conditional on the smoothing parameters, which are usually uncertain, but note that the statistic can also have low power if the rank, k, is too high relative to the EDF of the term!
se | array of standard error estimates for all parameter estimates.
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.
dev.expl | The proportion of the null deviance explained by the model.
edf | array of estimated degrees of freedom for the model terms.
residual.df | estimated residual degrees of freedom.
n | number of data.
gcv | minimized GCV score for the model, if GCV used.
ubre | minimized UBRE score for the model, if UBRE used.
scale | estimated (or given) scale parameter.
family | the family used.
formula | the original GAM formula.
dispersion | the scale parameter.
pTerms.df | the degrees of freedom associated with each parameteric term (excluding the constant).
pTerms.chi.sq | a Wald statistic for testing the null hypothesis that the each parametric term is zero.
pTerms.pv | p-values associated with the tests that each term is zero. For penalized fits these are approximate, being conditional on the smoothing parameters. The reference distribution is an appropriate chi-squared when the scale parameter is known, and is based on an F when it is not.
cov.unscaled | The estimated covariance matrix of the parameters (or estimators if freq=TRUE), divided by scale parameter.
cov.scaled | The estimated covariance matrix of the parameters (estimators if freq=TRUE).
p.table | significance table for parameters
s.table | significance table for smooths
p.Terms | significance table for parametric model terms

WARNING

The supplied p-values will often be underestimates if smoothing parameters have been estimated as part of model fitting.

Author(s)

Simon N. Wood simon.wood@r-project.org with substantial improvements by Henric Nilsson.

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

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

Wood, S.N. (2004b) On confidence intervals for GAMs based on penalized regression splines. Technical Report 04-12 Department of Statistics, University of Glasgow.

http://www.maths.bath.ac.uk/~sw283/

See Also

gam, predict.gam, gam.check, anova.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
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(summary(b)$edf)+1 ## get edf per smooth
b.d<-pmax(b.d,3) # can't have basis dimension less than 3!
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))
abline(0,1,col=2)


Define tensor product smooths in GAM formulae

Description

Function used in definition of tensor product smooth terms within gam model formulae. The function does not evaluate a smooth - it exists purely to help set up a model using tensor product based smooths.

Usage

te(..., k=NA,bs="cr",m=0,d=NA,by=NA,fx=FALSE,mp=TRUE,np=TRUE)

Arguments

... | a list of variables that are the covariates that this smooth is a function of.
k | the dimension(s) of the bases used to represent the smooth term. If not supplied then set to 5^d. If supplied as a single number then this basis dimension is used for each basis. If supplied as an array then the elements are the dimensions of the component (marginal) bases of the tensor product. See choose.k for further information.
bs | array (or single character string) specifying the type for each marginal basis. "cr" for cubic regression spline; "cs" for cubic regression spline with shrinkage; "cc" for periodic/cyclic cubic regression spline; "tp" for thin plate regression spline; "ts" for t.p.r.s. with extra shrinkage. User defined bases can also be used here (see smooth.construct for an example). If only one identifier is given then this is used for all bases.
m | The order of the penalty for each t.p.r.s. term (e.g. 2 for normal cubic spline penalty with 2nd derivatives). If a single number is given then it is used for all terms. 0 autoinitializes. m is ignored for the "cr" and "cc" bases.
d | array of marginal basis dimensions. For example if you want a smooth for 3 covariates made up of a tensor product of a 2 dimensional t.p.r.s. basis and a 1-dimensional basis, then set d=c(2,1).
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 by variable would usually be the dummy variable coding one level of the factor. See the examples below.
fx | indicates whether the term is a fixed d.f. regression spline (TRUE) or a penalized regression spline (FALSE).
mp | TRUE to use multiple penalties for the smooth. FALSE to use only a single penalty: single penalties are not recommended - they tend to allow only rather wiggly models.
np | TRUE to use the `normal parameterization' for a tensor product smooth. This parameterization represents any 1-d marginal smooths using a parameterization where the parameters are function values at `knots' spread evenly through the data. The parameterization makes the penalties easily interpretable, however it can reduce numerical stability.

Details

Smooths of several covariates can be constructed from tensor products of the bases used to represent smooths of one (or sometimes more) of the covariates. To do this `marginal' bases are produced with associated model matrices and penalty matrices, and these are then combined in the manner described in tensor.prod.model.matrix and tensor.prod.penalties, to produce a single model matrix for the smooth, but multiple penalties (one for each marginal basis). The basis dimension of the whole smooth is the product of the basis dimensions of the marginal smooths.

An option for operating with a single penalty (The Kronecker product of the marginal penalties) is provided, but it is rarely of practical use: the penalty is typically so rank deficient that even the smoothest resulting model will have rather high estimated degrees of freedom.

Tensor product smooths are especially useful for representing functions of covariates measured in different units, although they are typically not quite as nicely behaved as t.p.r.s. smooths for well scaled covariates.

Note also that GAMs constructed from lower rank tensor product smooths are nested within GAMs constructed from higher rank tensor product smooths if the same marginal bases are used in both cases (the marginal smooths themselves are just special cases of tensor product smooths.)

The `normal parameterization' (np=TRUE) re-parameterizes the marginal smooths of a tensor product smooth so that the parameters are function values at a set of points spread evenly through the range of values of the covariate of the smooth. This means that the penalty of the tensor product associated with any particular covariate direction can be interpreted as the penalty of the appropriate marginal smooth applied in that direction and averaged over the smooth. Currently this is only done for marginals of a single variable. This parameterization can reduce numerical stability for when used with marginal smooths other than "cc", "cr" and "cs": if this causes problems, set np=FALSE.

The function does not evaluate the variable arguments.

Value

A class tensor.smooth.spec object defining a tensor product smooth to be turned into a basis and penalties by the smooth.construct.tensor.smooth.spec function.
The returned object contains the following items:

margin | A list of smooth.spec objects of the type returned by s, defining the basis from which the tensor product smooth is constructed.
term | An array of text strings giving the names of the covariates that the term is a function of.
by | is the name of any by variable as text ("NA" for none).
fx | logical array with element for each penalty of the term (tensor product smooths have multiple penalties). TRUE if the penalty is to be ignored, FALSE, otherwise.
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 s() with all the options spelled out explicitly.
label | A suitable text label for this smooth term.
dim | The dimension of the smoother - i.e. the number of covariates that it is a function of.
mp | TRUE is multiple penalties are to be used (default).
np | TRUE to re-parameterize 1-D marginal smooths in terms of function values (defualt).

Author(s)

Simon N. Wood simon.wood@r-project.org

References

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for Generalized Additive Mixed Models. Biometrics

http://www.maths.bath.ac.uk/~sw283/

See Also

s,gam,gamm

Examples

# following shows how tensor pruduct deals nicely with
# badly scaled covariates (range of x 5% of range of z )
test1<-function(x,z,sx=0.3,sz=0.4)
{ x<-x*20
  (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)/20;z<-runif(n);
xs<-seq(0,1,length=30)/20;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)
f <- test1(x,z)
y <- f + rnorm(n)*0.2
b1<-gam(y~s(x,z))
persp(xs,zs,truth);title("truth")
vis.gam(b1);title("t.p.r.s")
b2<-gam(y~te(x,z))
vis.gam(b2);title("tensor product")
b3<-gam(y~te(x,z,bs=c("tp","tp")))
vis.gam(b3);title("tensor product")
par(old.par)

test2<-function(u,v,w,sv=0.3,sw=0.4)
{ ((pi**sv*sw)*(1.2*exp(-(v-0.2)^2/sv^2-(w-0.3)^2/sw^2)+
  0.8*exp(-(v-0.7)^2/sv^2-(w-0.8)^2/sw^2)))*(u-0.5)^2*20
}
n <- 500
v <- runif(n);w<-runif(n);u<-runif(n)
f <- test2(u,v,w)
y <- f + rnorm(n)*0.2
# tensor product of a 2-d thin plate regression spline and 1-d cr spline
b <- gam(y~te(v,w,u,k=c(30,5),d=c(2,1),bs=c("tp","cr")))
op <- par(mfrow=c(2,2))
vis.gam(b,cond=list(u=0),color="heat",zlim=c(-0.2,3.5))
vis.gam(b,cond=list(u=.33),color="heat",zlim=c(-0.2,3.5))
vis.gam(b,cond=list(u=.67),color="heat",zlim=c(-0.2,3.5))
vis.gam(b,cond=list(u=1),color="heat",zlim=c(-0.2,3.5))
par(op)


Utility functions for constructing tensor product smooths

Description

Produce model matrices or penalty matrices for a tensor product smooth from the model matrices or penalty matrices for the marginal bases of the smooth.

Usage

tensor.prod.model.matrix(X)
tensor.prod.penalties(S)

Arguments

X | a list of model matrices for the marginal bases of a smooth
S | a list of penalties for the marginal bases of a smooth.

Details

If X[[1]], X[[2]] ... X[[m]] are the model matrices of the marginal bases of a tensor product smooth then the ith row of the model matrix for the whole tensor product smooth is given by X[[1]][i,]%x%X[[2]][i,]%x% ... X[[m]][i,], where %x% is the Kronecker product. Of course the routine operates column-wise, not row-wise!

If S[[1]], S[[2]] ... S[[m]] are the penalty matrices for the marginal bases, and I[[1]], I[[2]] ... I[[m]] are corresponding identity matrices, each of the same dimension as its corresponding penalty, then the tensor product smooth has m associate penalties of the form:

S[[1]]%x%I[[2]]%x% ... I[[m]],

I[[1]]%x%S[[2]]%x% ... I[[m]]

...

I[[1]]%x%I[[2]]%x% ... S[[m]].

Of course it's important that the model matrices and penalty matrices are presented in the same order when constructing tensor product smooths.

Value

Either a single model matrix for a tensor product smooth, or a list of penalty terms for a tensor product smooth.

Author(s)

Simon N. Wood simon.wood@r-project.org

References

http://www.maths.bath.ac.uk/~sw283/

See Also

te, smooth.construct.tensor.smooth.spec

Examples

X <- list(matrix(1:4,2,2),matrix(5:10,2,3))
tensor.prod.model.matrix(X)

S<-list(matrix(c(2,1,1,2),2,2),matrix(c(2,1,0,1,2,1,0,1,2),3,3))
tensor.prod.penalties(S)


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)

Arguments

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 x (in arbitrary order).

Author(s)

Simon N. Wood simon.wood@r-project.org

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)


Extract parameter (estimator) covariance matrix from GAM fit

Description

Extracts the Bayesian posterior covariance matrix of the parameters or frequentist covariance matrix of the parameter estimators from a fitted gam object.

Usage

## S3 method for class 'gam':
vcov(object, freq = TRUE, dispersion = NULL, ...)

Arguments

object | fitted model object of class gam as produced by gam().
freq | TRUE to return the frequentist covariance matrix of the parameter estimators, FALSE to return the Bayesian posterior covariance matrix of the parameters.
dispersion | a value for the dispersion parameter: not normally used.
... | other arguments, currently ignored.

Details

Basically, just extracts object$Ve or object$Vp from a gamObject.

Value

A matrix corresponding to the estimated frequentist covariance matrix of the model parameter estimators/coefficients, or the estimated posterior covariance matrix of the parameters, depending on the argument freq.

Author(s)

Henric Nilsson. Maintained by Simon N. Wood simon.wood@r-project.org

See Also

gam

Examples

n <- 100
x <- runif(n)
y <- sin(x*2*pi) + rnorm(n)*.2
mod <- gam(y~s(x,bs="cc",k=10),knots=list(x=seq(0,1,length=10)))
diag(vcov(mod))


Visualization of GAM objects

Description

Produces perspective or contour plot views of gam model predictions, fixing all but the values in view to the values supplied in cond.

Usage

vis.gam(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,
        color="heat",contour.col=NULL,se=-1,type="link",
        plot.type="persp",zlim=NULL,nCol=50,...)

Arguments

x | a gam object, produced by gam()
view | an array containing the names of the two main effect terms to be displayed on the x and y dimensions of the plot. If omitted the first two suitable terms will be used. Names must be names from names(x$model). i.e. terms are used, rather than variables.
cond | a named list of the values to use for the other predictor terms (not in view). Terms omitted from this list will have their values set to their mean for continuous variables, or first level for factors. Names must correspond to names(x$model).
n.grid | The number of grid nodes in each direction used for calculating the plotted surface.
too.far | plot grid nodes that are too far from the points defined by the variables given in view can be excluded from the plot. too.far determines what is too far. The grid is scaled into the unit square along with the view variables and then grid nodes more than too.far from the predictor variables are excluded.
col | The colours for the facets of the plot. If this is NA then if se>0 the facets are transparent, otherwise the colour scheme specified in color is used. If col is not NA then it is used as the facet colour.
color | the colour scheme to use for plots when se<=0. One of "topo", "heat", "cm", "terrain", "gray" or "bw". Schemes "gray" and "bw" also modify the colors used when se>0.
contour.col | sets the colour of contours when using plot.type="contour". Default scheme used if NULL.
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 se standard errors, one at the predicted values and one at the predicted values plus se standard errors.
type | "link" to plot on linear predictor scale and "response" to plot on the response scale.
plot.type | one of "contour" or "persp".
zlim | a two item array giving the lower and upper limits for the z-axis scale. NULL to choose automatically.
nCol | The number of colors to use in color schemes.
... | other options to pass on to persp, image or contour. In particular ticktype="detailed" will add proper axes labelling to the plots.

Details

The x and y limits are determined by the ranges of the terms named in view. If 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 +/- se standard errors. Parts of the x-y plane too far from data can be excluded by setting too.far

All options to the underlying graphics functions can be reset by passing them as extra arguments ...: such supplied values will always over-ride the defualt values used by vis.gam.

Value

Simply produces a plot.

Author(s)

Simon Wood simon.wood@r-project.org

Based on an original idea and design by Mike Lonergan.

See Also

persp and 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=.75),theta=210,phi=40,
        too.far=.07)
# contour examples....
vis.gam(g, view=c("x1","x2"),plot.type="contour",color="heat")
vis.gam(g, view=c("x1","x2"),plot.type="contour",color="terrain")
vis.gam(g, view=c("x1","x2"),plot.type="contour",color="topo")
vis.gam(g, view=c("x1","x2"),plot.type="contour",color="cm")

# ..... 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