Movatterモバイル変換

Mixed GAM Computation Vehicle with GCV/AIC/REML/NCV smoothness estimation and GAMMs by REML/PQL

Description

mgcv provides functions for generalized additive modelling (gam andbam) andgeneralized additive mixed modelling (gamm, andrandom.effects), including location scale and shape extensions. Shape constrained models can be estimated usingscasm. The term GAM is taken to include any model dependent on unknown smooth functions of predictors and estimated by quadratically penalized (possibly quasi-) likelihood maximization. Available distributions are covered infamily.mgcv and available smooths insmooth.terms.

Particular features of the package are facilities for automatic smoothness selection (Wood, 2004, 2011), and the provision of a variety of smooths of more than one variable. User defined smooths can be added. A Bayesian approach to confidence/credible interval calculation isprovided. Linear functionals of smooths, penalization of parametric model terms and linkage of smoothing parameters are all supported. Lower level routines for generalized ridge regression and penalized linearly constrained least squares are also available. In addition to the main modelling functions,jagam provided facilities to ease the set up of models for use with JAGS, whileginla provides marginal inference via a version of Integrated Nested Laplace Approximation.

Details

mgcv provides generalized additive modelling functionsgam,predict.gam andplot.gam, which are very similarin use to the S functions of the same name designed by Trevor Hastie (with some extensions). However the underlying representation and estimation of the models is based on apenalized regression spline approach, with automatic smoothness selection. Anumber of other functions such assummary.gam andanova.gam are also provided, for extracting information from a fittedgamObject.

Use ofgam is much like use ofglm, except thatwithin agam model formula, isotropic smooths of any number of predictors can be specified usings terms, while scale invariant smooths of any number ofpredictors can be specified usingte,ti ort2 terms.smooth.terms provides an overview of the built in smooth classes, andrandom.effects should be refered to for an overview of random effects terms (see alsomrf for Markov random fields). Estimation is bypenalized likelihood or quasi-likelihood maximization, with smoothnessselection by GCV, GACV, gAIC/UBRE,NCV or (RE)ML. Seegam,gam.models,linear.functional.terms andgam.selection for some discussion of model specification andselection. For detailed control of fitting seegam.convergence,gam argumentsmethod andoptimizer andgam.control. For checking andvisualization seegam.check,choose.k,vis.gam andplot.gam.While a number of types of smoother are built into the package, it is alsoextendable with user defined smooths, seesmooth.construct, for example.

A Bayesian approach to smooth modelling is used to derive standard errors onpredictions, and hence credible intervals (see Marra and Wood, 2012). The Bayesian covariance matrix forthe model coefficients is returned inVp of thegamObject. Seepredict.gam for examples of howthis can be used to obtain credible regions for any quantity derived from thefitted model, either directly, or by direct simulation from the posteriordistribution of the model coefficients. Approximate p-values can also be obtained for testing individual smooth terms for equality to the zero function, using similar ideas (see Wood, 2013a,b). Frequentistapproximations can be used for hypothesis testing based model comparison. Seeanova.gam andsummary.gam for more on hypothesis testing.

For large datasets (that is large n) seebam which is a version ofgam with a much reduced memory footprint.bam(...,discrete=TRUE) offers the very efficient methods of Wood et al. (2017) and Li and Wood (2020).

The package also provides a generalized additive mixed modelling function,gamm, based on a PQL approach andlme from thenlme library (for anlme4 based version, see packagegamm4).gamm is particularly usefulfor modelling correlated data (i.e. where a simple independence model for theresidual variation is inappropriate). Seerandom.effects for including random effects in models estimated bygam,bam orscasm.

In addition, low level routinemagiccan fit models to data with a known correlation structure.

Some underlying GAM fitting methods are available as low level fittingfunctions: seemagic. But there is little functionality that can not be more conventiently accessed viagam . Penalized weighted least squares with linear equality and inequality constraints is provided bypcls.

For a complete list of functions typelibrary(help=mgcv). See alsomgcv.FAQ.

Author(s)

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

with contributions and/or help from Natalya Pya, Thomas Kneib, Kurt Hornik, Mike Lonergan, Henric Nilsson,Fabian Scheipl and Brian Ripley.

Polish translation - Lukasz Daniel; German translation - Chris Leick, Detlef Steuer; French Translation - Philippe Grosjean

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

Part funded by EPSRC: EP/K005251/1

References

https://webhomes.maths.ed.ac.uk/~swood34/

These provide details for the underlying mgcv methods, and fuller references to the large literature on which the methods are based.

Wood, S.N. (2025) Generalized Additive Models. Annual Review ofStatistics and Its Applications 12:497-526doi:10.1146/annurev-statistics-112723-034249

Wood, S. N. (2020) Inference and computation with generalizedadditive models and their extensions. Test 29(2): 307-339.doi:10.1007/s11749-020-00711-5

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models (with discussion).Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

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

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals for Generalized AdditiveModel Components. Scandinavian Journal of Statistics, 39(1), 53-74.

Wood, S.N. (2013a) A simple test for random effects in regression models. Biometrika 100:1005-1010doi:10.1093/biomet/ast038

Wood, S.N. (2013b) On p-values for smooth components of an extended generalized additive model. Biometrika 100:221-228doi:10.1093/biomet/ass048

Wood, S.N. (2017)Generalized Additive Models: an introduction with R (2nd edition),CRCdoi:10.1201/9781315370279

Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association. 112(519):1199-1210doi:10.1080/01621459.2016.1195744

Li, Z & S.N. Wood (2020) Faster model matrix crossproducts for large generalized linear models with discretized covariates. Statistics and Computing. 30:19-25doi:10.1007/s11222-019-09864-2

Development of mgcv version 1.8 was part funded by EPSRC grants EP/K005251/1 and EP/I000917/1.

Examples

## see examples for gam, bam and gamm

Level 5 fractional factorial designs

Description

Computes level 5 fractional factorial designs for up to 120 factorsusing the agorithm of Sanchez and Sanchez (2005), and optionally central composite designs.

Usage

FFdes(size=5,ccd=FALSE)

Arguments

Details

Basically a translation of the code provided in the appendix of Sanchez and Sanchez (2005).

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Sanchez, S. M. & Sanchez, P. J. (2005) Very large fractional factorial and central composite designs.ACM Transactions on Modeling and Computer Simulation. 15: 362-377

Examples

  require(mgcv)  plot(rbind(0,FFdes(2,TRUE)),xlab="x",ylab="y",       col=c(2,1,1,1,1,4,4,4,4),pch=19,main="CCD")  FFdes(5)  FFdes(5,TRUE)

Neighbourhood Cross Validation

Description

NCV estimates smoothing parameters by optimizing the average ability of a model to predict subsets of data when subsets of data are omitted from fitting. Usually the predicted subset is a subset of the omitted subset. If both subsets are the same single datapoint, and the average is over all datapoints, then NCV is leave-one-out cross validation. QNCV is a quadratic approximation to NCV, guaranteed finite for any family link combination.

In detail, suppose that a model is estimated by minimizing a penalized loss

\sum_i D(y_i,\theta_i) + \sum_j \lambda_j \beta^{\sf T} {S}_j \beta

whereD is a loss (such as a negative log likelihood), dependent on responsey_i and parameter vector\theta_i, which in turn depends on covariates via one or more smooth linear predictors with coefficients\beta. The quadratic penalty terms penalize model complexity:S_j is a known matrix and\lambda_j an unknown smoothing parameter. Given smoothing parameters the penalized loss is readily minimized to estimate\beta.

The smoothing parameters also have to be estimated. To this end, choosek = 1,\ldots,m subsets\alpha(k)\subset \{1,\ldots,n\} and\delta(k)\subset \{1,\ldots,n\}. Usually\delta(k) is a subset of (or equal to)\alpha(k). Let\theta_i^{\alpha(k)} denote the estimate of\theta_i when the points indexed by\alpha(k) are omitted from fitting. Then the NCV criterion

V = \sum_{k=1}^m \sum_{i \in \delta(k)} D(y_i,\theta_i^{\alpha(k)})

is minimized w.r.t. the smoothing parameters,\lambda_j. Ifm=n and\alpha(k)=\delta(k)=k then ordinary leave-one-out cross validation is recovered. This formulation covers many of the variants of cross validation reviewed in Arlot and Celisse (2010), for example.

Except for a quadratic loss,V can not be computed exactly, but it can be computed toO(n^{-2}) accuracy (fixed basis size), by taking single Newton optimization steps from the full data\beta estimates to the equivalent when each\alpha(k) is dropped. This is whatmgcv does. The Newton steps require update of the full model Hessian to the equivalent when each datum is dropped. This can be achieved atO(p^2) cost, wherep is the dimension of\beta. Hence, for example, the ordinary cross validation criterion is computable at theO(np^2) cost of estimating the model given smoothing parameters.

The NCV score computed in this way is optimized using a BFGS quasi-Newton method, adapted to the case in which smoothing parameters tending to infinity may cause indefiniteness.

Forbam(...,discrete=TRUE) NCV can be used to estimate the smoothing parameters of the working penalized weighted linear model. This is typically very expensive relative to REML estimation, but the cost can be mitigated by basing the NCV criterion on a sample (usually random of the neighbourhoods). Note that in the case in which leave-out-neighbourhood CV is used to reduce the effects of autocorrelation, it is important to supply a neighbourhood for each observation, and to supply asample element innei. This ensures that the parameter covariance matrix is estimated correctly.

Spatial and temporal short range autocorrelation

A routine applied problem is that smoothing parameters tend to be underestimated in the presence of un-modelled short range autocorrelation, as the smooths try to fit the local excursions in the data caused by the local autocorrelation. Cross validation will tend to 'fit the noise' when there is autocorellation, since a model that fits the noise in the data correlated with an omitted datum, will also tend to closely fit the noise in the omitted datum, because of the correlation. That is autocorrelation works against the avoidance of overfit that cross validation seeks to achieve.

For short range autocorrelation the problems can be avoided, or at least mitigated, by predicting each datum when all the data in its ‘local’ neighbourhood are omitted. The neighbourhoods being constructed in order that un-modelled correlation is minimized between the point of interest and points outside its neighbourhood. That is we setm=n,\delta(k)=k and\alpha(k) = {\tt nei}(k), wherenei(k) are the indices of the neighbours of pointk. This approach has been known for a long time (e.g. Chu and Marron, 1991; Robert et al. 2017), but was previously rather too expensive for regular use for smoothing parameter estimation.

Specifying the neighbourhoods

The neighbourhood subsets\alpha(k) and\delta(k) have to be supplied togam, and thenei argument does this. It is a list with the following arguments.

• a is the vector of indices to be dropped for each neighbourhood.

• ma gives the end of each neighbourhood. Sonei$a[(nei$ma[j-1]+1):nei$ma[j]] gives the points dropped for the neighbourhoodj: that is\alpha(j).

• d is the vector of indices of points to predict.

• md gives the corresponding endpointsmd. Sonei$d[(nei$md[j-1]+1):nei$md[j]] indexes the points to predict for neighbourhood j: that is\delta(j).

• sample is an optional element used bybam. If it is a single number it gives the number of neighbourhoods to randomly sample to construct the NCV criterion. If it is a vector then it should contain the indices of the neighbourhoods to use.

• jackknife is an optional element used bygam. If supplied andTRUE then variance estimates are based on the raw Jackkife estimate, ifFALSE then on the standard Bayesian results. If not supplied (usual) then an estimator accounting for the neighbourhood structure is used, that largely accounts for any correlation present within neighbourhoods.jackknife is ignored if NCV is being calculated for a model where another method is used for smoothing parameter selection.

Ifnei==NULL (ora orma are missing) then leave-one-out cross validation is used. Ifnei is supplied but NCV is not selected as the smoothing parameter estimation method, then it is simply computed (but not optimized).

Numerical issues

If a model is specified in which some coefficient values,\beta, have non-finite likelihood then the NCV criterion computed with single Newton steps could also be non-finite. A simple fix replaces the NCV criterion with a quadratic approximation to the criterion around the full data fit. The quadratic approximation is always finite. This 'QNCV' is essential for some families, such asgevlss. QNCV is not needed forbam estimation.

Although the leading order cost of NCV is the same as REML or GCV, the actual cost is higher because the dominant operations costs are in matrix-vector, rather than matrix-matrix, operations, so BLAS speed ups are small. However multi-core computing is worthwhile for NCV. See the optionncv.threads ingam.control.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Chu and Marron (1991) Comparison of two bandwidth selectors with dependent errors. The Annals of Statistics. 19, 1906-1918

Arlot, S. and A. Celisse (2010). A survey of cross-validation procedures for model selection. Statistics Surveys 4, 40-79

Roberts et al. (2017) Cross-validation strategies for data with temporal,spatial, hierarchical, or phylogenetic structure. Ecography 40(8), 913-929.

Wood S.N. (2026) On Neighbourhood Cross Validation. Statistical Science (in press).https://arxiv.org/abs/2404.16490

Examples

require(mgcv)nei.cor <- function(h,n) { ## construct nei structure  nei <- list(md=1:n,d=1:n)  nei$ma <- cumsum(c((h+1):(2*h+1),rep(2*h+1,n-2*h-2),(2*h+1):(h+1)))  a0 <- rep(0,0); if (h>0) for (i in 1:h) a0 <- c(a0,1:(h+i))  a1 <- n-a0[length(a0):1]+1  nei$a <- c(a0,1:(2*h+1)+rep(0:(n-2*h-1),each=2*h+1),a1)  nei}set.seed(1)n <- 500;sig <- .6x <- 0:(n-1)/(n-1)f <- sin(4*pi*x)*exp(-x*2)*5/2e <- rnorm(n,0,sig)for (i in 2:n) e[i] <- 0.6*e[i-1] + e[i]y <- f + e ## autocorrelated datanei <- nei.cor(4,n) ## construct neighbourhoods to mitigate b0 <- gam(y~s(x,k=40)) ## GCV based fitgc <- gam.control(ncv.threads=2)b1 <- gam(y~s(x,k=40),method="NCV",nei=nei,control=gc)## use "QNCV", which is identical here...b2 <- gam(y~s(x,k=40),method="QNCV",nei=nei,control=gc)## plot GCV and NCV based fits...f <- f - mean(f)par(mfrow=c(1,2))plot(b0,rug=FALSE,scheme=1);lines(x,f,col=2)plot(b1,rug=FALSE,scheme=1);lines(x,f,col=2)

Prediction methods for smooth terms in a GAM

Description

Takessmooth objects produced bysmooth.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 functionPredictMat.

Usage

Predict.matrix(object,data)Predict.matrix2(object,data)

Arguments

.

Details

Smooth terms in a GAM formula are turned into smooth specification objects of classxx.smooth.spec during processing of the formula. Each of these objects isconverted to a smooth object using an appropriatesmooth.construct function. ThePredict.matrix functions are used to obtain the matrix that will map the parameters associated with a smooth term tothe predicted values for the term at new covariate values.

Note that new smooth classes can be added by writing a newsmooth.construct method function and a correspondingPredict.matrix method function: see the example code provided forsmooth.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 inobject. If the smooth classis one which generates offsets the corresponding offset is returned asattribute"offset" of the matrix.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

See Also

gam,gamm,smooth.construct,PredictMat

Examples

# See smooth.construct examples

Predict matrix method functions

Description

The various built in smooth classes for use withgam have associatePredict.matrix method functions to enable prediction from the fitted model.

Usage

## S3 method for class 'cr.smooth'Predict.matrix(object, data)## S3 method for class 'cs.smooth'Predict.matrix(object, data)## S3 method for class 'cyclic.smooth'Predict.matrix(object, data)## S3 method for class 'pspline.smooth'Predict.matrix(object, data)## S3 method for class 'tensor.smooth'Predict.matrix(object, data)## S3 method for class 'tprs.smooth'Predict.matrix(object, data)## S3 method for class 'ts.smooth'Predict.matrix(object, data)## S3 method for class 't2.smooth'Predict.matrix(object, data)

Arguments

.

Details

The Predict matrix function is not normally called directly, but is rather used internally bypredict.gam etc. to predict from a fittedgam model. SeePredict.matrix for more details, or the specificsmooth.construct pages for details on a particular smooth class.

Value

A matrix mapping the coeffients for the smooth term to its values at the supplied data values.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

Examples

## see smooth.construct

Prediction matrix for soap film smooth

Description

Creates a prediction matrix for a soap film smooth object,mapping the coefficients of the smooth to the linear predictor component forthe smooth. This is thePredict.matrix method function required bygam.

Usage

## S3 method for class 'soap.film'Predict.matrix(object,data)## S3 method for class 'sf.wiggly'Predict.matrix(object,data)## S3 method for class 'sf.film'Predict.matrix(object,data)

Arguments

Details

The smooth object will be largely what is returned fromsmooth.construct.so.smooth.spec, although elementsX andS are not needed, and need not be present, of course.

Value

A matrix. This may have an"offset" attribute corresponding tothe contribution from any known boundary conditions on the smooth.

See Also

smooth.construct.so.smooth.spec

Examples

## This is a lower level example. The basis and ## penalties are obtained explicitly ## and `magic' is used as the fitting routine...require(mgcv)set.seed(66)## create a boundary...fsb <- list(fs.boundary())## create some internal knots...knots <- data.frame(x=rep(seq(-.5,3,by=.5),4),                    y=rep(c(-.6,-.3,.3,.6),rep(8,4)))## Simulate some fitting data, inside boundary...n<-1000x <- runif(n)*5-1;y<-runif(n)*2-1z <- fs.test(x,y,b=1)ind <- inSide(fsb,x,y) ## remove outsidersz <- z[ind];x <- x[ind]; y <- y[ind] n <- length(z)z <- z + rnorm(n)*.3 ## add noise## plot boundary with knot and data locationsplot(fsb[[1]]$x,fsb[[1]]$y,type="l");points(knots$x,knots$y,pch=20,col=2)points(x,y,pch=".",col=3);## set up the basis and penalties...sob <- smooth.construct2(s(x,y,bs="so",k=40,xt=list(bnd=fsb,nmax=100)),              data=data.frame(x=x,y=y),knots=knots)## ... model matrix is element `X' of sob, penalties matrices ## are in list element `S'.## fit using `magic'um <- magic(z,sob$X,sp=c(-1,-1),sob$S,off=c(1,1))beta <- um$b## produce plots...par(mfrow=c(2,2),mar=c(4,4,1,1))m<-100;n<-50 xm <- seq(-1,3.5,length=m);yn<-seq(-1,1,length=n)xx <- rep(xm,n);yy<-rep(yn,rep(m,n))## plot truth...tru <- matrix(fs.test(xx,yy),m,n) ## truthimage(xm,yn,tru,col=heat.colors(100),xlab="x",ylab="y")lines(fsb[[1]]$x,fsb[[1]]$y,lwd=3)contour(xm,yn,tru,levels=seq(-5,5,by=.25),add=TRUE)## Plot soap, by first predicting on a fine grid...## First get prediction matrix...X <- Predict.matrix2(sob,data=list(x=xx,y=yy))## Now the predictions...fv <- X%*%beta## Plot the estimated function...image(xm,yn,matrix(fv,m,n),col=heat.colors(100),xlab="x",ylab="y")lines(fsb[[1]]$x,fsb[[1]]$y,lwd=3)points(x,y,pch=".")contour(xm,yn,matrix(fv,m,n),levels=seq(-5,5,by=.25),add=TRUE)## Plot TPRS...b <- gam(z~s(x,y,k=100))fv.gam <- predict(b,newdata=data.frame(x=xx,y=yy))names(sob$sd$bnd[[1]]) <- c("xx","yy","d")ind <- inSide(sob$sd$bnd,xx,yy)fv.gam[!ind]<-NAimage(xm,yn,matrix(fv.gam,m,n),col=heat.colors(100),xlab="x",ylab="y")lines(fsb[[1]]$x,fsb[[1]]$y,lwd=3)points(x,y,pch=".")contour(xm,yn,matrix(fv.gam,m,n),levels=seq(-5,5,by=.25),add=TRUE)

Find rank of upper triangular matrix

Description

Finds rank of upper triangular matrix R, by estimating conditionnumber of upperrank byrank block, and reducingrank until this is acceptably low. Assumes R has been computed by a method that uses pivoting, usually pivoted QR or Choleski.

Usage

Rrank(R,tol=.Machine$double.eps^.9)

Arguments

Details

The method is based on Cline et al. (1979) as described in Golub and van Loan (1996).

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Cline, A.K., C.B. Moler, G.W. Stewart and J.H. Wilkinson (1979) An estimate for the condition number of a matrix. SIAM J. Num. Anal. 16, 368-375

Golub, G.H, and C.F. van Loan (1996) Matrix Computations 3rd ed. Johns Hopkins University Press, Baltimore.

Examples

  set.seed(0)  n <- 10;p <- 5  x <- runif(n*(p-1))  X <- matrix(c(x,x[1:n]),n,p)  qrx <- qr(X,LAPACK=TRUE)  Rrank(qr.R(qrx))

Re-parametrizing model matrix X

Description

INTERNAL routine to apply initial Sl re-parameterization to model matrix X,or, ifinverse==TRUE, to apply inverse re-parametrization to parameter vector or covariance matrix.

Usage

Sl.inirep(Sl,X,l,r,nt=1)Sl.initial.repara(Sl, X, inverse = FALSE, both.sides = TRUE, cov = TRUE,  nt = 1)

Arguments

Value

A re-parametrized version ofX.

Author(s)

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

Applying re-parameterization from log-determinant of penalty matrix tomodel matrix.

Description

INTERNAL routine to apply re-parameterization from log-determinant of penalty matrix,ldetS tomodel matrix,X, blockwise.

Usage

Sl.repara(rp, X, inverse = FALSE, both.sides = TRUE)

Arguments

Value

A re-parametrized version ofX.

Author(s)

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

Setting up a list representing a block diagonal penalty matrix

Description

INTERNAL function for setting up a list representing a block diagonal penalty matrixfrom the object produced bygam.setup.

Usage

Sl.setup(G,cholesky=FALSE,no.repara=FALSE,sparse=FALSE,keepS=FALSE)

Arguments

Value

A list with an element for each block. For block, b,Sl[[b]] is a list with the following elements

• repara: should re-parameterization be applied to model matrix, etc?UsuallyFALSE if non-linear in coefficients.

• start, stop: such thatstart:stop are the indexes of the parameters of this block.

• S: a list of penalty matrices for the block (dim = stop-start+1)Iflength(S)==1 then this will be an identity penalty.Otherwise it is a multiple penalty, and anrS list of squareroot penalty matrices will be added.S (ifrepara==TRUE) andrS (always)will be projected into range space of total penalty matrix.

• rS: square root of penalty matrices if multiple penalties are used.

• D: a reparameterization matrix for the block. Applies to cols/params instart:stop.If numeric thenX[,start:stop]%*%diag(D) is re-parametrization ofX[,start:stop],andb.orig = D*b.repara (whereb.orig is the original parameter vector).If matrix thenX[,start:stop]%*%D is re-parametrization ofX[,start:stop],andb.orig = D%*%b.repara (whereb.orig is the original parameter vector).

• S0: original penalties if requested.

Author(s)

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

GAM Tweedie families

Description

Tweedie families, designed for use withgam from themgcv library.Restricted to variance function powers between 1 and 2. A useful alternative toquasi when afull likelihood is desirable.Tweedie is for use with fixedp.tw is for use whenpis to be estimated during fitting. For fixedp between 1 and 2 the Tweedie is an exponential family distribution with variance given by the mean to the powerp.

tw is only useable withgam andbam but notgamm.Tweedie works with all three.

Usage

Tweedie(p=1, link = power(0))tw(theta = NULL, link = "log",a=1.01,b=1.99)

Arguments

Details

A Tweedie random variable with 1<p<2 is a sum ofN gamma random variables whereN has a Poisson distribution. The p=1 case is a generalization of a Poisson distribution and is a discrete distribution supported on integer multiples of the scale parameter. For 1<p<2 the distribution is supported on the positive reals with a point mass at zero. p=2 is a gamma distribution. As p gets very close to 1 the continuous distribution begins to converge on the discretely supported limit at p=1, and is therefore highly multimodal. SeeldTweedie for more on this behaviour.

Tweedie is based partly on thepoisson family, and partly ontweedie from thestatmod package. It includes extra components to work with allmgcv GAM fitting methods as well as anaic function.

The Tweedie density involves a normalizing constant with no closed form, so this is evaluated using the series evaluation method of Dunn and Smyth (2005), with extensions to also compute the derivatives w.r.t.p and the scale parameter. Without restrictingp to (1,2) the calculation of Tweedie densities is more difficult, and there does not currently seem to be an implementation which offers any benefit overquasi. If you need this case then thetweedie package is the place to start.

Value

ForTweedie, an object inheriting from classfamily, with additional elements

Fortw, an object of classextended.family.

WARNINGS

Asp and the scale parameter tend to 1 the Tweedie distribution tends towards a Poisson, however it remains continuously supported, meaning that the corresponding log likelihood does not converge on the Poisson log likelihood. For example, AIC can be radically different for Poisson and Tweedie models for which the fit is almost identical.

Author(s)

Simon N. Woodsimon.wood@r-project.org.

References

Dunn, P.K. and G.K. Smyth (2005) Series evaluation of Tweedie exponential dispersion model densities. Statistics and Computing 15:267-280

Tweedie, M. C. K. (1984). An index which distinguishes betweensome important exponential families. Statistics: Applications andNew Directions. Proceedings of the Indian Statistical InstituteGolden Jubilee International Conference (Eds. J. K. Ghosh and J.Roy), pp. 579-604. Calcutta: Indian Statistical Institute.

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

See Also

ldTweedie,rTweedie

Examples

library(mgcv)set.seed(3)n<-400## Simulate data...dat <- gamSim(1,n=n,dist="poisson",scale=.2)dat$y <- rTweedie(exp(dat$f),p=1.3,phi=.5) ## Tweedie response## Fit a fixed p Tweedie, with wrong link ...b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=Tweedie(1.25,power(.1)),         data=dat)plot(b,pages=1)print(b)## Same by approximate REML...b1 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=Tweedie(1.25,power(.1)),          data=dat,method="REML")plot(b1,pages=1)print(b1)## estimate p as part of fittingb2 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=tw(),          data=dat,method="REML")plot(b2,pages=1)print(b2)rm(dat)

Internal functions for discretized model matrix handling

Description

Routines for computing with discretized model matrices as described in Wood et al. (2017) and Li and Wood (2019).

Usage

XWXd(X,w,k,ks,ts,dt,v,qc,nthreads=1,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,     lt=NULL,rt=NULL)XWyd(X,w,y,k,ks,ts,dt,v,qc,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,lt=NULL)Xbd(X,beta,k,ks,ts,dt,v,qc,drop=NULL,lt=NULL)diagXVXd(X,V,k,ks,ts,dt,v,qc,drop=NULL,nthreads=1,lt=NULL,rt=NULL)ijXVXd(i,j,X,V,k,ks,ts,dt,v,qc,drop=NULL,nthreads=1,lt=NULL,rt=NULL)

Arguments

Details

These functions are really intended to be internal, but are exported so that they can be used in the initialization code of families without problem. They are primarily used bybam to implement the methods given in the references.XWXd producesX^TWX,XWy producesX^TWy,Xbd producesX\beta anddiagXVXd produces the diagonal ofXVX^T, whileijXVXd evaluates the scatteredi,j elements indexed ini andj.

The"lpip" attribute ofX is a list of the coefficient indices for each term. Required if subsetting vialt andrt.

X can be a list of sparse matrices of class"dgCMatrix", in which case reverse indices are needed, mapping stored matrix rows to rows in the full matrix (that is the reverse ofk which maps full matrix rows to the stored unique matrix rows).r is the same dimension ask whileoff is a list with as many elements ask has columns.r andoff are supplied as attributes toX . For simplicity letr andoff denote a single column and element corresponding to each other: thenr[off[j]:(off[j+1]-1)] contains the rows of the full matrix corresponding to rowj of the stored matrix. The reverse indices are essential for efficient computation with sparse matrices. See the example code for how to create them efficiently from the forward index matrix,k.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association. 112(519):1199-1210doi:10.1080/01621459.2016.1195744

Li, Z & S.N. Wood (2019) Faster model matrix crossproducts for large generalized linear models with discretized covariates. Statistics and Computing.doi:10.1007/s11222-019-09864-2

Examples

  library(mgcv);library(Matrix)  ## simulate some data creating a marginal matrix sequence...  set.seed(0);n <- 4000  dat <- gamSim(1,n=n,dist="normal",scale=2)  dat$x4 <- runif(n)  dat$y <- dat$y + 3*exp(dat$x4*15-5)/(1+exp(dat$x4*15-5))  dat$fac <- factor(sample(1:20,n,replace=TRUE))  G <- gam(y ~ te(x0,x2,k=5,bs="bs",m=1)+s(x1)+s(x4)+s(x3,fac,bs="fs"),           fit=FALSE,data=dat,discrete=TRUE)  p <- ncol(G$X)  ## create a sparse version...  Xs <- list(); r <- G$kd*0; off <- list()  for (i in 1:length(G$Xd)) Xs[[i]] <- as(G$Xd[[i]],"dgCMatrix")  for (j in 1:nrow(G$ks)) { ## create the reverse indices...    nr <- nrow(Xs[[j]]) ## make sure we always tab to final stored row     for (i in G$ks[j,1]:(G$ks[j,2]-1)) {      r[,i] <- (1:length(G$kd[,i]))[order(G$kd[,i])]      off[[i]] <- cumsum(c(1,tabulate(G$kd[,i],nbins=nr)))-1    }  }  attr(Xs,"off") <- off;attr(Xs,"r") <- r   par(mfrow=c(2,3))  beta <- runif(p)  Xb0 <- Xbd(G$Xd,beta,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  Xb1 <- Xbd(Xs,beta,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  range(Xb0-Xb1);plot(Xb0,Xb1,pch=".")  bb <- cbind(beta,beta+runif(p)*.3)  Xb0 <- Xbd(G$Xd,bb,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  Xb1 <- Xbd(Xs,bb,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  range(Xb0-Xb1);plot(Xb0,Xb1,pch=".")  p <- length(beta) ## extract full model matrix...  X <- matrix(Xbd(G$Xd,diag(p),G$kd,G$ks,G$ts,G$dt,G$v,G$qc),ncol=p)  w <- runif(n)  XWy0 <- XWyd(G$Xd,w,y=dat$y,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)   XWy1 <- XWyd(Xs,w,y=dat$y,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  range(XWy1-XWy0);plot(XWy1,XWy0,pch=".")  yy <- cbind(dat$y,dat$y+runif(n)-.5)  XWy0 <- XWyd(G$Xd,w,y=yy,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)   XWy1 <- XWyd(Xs,w,y=yy,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  range(XWy1-XWy0);plot(XWy1,XWy0,pch=".")  A <- XWXd(G$Xd,w,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  B <- XWXd(Xs,w,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  D <- crossprod(X,w*X) ## direct computation  range(A-D)  range(A-B);plot(A,B,pch=".")  ## compute some cross product terms only...  A <- XWXd(G$Xd,w,G$kd,G$ks,G$ts,G$dt,G$v,G$qc,lt=1:3,rt=4:5)  range(A-D[1:nrow(A),(nrow(A)+1):ncol(D)])  V <- crossprod(matrix(runif(p*p),p,p))  ii <- c(20:30,100:200)  jj <- c(50:90,150:160)  V[ii,jj] <- 0;V[jj,ii] <- 0  d1 <- diagXVXd(G$Xd,V,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  Vs <- as(V,"dgCMatrix")  d2 <- diagXVXd(Xs,Vs,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)  range(d1-d2);plot(d1,d2,pch=".")

Approximate hypothesis tests related to GAM fits

Description

Performs hypothesis tests relating to one or more fittedgam objects. For a single fittedgam object, Wald tests ofthe significance of each parametric and smooth term are performed, so interpretation is analogous todrop1 rather thananova.lm (i.e. it's like type III ANOVA, rather than a sequential type I ANOVA). Otherwise the fitted models are compared using an analysis of deviance table or GLRT test: this latter approach should not be use to test the significance of terms which can be penalized to zero. Models to be compared should be fitted to the same data using the same smoothing parameter selection method.

Usage

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

Arguments

Details

If more than one fitted model is provided thananova.glm isused, with the difference in model degrees of freedom being taken as the difference in effective degress of freedom (when possible this is a smoothing parameter uncertainty corrected version).For extended and general families this is set so that a GLRT test is used. The p-values resulting from the multi-model case are only approximate, and must be used with care. The approximation is most accurate when the comparison relates to unpenalized terms, or smoothers with a null space of dimension greater than zero.(Basically we require that the difference terms could be well approximated by unpenalized terms with degrees of freedom approximately the effective degrees of freedom). In simulations the p-values are usually slightly too low. For terms with a zero-dimensional null space (i.e. those which can be penalized to zero) the approximation is often very poor, and significance can be greatly overstated: i.e. p-values are often substantially too low. This case applies to random effect terms.

Note also that in the multi-model call toanova.gam, it is quite possible for a model with more terms to end up with lower effective degrees of freedom, but better fit, than the notionally null model with fewer terms. In such cases it is very rare that it makes sense to perform any sort of test, since there is then no basis on which to accept the notional null model.

If only one model is provided then the significance of each model termis assessed using Wald like tests, conditional on the smoothing parameter estimates: seesummary.gam and Wood (2013a,b) for details. The p-values provided here are better justified than in the multi model case, and have close to the correct distribution under the null, unless smoothing parameters are poorly identified. ML or REML smoothing parameter selection leads to the best results in simulations as they tend to avoid occasional severe undersmoothing. In replication of the full simulation study of Scheipl et al. (2008) the tests give almost indistinguishable power to the method recommended there, but slightly too low p-values under the null in their section 3.1.8 test for a smooth interaction (the Scheipl et al. recommendation is not used directly, because it only applies in the Gaussian case, and requires model refits, but it is available in packageRLRsim).

In the single model caseprint.anova.gam is used as the printing method.

By default the p-values for parametric model terms are also based on Wald tests using the Bayesian covariance matrix for the coefficients. This is appropriate when there are "re" terms present, and is otherwise rather similar to the results using the frequentist covariance matrix (freq=TRUE), since the parametric terms themselves are usually unpenalized. Default P-values for parameteric terms that are penalized using theparaPen argument will not be good.

Value

In the multi-model caseanova.gam produces output identical toanova.glm, which it in fact uses.

In the single model case an object of classanova.gam is produced,which is in fact an object returned fromsummary.gam.

print.anova.gam simply produces tabulated output.

WARNING

If models 'a' and 'b' differ only in terms with no un-penalized components (such as random effects) then p values from anova(a,b) are unreliable, and usually much too low.

Default P-values will usually be wrong for parametric terms penalized using ‘paraPen’: use freq=TRUEto obtain better p-values when the penalties are full rank and represent conventional random effects.

For a single model, interpretation is similar to drop1, not anova.lm.

Author(s)

Simon N. Woodsimon.wood@r-project.org with substantialimprovements by Henric Nilsson.

References

Scheipl, F., Greven, S. and Kuchenhoff, H. (2008) Size and power of tests for a zero random effect variance or polynomial regression in additive and linear mixed models. Comp. Statist. Data Anal. 52, 3283-3299

Wood, S.N. (2013a) On p-values for smooth components of an extended generalized additive model. Biometrika 100:221-228doi:10.1093/biomet/ass048

Wood, S.N. (2013b) A simple test for random effects in regression models. Biometrika 100:1005-1010doi:10.1093/biomet/ast038

See Also

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

Examples

library(mgcv)set.seed(0)dat <- gamSim(5,n=200,scale=2)b<-gam(y ~ x0 + s(x1) + s(x2) + s(x3),data=dat)anova(b)b1<-gam(y ~ x0 + s(x1) + s(x2),data=dat)anova(b,b1,test="F")

Generalized additive models for very large datasets

Description

Fits a generalized additive model (GAM) to a very largedata set, the term ‘GAM’ being taken to include any quadratically penalized GLM (the extended familieslisted infamily.mgcv can also be used). The degree of smoothness of model terms is estimated as part offitting. In use the function is much likegam, except that the numerical methodsare designed for datasets containing upwards of several tens of thousands of data (see Wood, Goude and Shaw, 2015). The advantage ofbam is much lower memory footprint thangam, but it can also be much faster, for large datasets.bam can also compute on a cluster set up by theparallel package.

An alternative fitting approach (Wood et al. 2017, Li and Wood, 2019) is provided by thediscrete==TRUE method. In this case a method based on discretization of covariate values and C code level parallelization (controlled by thenthreads argument instead of thecluster argument) is used. This extends both the data set and model size that are practical. Number of response data can not exceed.Machine$integer.max.

Usage

bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,    na.action=na.omit, offset=NULL,method="fREML",control=list(),    select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL,    paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE,    cluster=NULL,nthreads=1,gc.level=0,use.chol=FALSE,samfrac=1,    coef=NULL,drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,    in.out=NULL,nei=NULL,...)

Arguments

Details

Whendiscrete=FALSE,bam operates by first setting up the basis characteristics for the smooths, using a representative subsample of the data. Then the model matrix is constructed in blocks usingpredict.gam. For each block the factor R, from the QR decomposition of the whole model matrix is updated, along with Q'y. and the sum of squares of y. At the end of block processing, fitting takes place, without the need to ever form the whole model matrix.

In the generalized case, the same trick is used with the weighted model matrix and weighted pseudodata, at each step of the PIRLS.Smoothness selection is performed on the working model at each stage (performance oriented iteration), to maintain the small memory footprint. This is trivial to justify in the case of GCV or Cp/UBRE/AIC based model selection, and for REML/ML is justified via the asymptotic multivariate normality of Q'z where z is the IRLS pseudodata.

For full method details see Wood, Goude and Shaw (2015).

Note that POI is not as stable as the default nested iteration used withgam, but that for very large, information rich,datasets, this is unlikely to matter much.

Note also that it is possible to spend most of the computational time on basis evaluation, if an expensive basis is used. In practice this means that the default"tp" basis should be avoided: almost any other basis (e.g."cr" or"ps") can be used in the 1D case, and tensor product smooths (te) are typically much less costly in the multi-dimensional case.

Ifcluster is provided as a cluster set up usingmakeCluster (ormakeForkCluster) from theparallel package, then the rate limiting QR decomposition of the model matrix is performed in parallel using this cluster. Note that the speed ups are often not that great. On a multi-core machine it is usually best to set the cluster size to the number of physical cores, which is often less than what is reported bydetectCores. Using more than the number of physical cores can result in no speed up at all (or even a slow down). Note that a highly parallel BLAS may negate all advantage from using a cluster of cores. Computing in parallel of course requires more memory than computing in series. See examples.

Whendiscrete=TRUE the covariate data are first discretized. Discretization takes place on a smooth by smooth basis, or in the case of tensor product smooths (or any smooth that can be represented as such, such as random effects), separately for each marginal smooth. The required spline bases are then evaluated at the discrete values, and stored, along with index vectors indicating which original observation they relate to. Fitting is by a version of performance oriented iteration/PQL using REML smoothing parameter selection on each iterative working model (as for the default method). The iteration is based on the derivatives of the REML score, without computing the score itself, allowing the expensive computations to be reduced to one parallel block Cholesky decomposition per iteration (plus two basic operations of equal cost, but easily parallelized). Unlike standard POI/PQL, only one step of the smoothing parameter update for the working model is taken at each step (rather than iterating to the optimal set of smoothing parameters for each working model). At each step a weighted model matrix crossproduct of the model matrix is required - this is efficiently computed from the pre-computed basis functions evaluated at the discretized covariate values. Efficient computation with tensor product terms means that some terms within a tensor product may be re-ordered for maximum efficiency. See Wood et al (2017) and Li and Wood (2019) for full details.

Withdiscrete=TRUE NCV can be used for smoothing parameter estimation. SeeNCV for details of how to set up neighbourhoods using thenei argument. NCV is applied to the working penalized linear model using in iterative fitting. The computation is more efficient than that used bygam provided neighbourhoods are small, but is still costly. It is not efficient for k-fold cross validation, for example. The small neighbourhood restriction is because a matrix of side length given by the neibourhood side needs to be inverted for each neighbourhood. Because computational cost is typically much higher than for REML, for large data/models it is usually worth basing the NCV criterion on a sub-sample of the full data set. Supplying asample element of thenei list is the way to do this. SeeNCV.

Whendiscrete=TRUE parallel computation is controlled using thenthreads argument. For this method no cluster computation is used, and theparallel package is not required. Note that actual speed up from parallelization depends on the BLAS installed and your hardware. With the (R default) reference BLAS using several threads can make a substantial difference, but with a single threaded tuned BLAS, such as openblas, the effect is less marked (since cache use is typically optimized for one thread, and is then sub optimal for several). However the tuned BLAS is usually much faster than using the reference BLAS, however many threads you use. If you have a multi-threaded BLAS installed then you should leaventhreads at 1, since calling a multi-threaded BLAS from multiple threads usually slows things down: the only exception to this is that you might choose to form discrete matrix cross products (the main cost in the fitting routine) in a multi-threaded way, but use single threaded code for other computations: this can be achieved by e.g.nthreads=c(2,1), which would use 2 threads for discrete inner products, and 1 for most code calling BLAS. Not that the basic reason that multi-threaded performance is often disappointing is that most computers are heavily memory bandwidth limited, not flop rate limited. It is hard to get data to one core fast enough, let alone trying to get data simultaneously to several cores.

discrete=TRUE will often produce identical results to the methods without discretization, since covariates often only take a modest number of discrete values anyway, so no approximation at all is involved in the discretization process. Even when some approximation is involved, the differences are often very small as the algorithms discretize marginally whenever possible. For example each margin of a tensor product smooth is discretized separately, rather than discretizing onto a grid of covariate values (for an equivalent isotropic smooth we would have to discretize onto a grid). The marginal approach allows quite fine scale discretization and hence very low approximation error. Note that when using the smoothid mechanism to link smoothing parameters, the discrete method cannot force the linked bases to be identical, so some differences to the none discrete methods will be noticable.

The extended families given infamily.mgcv can also be used. The extra parameters of these are estimated by maximizing the penalized likelihood, rather than the restricted marginal likelihood as ingam. So estimates may differ slightly from those returned bygam. Estimation is accomplished by a Newton iteration to find the extra parameters (e.g. the theta parameter of the negative binomial or the degrees of freedom and scale of the scaled t) maximizing the log likelihood given the model coefficients at each iteration of the fitting procedure.

Value

An object of class"gam" as described ingamObject.

WARNINGS

Formethod set to"fREML","REML" or"ML" the returned smoothness criterion value is for the model itself, otherwise for the working model. The former is comparable togam (but typically a little larger, since it was the working model version that was actually optimized), the latter not.

The routine may be slower than optimal if the default"tp" basis is used.

This routine is less stable than ‘gam’ for the same dataset.

Withdiscrete=TRUE,te terms are efficiently computed, butt2 are not.

Anything close to the maximum n of.Machine$integer.max will need a very large amount of RAM and probablygc.level=1.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N., Goude, Y. & Shaw S. (2015) Generalized additive models for large datasets. Journal of the Royal Statistical Society, Series C 64(1): 139-155.doi:10.1111/rssc.12068

Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association. 112(519):1199-1210doi:10.1080/01621459.2016.1195744

Li, Z & S.N. Wood (2020) Faster model matrix crossproducts for large generalized linear models with discretized covariates. Statistics and Computing. 30:19-25doi:10.1007/s11222-019-09864-2

See Also

mgcv.parallel,mgcv-package,gamObject,gam.models,smooth.terms,linear.functional.terms,s,tepredict.gam,plot.gam,summary.gam,gam.side,gam.selection,gam.controlgam.check,linear.functional.termsnegbin,magic,vis.gam

Examples

library(mgcv)## See help("mgcv-parallel") for using bam in parallel## Sample sizes are small for fast run times.set.seed(3)dat <- gamSim(1,n=25000,dist="normal",scale=20)bs <- "cr";k <- 12b <- bam(y ~ s(x0,bs=bs)+s(x1,bs=bs)+s(x2,bs=bs,k=k)+           s(x3,bs=bs),data=dat)summary(b)plot(b,pages=1,rug=FALSE)  ## plot smooths, but not rugplot(b,pages=1,rug=FALSE,seWithMean=TRUE) ## `with intercept' CIs ba <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+            s(x3,bs=bs,k=k),data=dat,method="GCV.Cp") ## use GCVsummary(ba)## A Poisson example...k <- 15dat <- gamSim(1,n=21000,dist="poisson",scale=.1)system.time(b1 <- bam(y ~ s(x0,bs=bs)+s(x1,bs=bs)+s(x2,bs=bs,k=k),            data=dat,family=poisson()))b1## Similar using faster discrete method... system.time(b2 <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+            s(x3,bs=bs,k=k),data=dat,family=poisson(),discrete=TRUE))b2

Update a strictly additive bam model for new data.

Description

Gaussian with identity link models fitted bybam can be efficiently updated as new data becomes available,by simply updating the QR decomposition on which estimation is based, and re-optimizing the smoothing parameters, startingfrom the previous estimates. This routine implements this.

Usage

bam.update(b,data,chunk.size=10000)

Arguments

Details

bam.update updates the QR decomposition of the (weighted) model matrix of the GAM represented byb to take account of the new data. The orthogonal factor multiplied by the response vector is also updated. Given these updates the model and smoothing parameters can be re-estimated, as if the whole dataset (original and the new data) had been fitted in one go. The function will use the same AR1 model for the residuals as that employed in the original model fit (seerho parameter ofbam).

Note that there may be small numerical differences in fit between fitting the data all at once, and fitting in stages by updating, if the smoothing bases used have any of their details set with reference to the data (e.g. default knot locations).

Value

An object of class"gam" as described ingamObject.

WARNINGS

AIC computation does not currently take account of AR model, if used.

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

mgcv-package,bam

Examples

library(mgcv)## following is not *very* large, for obvious reasons...set.seed(8)n <- 5000dat <- gamSim(1,n=n,dist="normal",scale=5)dat[c(50,13,3000,3005,3100),]<- NAdat1 <- dat[(n-999):n,]dat0 <- dat[1:(n-1000),]bs <- "ps";k <- 20method <- "GCV.Cp"b <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+           s(x3,bs=bs,k=k),data=dat0,method=method)b1 <- bam.update(b,dat1)b2 <- bam.update(bam.update(b,dat1[1:500,]),dat1[501:1000,]) b3 <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+           s(x3,bs=bs,k=k),data=dat,method=method)b1;b2;b3## example with AR1 errors...e <- rnorm(n)for (i in 2:n) e[i] <- e[i-1]*.7 + e[i]dat$y <- dat$f + e*3dat[c(50,13,3000,3005,3100),]<- NAdat1 <- dat[(n-999):n,]dat0 <- dat[1:(n-1000),]b <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+           s(x3,bs=bs,k=k),data=dat0,rho=0.7)b1 <- bam.update(b,dat1)summary(b1);summary(b2);summary(b3)

Choleski decomposition of a band diagonal matrix

Description

Computes Choleski decomposition of a (symmetric positive definite) band-diagonal matrix,A.

Usage

bandchol(B)

Arguments

Details

Callsdpbtrf fromLAPACK. The point of this is that it hasO(k^2n) computational cost, rather than theO(n^3) required by dense matrix methods.

Value

LetR be the factor such thatt(R)%*%R = A.R is upper triangular and if the rows ofB contained the diagonals ofA on entry, then what is returned is an n by k matrix containing the diagonals ofR, packed asB was packed on entry. IfB was square on entry, thenR is returned directly. See examples.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Anderson, E., Bai, Z., Bischof, C., Blackford, S., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. and Sorensen, D., 1999. LAPACK Users' guide (Vol. 9). Siam.

Examples

require(mgcv)## simulate a banded diagonal matrixn <- 7;set.seed(8)A <- matrix(0,n,n)sdiag(A) <- runif(n);sdiag(A,1) <- runif(n-1)sdiag(A,2) <- runif(n-2)A <- crossprod(A) ## full matrix form...bandchol(A)chol(A) ## for comparison## compact storage form...B <- matrix(0,3,n)B[1,] <- sdiag(A);B[2,1:(n-1)] <- sdiag(A,1)B[3,1:(n-2)] <- sdiag(A,2)bandchol(B)

Censored Box-Cox Gaussian family

Description

Family for use withgam orbam, implementing regression for censoreddata that can be modelled as Gaussian after Box-Cox transformation. Ify>0 is the response and

z= \left \{ \begin{array}{ll} (y^\lambda-1)/\lambda & \lambda \ne 0 \\ \log(y) & \lambda =0\end{array} \right .

with mean\mu and standard deviationw^{-1/2}\exp(\theta),thenw^{1/2}(z-\mu)\exp(-\theta) follows anN(0,1) distribution. That is

z \sim N(\mu,e^{2\theta}w^{-1}).

\theta is a single scalar for all observations, and similarly Box-Cox parameter\lambda is a single scalar. Observations may be left, interval or right censored or uncensored.

Note that the regression model here specifies the mean of the Box-Cox transformed response, not the mean of the response itself: this is rather different to the usual GLM approach.

Usage

bcg(theta=NULL,link="identity")

Arguments

Details

If the family is used with a vector response, then it is assumed that there is no censoring. If there is censoring then the response should be supplied as a two column matrix. The first column is always numeric. Entries in the second column are as follows.

• If an entry is identical to the corresponding first column entry, then it is an uncensored observation.

• If an entry is numeric and different to the first column entry then there is interval censoring. The first column entry is the lower interval limit and the second column entry is the upper interval limit.y is only known to be between these limits.

• If the second column entry is0 then the observation is left censored at the value of the entry in the first column. It is only known thaty is less than or equal to the first column value.

• If the second column entry isInf then the observation is right censored at the value of the entry in the first column. It is only known thaty is greater than or equal to the first column value.

Any mixture of censored and uncensored data is allowed, but be aware that data consisting only of right and/or left censored data contain very little information.

Value

An object of classextended.family.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

See Also

cnorm,cpois,clog

Examples

library(mgcv)set.seed(3) ## Simulate some gamma data?dat <- gamSim(1,n=400,dist="normal",scale=1)dat$f <- dat$f/4 ## true linear predictor Ey <- exp(dat$f);scale <- .5 ## mean and GLM scale parameterdat$y <- rgamma(Ey*0,shape=1/scale,scale=Ey*scale)## Just Box-Cox no censoring...b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=bcg,data=dat)summary(b)plot(b,pages=1,scheme=1)## try various censoring...yb <- cbind(dat$y,dat$y)ind <- 1:100yb[ind,2] <- yb[ind,1] + runif(100)*3yb[51:100,2] <- 0 ## left censoredyb[101:140,2] <- Inf ## right censoredb <- gam(yb~s(x0)+s(x1)+s(x2)+s(x3),family=bcg,data=dat)summary(b)plot(b,pages=1,scheme=1)

GAM beta regression family

Description

Family for use withgam orbam, implementing regression for beta distributed data on (0,1).A linear predictor controls the mean,\mu of the beta distribution, while the variance is then\mu(1-\mu)/(1+\phi), with parameter\phi being estimated during fitting, alongside the smoothing parameters.

Usage

betar(theta = NULL, link = "logit",eps=.Machine$double.eps*100)

Arguments

Details

These models are useful for proportions data which can not be modelled as binomial. Note the assumption that data are in (0,1), despite the fact that for some parameter values 0 and 1 are perfectly legitimate observations. The restriction is needed to keep the log likelihood bounded for all parameter values. Any data exactly at 0 or 1 are reset to be just above 0 or just below 1 using theeps argument (in fact any observation<eps is reset toeps and any observation>1-eps is reset to1-eps). Note the effect of this resetting. If\mu\phi>1 then impossible 0s are replaced with highly improbableeps values. If the inequality is reversed then 0s with infinite probability density are replaced witheps values having high finite probability density. The equivalent condition for 1s is(1-\mu)\phi>1. Clearly all types of resetting are somewhat unsatisfactory, and care is needed if data contain 0s or 1s (often it makes sense to manually reset the 0s and 1s in a manner that somehow reflects the sampling setup).

Value

An object of classextended.family.

WARNINGS

Do read the details section if your data contain 0s and or 1s.

Author(s)

Natalya Pya (nat.pya@gmail.com) and Simon Wood (s.wood@r-project.org)

Examples

library(mgcv)## Simulate some beta data...set.seed(3);n<-400dat <- gamSim(1,n=n)mu <- binomial()$linkinv(dat$f/4-2)phi <- .5a <- mu*phi;b <- phi - a;dat$y <- rbeta(n,a,b) bm <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=betar(link="logit"),data=dat)bmplot(bm,pages=1)

BLAS thread safety

Description

Most BLAS implementations are thread safe, but some versions of OpenBLAS, for example, are not. This routine is a diagnostic helper function, which you will never need if you don't setnthreads>1, and even then are unlikely to need.

Usage

blas.thread.test(n=1000,nt=4)

Arguments

Details

While single threaded OpenBLAS 0.2.20 was thread safe, versions 0.3.0-0.3.6 are not, and from version 0.3.7 thread safety of the single threaded OpenBLAS requires making it with the optionUSE_LOCKING=1. The reference BLAS is thread safe, as are MKL and ATLAS. This routine repeatedly calls the BLAS from multi-threaded code and is sufficient to detect the problem in single threaded OpenBLAS 0.3.x.

A multi-threaded BLAS is often no faster than a single-threaded BLAS, while judicious use of threading in the code calling the BLAS can still deliver a modest speed improvement. For this reason it is often better to use a single threaded BLAS and thenthreads options tobam orgam. Forbam(...,discrete=TRUE) using several threads can be a substantial benefit, especially with the reference BLAS.

The MKL BLAS is mutlithreaded by default. Under linux setting environment variableMKL_NUM_THREADS=1 before starting R gives single threaded operation.

Author(s)

Simon N. Woodsimon.wood@r-project.org

Reporting mgcv bugs.

Description

mgcv works largely because many people have reported bugs over the years. If you find something that looks like a bug, please report it, so that the package can be improved.mgcv does not have a large development budget, so it is a big help if bug reports follow the following guidelines.

The ideal report consists of an email tosimon.wood@r-project.org with a subject line includingmgcv somewhere, containing

1. The results of runningsessionInfo in the R session where the problem occurs. This provides platform details, R and package version numbers, etc.

2. A brief description of the problem.

3. Short cut and paste-able code that produces the problem, including the code for loading/generating the data (using standard R functions likeload,read.table etc).

4. Any required data files. If you send real data it will only be used for the purposes of de-bugging.

Of course if you have dug deeper and have an idea of what is causing the problem, that is also helpful to know, as is any suggested code fix. (Don't send a fixed package .tar.gz file, however - I can't use this).

Author(s)

Simon N. Woodsimon.wood@r-project.org

Evaluate cyclic B spline basis

Description

UsessplineDesign to set up the model matrix for a cyclic B-spline basis.

Usage

cSplineDes(x, knots, ord = 4, derivs=0)

Arguments

Details

The routine is a wrapper that sets up a B-spline basis, where the basis functions wrap at the first and last knot locations.

Value

A matrix withlength(x) rows andlength(knots)-1 columns.

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

cyclic.p.spline

Examples

 require(mgcv) ## create some x's and knots... n <- 200 x <- 0:(n-1)/(n-1);k<- 0:5/5 X <- cSplineDes(x,k) ## cyclic spline design matrix ## plot evaluated basis functions... plot(x,X[,1],type="l"); for (i in 2:5) lines(x,X[,i],col=i) ## check that the ends match up.... ee <- X[1,]-X[n,];ee  tol <- .Machine$double.eps^.75 if (all.equal(ee,ee*0,tolerance=tol)!=TRUE)    stop("cyclic spline ends don't match!")  ## similar with uneven data spacing... x <- sort(runif(n)) + 1 ## sorting just makes end checking easy k <- seq(min(x),max(x),length=8) ## create knots X <- cSplineDes(x,k) ## get cyclic spline model matrix   plot(x,X[,1],type="l"); for (i in 2:ncol(X)) lines(x,X[,i],col=i) ee <- X[1,]-X[n,];ee ## do ends match?? tol <- .Machine$double.eps^.75 if (all.equal(ee,ee*0,tolerance=tol)!=TRUE)    stop("cyclic spline ends don't match!")

Deletion and rank one Cholesky factor update

Description

Given a Cholesky factor,R, of a matrix,A,choldrop finds the Cholesky factor ofA[-k,-k],wherek is an integer.cholup finds the factor ofA + uu^T (update) orA - uu^T (downdate).

Usage

choldrop(R,k)cholup(R,u,up)

Arguments

Details

First considercholdrop. IfR is upper triangular thent(R[,-k])%*%R[,-k] == A[-k,-k], butR[,-k] has elements on the first sub-diagonal, from its kth column onwards. To get from this to a triangular Cholesky factor ofA[-k,-k] we can apply a sequence of Givens rotations from the left to eliminate the sub-diagonal elements. The routine does this. IfR is a lower triangular factor then Givens rotations from the right are needed to remove the extra elements. Ifn is the dimension ofR then the update hasO(n^2) computational cost.

cholup (which assumesR is upper triangular) updates based on the observation that R^TR + uu^T = [u,R^T][u,R^T]^T = [u,R^T]Q^TQ[u,R^T]^T, and therefore we can constructQ so thatQ[u,R^T]^T=[0,R_1^T]^T, whereR_1 is the modified factor.Q is constructed from a sequence of Givens rotations in order to zero the elements ofu. Downdating is similar except that hyperbolic rotations have to be used in place of Givens rotations — see Golub and van Loan (2013, section 6.5.4) for details. Downdating only works ifA - uu^T is positive definite. Again the computational cost isO(n^2).

Note that the updates are vector oriented, and are hence not susceptible to speed up by use of an optimized BLAS. The updates are set up to be relatively Cache friendly, in that in the upper triangular case successive Givens rotations are stored for sequential application column-wise, rather than being applied row-wise as soon as they are computed. Even so, the upper triangular update is slightly slower than the lower triangular update.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Golub GH and CF Van Loan (2013) Matrix Computations (4th edition) Johns Hopkins

Examples

  require(mgcv)  set.seed(0)  n <- 6  A <- crossprod(matrix(runif(n*n),n,n))  R0 <- chol(A)  k <- 3  Rd <- choldrop(R0,k)  range(Rd-chol(A[-k,-k]))  Rd;chol(A[-k,-k])    ## same but using lower triangular factor A = LL'  L <- t(R0)  Ld <- choldrop(L,k)  range(Ld-t(chol(A[-k,-k])))  Ld;t(chol(A[-k,-k]))  ## Rank one update example  u <- runif(n)  R <- cholup(R0,u,TRUE)  Ru <- chol(A+u %*% t(u)) ## direct for comparison  R;Ru  range(R-Ru)  ## Downdate - just going back from R to R0  Rd <-  cholup(R,u,FALSE)  R0;Rd  range(R0-Rd)

Basis dimension choice for smooths

Description

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

Penalized regression smoothers gain computational efficiency by virtue ofbeing defined using a basis of relatively modest size,k. When settingup models in themgcv package, usings orteterms in a model formula,k must be chosen: the defaults are essentially arbitrary.

In practicek-1 (ork) sets the upper limit on the degrees of freedomassociated with ans smooth (1 degree of freedom is usually lost to the identifiabilityconstraint on the smooth). Forte smooths the upper limit of thedegrees of freedom is given by the product of thek values provided foreach marginal smooth less one, for the constraint. However the actualeffective degrees of freedom are controlled by the degree of penalizationselected during fitting, by GCV, AIC, REML or whatever is specified. The exceptionto this is if a smooth is specified using thefx=TRUE option, in whichcase it is unpenalized.

So, exact choice ofk is not generally critical: it should be chosen tobe large enough that you are reasonably sure of having enough degrees offreedom to represent the underlying ‘truth’ reasonably well, but small enoughto 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 ofk informally. If the effective degrees of freedom for a model term areestimated to be much less thank-1 then this is unlikely to be veryworthwhile, but as the EDF approachk-1, checking can be important. Auseful general purpose approach goes as follows: (i) fit your model andextract the deviance residuals; (ii) for each smooth term in your model, fitan equivalent, single, smooth to the residuals, using a substantiallyincreasedk to see if there is pattern in the residuals that couldpotentially be explained by increasingk. Examples are provided below.

The obvious, but more costly, alternative is simply to increase the suspectk and refit the original model. If there are no statistically important changes as a result of doing this, thenk was large enough. (Change in the smoothness selection criterion, and/or the effective degrees of freedom, whenk is increased, provide the obvious numerical measures for whether the fit has changed substantially.)

gam.check runs a simple simulation based check on the basis dimensions, which can help to flag up terms for whichk is too low. Grossly too smallkwill also be visible from partial residuals available withplot.gam.

One scenario that can cause confusion is this: a model is fitted withk=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 withk=20and the EDF increases to 8.7 - what is happening - how come the EDF was not8.7 the first time around? The explanation is that the function space withk=20 contains a larger subspace of functions with EDF 8.7 than did thefunction space withk=10: one of the functions in this larger subspacefits the data a little better than did any function in the smallersubspace. These subtleties seldom have much impact on the statisticalconclusions to be drawn from a model fit, however.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). CRC/Taylor & Francis.

Examples

## Simulate some data ....library(mgcv)set.seed(1) dat <- gamSim(1,n=400,scale=2)## 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),data=dat)plot(b,pages=1,residuals=TRUE) ## hint of a problem in s(x2)## the following suggests a problem with s(x2)gam.check(b)## Another approach (see below for more obvious method)....## 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,data=dat) ## finegam(rsd~s(x1,k=40,bs="cs"),gamma=1.4,data=dat) ## finegam(rsd~s(x2,k=40,bs="cs"),gamma=1.4,data=dat) ## `k' too lowgam(rsd~s(x3,k=40,bs="cs"),gamma=1.4,data=dat) ## fine## refit...b <- gam(y~s(x0,k=6)+s(x1,k=6)+s(x2,k=20)+s(x3,k=6),data=dat)gam.check(b) ## better## similar example with multi-dimensional smoothb1 <- gam(y~s(x0)+s(x1,x2,k=15)+s(x3),data=dat)rsd <- residuals(b1)gam(rsd~s(x0,k=40,bs="cs"),gamma=1.4,data=dat) ## finegam(rsd~s(x1,x2,k=100,bs="ts"),gamma=1.4,data=dat) ## `k' too lowgam(rsd~s(x3,k=40,bs="cs"),gamma=1.4,data=dat) ## finegam.check(b1) ## shows same problem## and a `te' exampleb2 <- gam(y~s(x0)+te(x1,x2,k=4)+s(x3),data=dat)rsd <- residuals(b2)gam(rsd~s(x0,k=40,bs="cs"),gamma=1.4,data=dat) ## finegam(rsd~te(x1,x2,k=10,bs="cs"),gamma=1.4,data=dat) ## `k' too lowgam(rsd~s(x3,k=40,bs="cs"),gamma=1.4,data=dat) ## finegam.check(b2) ## shows same problem## same approach works with other families in the original modeldat <- gamSim(1,n=400,scale=.25,dist="poisson")bp<-gam(y~s(x0,k=5)+s(x1,k=5)+s(x2,k=5)+s(x3,k=5),        family=poisson,data=dat,method="ML")gam.check(bp)rsd <- residuals(bp)gam(rsd~s(x0,k=40,bs="cs"),gamma=1.4,data=dat) ## finegam(rsd~s(x1,k=40,bs="cs"),gamma=1.4,data=dat) ## finegam(rsd~s(x2,k=40,bs="cs"),gamma=1.4,data=dat) ## `k' too lowgam(rsd~s(x3,k=40,bs="cs"),gamma=1.4,data=dat) ## finerm(dat) ## More obvious, but more expensive tactic... Just increase ## suspicious k until fit is stable.set.seed(0) dat <- gamSim(1,n=400,scale=2)## 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),         data=dat,method="REML")b ## edf for 3rd smooth is highest as proportion of k -- increase kb <- gam(y~s(x0,k=6)+s(x1,k=6)+s(x2,k=12)+s(x3,k=6),         data=dat,method="REML")b ## edf substantially up, -ve REML substantially downb <- gam(y~s(x0,k=6)+s(x1,k=6)+s(x2,k=24)+s(x3,k=6),         data=dat,method="REML")b ## slight edf increase and -ve REML changeb <- gam(y~s(x0,k=6)+s(x1,k=6)+s(x2,k=40)+s(x3,k=6),         data=dat,method="REML")b ## defintely stabilized (but really k around 20 would have been fine)

GAM censored logistic distribution family for log-logistic AFT models

Description

Family for use withgam orbam, implementing regression for censoredlogistic data. Ify is the response with mean\mu and scale parameters=w^{-1/2}\exp(\theta),theny has p.d.f.

\frac{\exp\{-(y-\mu)/s\}}{s[1+\exp\{-(y-\mu)/s\}]^2}.

\theta is a single scalar for all observations. Observations may be left, interval or right censored or uncensored.

Useful for log-logistic accelerated failure time (AFT) models, for example.

Usage

clog(theta=NULL,link="identity")

Arguments

Details

If the family is used with a vector response, then it is assumed that there is no censoring, and a regular Gaussian regression results. If there is censoring then the response should be supplied as a two column matrix. The first column is always numeric. Entries in the second column are as follows.

• If an entry is identical to the corresponding first column entry, then it is an uncensored observation.

• If an entry is numeric and different to the first column entry then there is interval censoring. The first column entry is the lower interval limit and the second column entry is the upper interval limit.y is only known to be between these limits.

• If the second column entry is-Inf then the observation is left censored at the value of the entry in the first column. It is only known thaty is less than or equal to the first column value.

• If the second column entry isInf then the observation is right censored at the value of the entry in the first column. It is only known thaty is greater than or equal to the first column value.

Any mixture of censored and uncensored data is allowed, but be aware that data consisting only of right and/or left censored data contain very little information.

Value

An object of classextended.family.

Author(s)

Chris Shen

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Examples

library(mgcv)######################################################### AFT model example for colon cancer survivial data...#######################################################library(survival) ## for datacol1 <- colon[colon$etype==1,] ## concentrate on single eventcol1$differ <- as.factor(col1$differ)col1$sex <- as.factor(col1$sex)## set up the AFT response... logt <- cbind(log(col1$time),log(col1$time))logt[col1$status==0,2] <- Inf ## right censoringcol1$logt <- -logt ## -ve conventional for AFT versus Cox PH comparison## fit the model...b <- gam(logt~s(age,by=sex)+sex+s(nodes)+perfor+rx+obstruct+adhere,         family=clog(),data=col1)plot(b,pages=1) ## ... compare this to ?cox.ph

GAM censored normal family for log-normal AFT and Tobit models

Description

Family for use withgam orbam, implementing regression for censorednormal data. Ify is the response with mean\mu and standard deviationw^{-1/2}\exp(\theta),thenw^{1/2}(y-\mu)\exp(-\theta) follows anN(0,1) distribution. That is

y \sim N(\mu,e^{2\theta}w^{-1}).

\theta is a single scalar for all observations. Observations may be left, interval or right censored or uncensored.

Useful for log-normal accelerated failure time (AFT) models, Tobit regression, and crudely rounded data, for example.

Usage

cnorm(theta=NULL,link="identity")

Arguments

Details

If the family is used with a vector response, then it is assumed that there is no censoring, and a regular Gaussian regression results. If there is censoring then the response should be supplied as a two column matrix. The first column is always numeric. Entries in the second column are as follows.

• If an entry is identical to the corresponding first column entry, then it is an uncensored observation.

• If an entry is numeric and different to the first column entry then there is interval censoring. The first column entry is the lower interval limit and the second column entry is the upper interval limit.y is only known to be between these limits.

• If the second column entry is-Inf then the observation is left censored at the value of the entry in the first column. It is only known thaty is less than or equal to the first column value.

• If the second column entry isInf then the observation is right censored at the value of the entry in the first column. It is only known thaty is greater than or equal to the first column value.

Any mixture of censored and uncensored data is allowed, but be aware that data consisting only of right and/or left censored data contain very little information.

Value

An object of classextended.family.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Examples

library(mgcv)######################################################### AFT model example for colon cancer survivial data...#######################################################library(survival) ## for datacol1 <- colon[colon$etype==1,] ## concentrate on single eventcol1$differ <- as.factor(col1$differ)col1$sex <- as.factor(col1$sex)## set up the AFT response... logt <- cbind(log(col1$time),log(col1$time))logt[col1$status==0,2] <- Inf ## right censoringcol1$logt <- -logt ## -ve conventional for AFT versus Cox PH comparison## fit the model...b <- gam(logt~s(age,by=sex)+sex+s(nodes)+perfor+rx+obstruct+adhere,         family=cnorm(),data=col1)plot(b,pages=1) ## ... compare this to ?cox.ph################################## A Tobit regression example...################################set.seed(3);n<-400dat <- gamSim(1,n=n)ys <- dat$y - 5 ## shift data down## truncate at zero, and set up response indicating this has happened...y <- cbind(ys,ys)y[ys<0,2] <- -Infy[ys<0,1] <- 0dat$yt <- yb <- gam(yt~s(x0)+s(x1)+s(x2)+s(x3),family=cnorm,data=dat)plot(b,pages=1)################################ A model for rounded data...##############################dat <- gamSim(1,n=n)y <- round(dat$y)y <- cbind(y-.5,y+.5) ## set up to indicate interval censoringdat$yi <- yb <- gam(yi~s(x0)+s(x1)+s(x2)+s(x3),family=cnorm,data=dat)plot(b,pages=1)

Reduced version of Columbus OH crime data

Description

By district crime data from Columbus OH, together with polygons describing district shape. Useful for illustrating use of simple Markov Random Field smoothers.

Usage

data(columb)data(columb.polys)

Format

columb is a 49 row data frame with the following columns

area

land area of district

home.value

housing value in 1000USD.

income

household income in 1000USD.

crime

residential burglaries and auto thefts per 1000 households.

open.space

measure of open space in district.

district

code identifying district, and matchingnames(columb.polys).

columb.polys contains the polygons defining the areas in the format described below.

Details

The data framecolumb relates to the districts whose boundaries are coded incolumb.polys.columb.polys[[i]] is a 2 column matrix, containing the vertices of the polygons defining the boundary of the ith district.columb.polys[[2]] has an artificial hole inserted to illustrate how holes in districts can be spefified. Different polygons defining the boundary of a district are separated by NA rows incolumb.polys[[1]], and a polygon enclosed within another is treated as a hole in that region (a hole should never come first).names(columb.polys) matchescolumb$district (order unimportant).

Source

The data are adapted from thecolumbus example in thespdep package, where the original source is given as:

Anselin, Luc. 1988. Spatial econometrics: methods and models. Dordrecht: Kluwer Academic, Table 12.1 p. 189.

Examples

## see ?mrf help files

GAM concurvity measures

Description

Produces summary measures of concurvity betweengam components.

Usage

concurvity(b,full=TRUE)

Arguments

Details

Concurvity occurs when some smooth term in a model could be approximated by one or more of the other smooth terms in the model. This is often the case when a smooth of space is included in a model, along with smooths of other covariates that also vary more or less smoothly in space. Similarly it tends to be an issue in models including a smooth of time, along with smooths of other time varying covariates.

Concurvity can be viewed as a generalization of co-linearity, and causes similar problems of interpretation. It can also make estimates somewhat unstable (so that they become sensitive to apparently innocuous modelling details, for example).

This routine computes three related indices of concurvity, all bounded between 0 and 1, with 0 indicating no problem, and 1 indicating total lack of identifiability. The three indices are all based on the idea that a smooth term, f, in the model can be decomposed into a part, g, that lies entirely in the space of one or more other terms in the model, and a remainder part that is completely within the term's own space. If g makes up a large part of f then there is a concurvity problem. The indices used are all based on the square of ||g||/||f||, that is the ratio of the squared Euclidean norms of the vectors of f and g evaluated at the observed covariate values.

The three measures are as follows

worst

This is the largest value that the square of ||g||/||f|| could take for any coefficient vector. This is a fairly pessimistic measure, as it looks at the worst case irrespective of data. This is the only measure that is symmetric.

observed

This just returns the value of the square of ||g||/||f|| according to the estimated coefficients. This could be a bit over-optimistic about the potential for a problem in some cases.

estimate

This is the squared F-norm of the basis for g divided by the F-norm of the basis for f. It is a measure of the extent to which the f basis can be explained by the g basis. It does not suffer from the pessimism or potential for over-optimism of the previous two measures, but is less easy to understand.

Value

Iffull=TRUE a matrix with one column for each term and one row for each of the 3 concurvity measures detailed below. Iffull=FALSE a list of 3 matrices, one for each of the three concurvity measures detailed below. Each row of the matrix relates to how the model terms depend on the model term supplying that rows name.

Author(s)

Simon N. Woodsimon.wood@r-project.org

Examples

library(mgcv)## simulate data with concurvity...set.seed(8);n<- 200f2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 *            (10 * x)^3 * (1 - x)^10t <- sort(runif(n)) ## first covariate## make covariate x a smooth function of t + noise...x <- f2(t) + rnorm(n)*3## simulate response dependent on t and x...y <- sin(4*pi*t) + exp(x/20) + rnorm(n)*.3## fit model...b <- gam(y ~ s(t,k=15) + s(x,k=15),method="REML")## assess concurvity between each term and `rest of model'...concurvity(b)## ... and now look at pairwise concurvity between terms...concurvity(b,full=FALSE)

Additive Cox Proportional Hazard Model

Description

Thecox.ph family implements the Cox Proportional Hazards model with Peto's correction for ties, optional stratification, and estimation by penalized partial likelihood maximization, for use withgam. In the model formula, event time is the response. Under stratification the response has two columns: timeand a numeric index for stratum. Theweights vector provides the censoring information (0 for censoring, 1 for event).cox.ph deals with the case in which each subject hasone event/censoring time and one row of covariate values. When each subject has several time dependentcovariates seecox.pht.

See example below for conditional logistic regression.

Usage

cox.ph(link="identity")

Arguments

Details

Used withgam to fit Cox Proportional Hazards models to survival data. The model formula will have event/censoring times on the left hand side and the linear predictor specification on the right hand side. Censoringinformation is provided by theweights argument togam, with 1 indicating an event and 0 indicating censoring.

Stratification is possible, allowing for different baseline hazards in different strata. In that case the response has two columns: the first is event/censoring time and the second is a numeric stratum index. See below for an example.

Prediction from the fitted model object (using thepredict method) withtype="response" will predict on the survivor function scale. This requires evaluation times to be provided as well as covariates (see example). Also see example codebelow for extracting the cumulative baseline hazard/survival directly. Thefitted.values stored in the model object aresurvival function estimates for each subject at their event/censoring time.

deviance,martingale,score, orschoenfeld residuals can be extracted. See Klein amd Moeschberger (2003) for descriptions. The score residuals are returned as a matrix of the same dimension as the model matrix, with a"terms" attribute, which is a list indicating which model matrix columns belong to which model terms. The score residuals are scaled. For parameteric terms this is by the standard deviation of associated model coefficient. For smooth terms the sub matrix of score residuals for the term is postmultiplied by the transposed Cholesky factor of the covariance matrix for the term's coefficients. This is a transformation that makes the coefficients approximately independent, as required to make plots of the score residuals against event time interpretable for checking the proportional hazards assumption (see Klein amd Moeschberger, 2003, p376). Penalization causes drift in the score residuals, which is also removed, to allow the residuals to be approximately interpreted as unpenalized score residuals. Schoenfeld and score residuals are computed by strata. See the examples for simple PH assuption checks by plotting score residuals, and Klein amd Moeschberger (2003, section 11.4) for details. Note that high correlation between terms can undermine these checks.

Estimation of model coefficients is by maximising the log-partial likelihood penalized by the smoothing penalties. See e.g. Hastie and Tibshirani, 1990, section 8.3. for the partial likelihood used (with Peto's approximation for ties), but note that optimization of the partial likelihood does not follow Hastie and Tibshirani. See Klein amd Moeschberger (2003) for estimation of residuals, the cumulative baseline hazard, survival function and associated standard errors (the survival standard error expression has a typo).

The percentage deviance explained reported for Cox PH models is based on the sum of squares of the deviance residuals, as the model deviance, and the sum of squares of the deviance residuals when the covariate effects are set to zero, as the null deviance. The same baseline hazard estimate is used for both.

This family deals efficiently with the case in which each subject has one event/censoring time and one row of covariate values. For studies in which there are multiple time varying covariate measures for each subject then the equivalent Poisson model should be fitted to suitable pseudodata usingbam(...,discrete=TRUE). Seecox.pht.

Value

An object inheriting from classgeneral.family.

References

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

Klein, J.P and Moeschberger, M.L. (2003) Survival Analysis: Techniques forCensored and Truncated Data (2nd ed.) Springer.

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

See Also

cox.pht,cnorm

Examples

library(mgcv)library(survival) ## for datacol1 <- colon[colon$etype==1,] ## concentrate on single eventcol1$differ <- as.factor(col1$differ)col1$sex <- as.factor(col1$sex)b <- gam(time~s(age,by=sex)+sex+s(nodes)+perfor+rx+obstruct+adhere,         family=cox.ph(),data=col1,weights=status)summary(b) plot(b,pages=1,all.terms=TRUE) ## plot effectsplot(b$linear.predictors,residuals(b))## plot survival function for patient j...np <- 300;j <- 6newd <- data.frame(time=seq(0,3000,length=np))dname <- names(col1)for (n in dname) newd[[n]] <- rep(col1[[n]][j],np)newd$time <- seq(0,3000,length=np)fv <- predict(b,newdata=newd,type="response",se=TRUE)plot(newd$time,fv$fit,type="l",ylim=c(0,1),xlab="time",ylab="survival")lines(newd$time,fv$fit+2*fv$se.fit,col=2)lines(newd$time,fv$fit-2*fv$se.fit,col=2)## crude plot of baseline survival...plot(b$family$data$tr,exp(-b$family$data$h),type="l",ylim=c(0,1),     xlab="time",ylab="survival")lines(b$family$data$tr,exp(-b$family$data$h + 2*b$family$data$q^.5),col=2)lines(b$family$data$tr,exp(-b$family$data$h - 2*b$family$data$q^.5),col=2)lines(b$family$data$tr,exp(-b$family$data$km),lty=2) ## Kaplan Meier## Checking the proportional hazards assumption via scaled score plots as## in Klein and Moeschberger Section 11.4 p374-376... ph.resid <- function(b,stratum=1) {## convenience function to plot scaled score residuals against time,## by term. Reference lines at 5% exceedance prob for Brownian bridge## (see KS test statistic distribution).  rs <- residuals(b,"score");term <- attr(rs,"term")  if (is.matrix(b$y)) {    ii <- b$y[,2] == stratum;b$y <- b$y[ii,1];rs <- rs[ii,]  }  oy <- order(b$y)  for (i in 1:length(term)) {    ii <- term[[i]]; m <- length(ii)    plot(b$y[oy],rs[oy,ii[1]],ylim=c(-3,3),type="l",ylab="score residuals",         xlab="time",main=names(term)[i])    if (m>1) for (k in 2:m) lines(b$y[oy],rs[oy,ii[k]],col=k);    abline(-1.3581,0,lty=2);abline(1.3581,0,lty=2)   }  }par(mfrow=c(2,2))ph.resid(b)## stratification example, with 2 randomly allocated strata## so that results should be similar to previous....col1$strata <- sample(1:2,nrow(col1),replace=TRUE) bs <- gam(cbind(time,strata)~s(age,by=sex)+sex+s(nodes)+perfor+rx+obstruct          +adhere,family=cox.ph(),data=col1,weights=status)plot(bs,pages=1,all.terms=TRUE) ## plot effects## baseline survival plots by strata...for (i in 1:2) { ## loop over strata## create index selecting elements of stored hazard info for stratum i...ind <- which(bs$family$data$tr.strat == i)if (i==1) plot(bs$family$data$tr[ind],exp(-bs$family$data$h[ind]),type="l",      ylim=c(0,1),xlab="time",ylab="survival",lwd=2,col=i) else      lines(bs$family$data$tr[ind],exp(-bs$family$data$h[ind]),lwd=2,col=i)lines(bs$family$data$tr[ind],exp(-bs$family$data$h[ind] +      2*bs$family$data$q[ind]^.5),lty=2,col=i) ## upper cilines(bs$family$data$tr[ind],exp(-bs$family$data$h[ind] -      2*bs$family$data$q[ind]^.5),lty=2,col=i) ## lower cilines(bs$family$data$tr[ind],exp(-bs$family$data$km[ind]),col=i) ## KM}## Simple simulated known truth example...ph.weibull.sim <- function(eta,gamma=1,h0=.01,t1=100) {   lambda <- h0*exp(eta)   n <- length(eta)  U <- runif(n)  t <- (-log(U)/lambda)^(1/gamma)  d <- as.numeric(t <= t1)  t[!d] <- t1  list(t=t,d=d)}n <- 500;set.seed(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)^10f3 <- function(x) 0*xf <- f0(x0) + f1(x1)  + f2(x2)g <- (f-mean(f))/5surv <- ph.weibull.sim(g)surv$x0 <- x0;surv$x1 <- x1;surv$x2 <- x2;surv$x3 <- x3b <- gam(t~s(x0)+s(x1)+s(x2,k=15)+s(x3),family=cox.ph,weights=d,data=surv)plot(b,pages=1)## Another one, including a violation of proportional hazards for## effect of x2...set.seed(2)h <- exp((f0(x0)+f1(x1)+f2(x2)-10)/5)t <- rexp(n,h);d <- as.numeric(t<20)## first with no violation of PH in the simulation...b <- gam(t~s(x0)+s(x1)+s(x2)+s(x3),family=cox.ph,weights=d)plot(b,pages=1)ph.resid(b) ## fine## Now violate PH for x2 in the simulation...ii <- t>1.5h1 <- exp((f0(x0)+f1(x1)+3*f2(x2)-10)/5)t[ii] <- 1.5 + rexp(sum(ii),h1[ii]);d <- as.numeric(t<20)b <- gam(t~s(x0)+s(x1)+s(x2)+s(x3),family=cox.ph,weights=d)plot(b,pages=1)ph.resid(b) ## The checking plot picks up the problem in s(x2) ## conditional logistic regression models are often estimated using the ## cox proportional hazards partial likelihood with a strata for each## case-control group. A dummy vector of times is created (all equal). ## The following compares to 'clogit' for a simple case. Note that## the gam log likelihood is not exact if there is more than one case## per stratum, corresponding to clogit's approximate method.library(survival);library(mgcv)infert$dumt <- rep(1,nrow(infert))mg <- gam(cbind(dumt,stratum) ~ spontaneous + induced, data=infert,          family=cox.ph,weights=case)ms <- clogit(case ~ spontaneous + induced + strata(stratum), data=infert,             method="approximate")summary(mg)$p.table[1:2,]; ms

Additive Cox proportional hazard models with time varying covariates

Description

Thecox.ph family only allows one set of covariate values per subject. If each subject has several time varying covariate measurements then it is still possible to fit a proportional hazards regression model, via an equivalent Poisson model. The recipe is provided by Whitehead (1980) and is equally valid in the smooth additive case. Its drawback is that the equivalent Poisson dataset can be quite large.

The trick is to generate an artificial Poisson observation for each subject in the risk set at each non-censored event time. The corresponding covariate values for each subject are whatever they are at the event time, while the Poisson response is zero for all subjects except those experiencing the event at that time (this corresponds to Peto's correction for ties). The linear predictor for the model must include an intercept for each event time (the cumulative sum of the exponential of these is the Breslow estimate of the baseline hazard).

Below is some example code employing this trick for thepbcseq data from thesurvival package. It usesbam for fitting with thediscrete=TRUE option for efficiency: there is some approximation involved in doing this, and the exact equivalent to what is done incox.ph is rather obtained by usinggam withmethod="REML" (taking many times the computational time for the example below). An alternative fits the model as a conditional logistic model using stratified Cox PH with event times as strata (see example). This would be identical in the unpenalized case, but smoothing parameter estimates can differ.

The functiontdpois in the example code uses crude piecewise constant interpolation for the covariates, in which the covariate value at an event time is taken to be whatever it was the previous time that it was measured. Obviously more sophisticated interpolation schemes might be preferable.

References

Whitehead (1980) Fitting Cox's regression model to survival data using GLIM. Applied Statistics 29(3):268-275

Examples

require(mgcv);require(survival)## First define functions for producing Poisson model data frameapp <- function(x,t,to) {## wrapper to approx for calling from apply...  y <- if (sum(!is.na(x))<1) rep(NA,length(to)) else       approx(t,x,to,method="constant",rule=2)$y  if (is.factor(x)) factor(levels(x)[y],levels=levels(x)) else y} ## apptdpois <- function(dat,event="z",et="futime",t="day",status="status1",                                                            ) {## dat is data frame. id is patient id; et is event time; t is## observation time; status is 1 for death 0 otherwise;## event is name for Poisson response.  if (event %in% names(dat)) warning("event name in use")  require(utils) ## for progress bar  te <- sort(unique(dat[[et]][dat[[status]]==1])) ## event times  sid <- unique(dat[[id]])  inter <- interactive()  if (inter) prg <- txtProgressBar(min = 0, max = length(sid), initial = 0,         char = "=",width = NA, title="Progress", style = 3)  ## create dataframe for poisson model data  dat[[event]] <- 0; start <- 1  dap <- dat[rep(1:length(sid),length(te)),]  for (i in 1:length(sid)) { ## work through patients    di <- dat[dat[[id]]==sid[i],] ## ith patient's data    tr <- te[te <= di[[et]][1]] ## times required for this patient    ## Now do the interpolation of covariates to event times...    um <- data.frame(lapply(X=di,FUN=app,t=di[[t]],to=tr))    ## Mark the actual event...    if (um[[et]][1]==max(tr)&&um[[status]][1]==1) um[[event]][nrow(um)] <- 1     um[[et]] <- tr ## reset time to relevant event times    dap[start:(start-1+nrow(um)),] <- um ## copy to dap    start <- start + nrow(um)    if (inter) setTxtProgressBar(prg, i)  }  if (inter) close(prg)  dap[1:(start-1),]} ## tdpois## The following typically takes a minute or less...## Convert pbcseq to equivalent Poisson form...pbcseq$status1 <- as.numeric(pbcseq$status==2) ## death indicatorpb <- tdpois(pbcseq) ## conversionpb$tf <- factor(pb$futime) ## add factor for event time## Fit Poisson model...b <- bam(z ~ tf - 1 + sex + trt + s(sqrt(protime)) + s(platelet)+ s(age)+s(bili)+s(albumin), family=poisson,data=pb,discrete=TRUE,nthreads=2)pb$dumt <- rep(1,nrow(pb)) ## dummy time## Fit as conditional logistic... b1 <- gam(cbind(dumt,tf) ~ sex + trt + s(sqrt(protime)) + s(platelet)+ s(age) + s(bili) + s(albumin),family=cox.ph,data=pb,weights=z)par(mfrow=c(2,3))plot(b,scale=0)## compute residuals...chaz <- tapply(fitted(b),pb$id,sum) ## cum haz by subjectd <- tapply(pb$z,pb$id,sum) ## censoring indicatormrsd <- d - chaz ## Martingaledrsd <- sign(mrsd)*sqrt(-2*(mrsd + d*log(chaz))) ## deviance## plot survivor function and s.e. band for subject 25te <- sort(unique(pb$futime)) ## event timesdi <- pbcseq[pbcseq$id==25,] ## data for subject 25pd <- data.frame(lapply(X=di,FUN=app,t=di$day,to=te)) ## interpolate to tepd$tf <- factor(te)X <- predict(b,newdata=pd,type="lpmatrix")eta <- drop(X%*%coef(b)); H <- cumsum(exp(eta))J <- apply(exp(eta)*X,2,cumsum)se <- diag(J%*%vcov(b)%*%t(J))^.5plot(stepfun(te,c(1,exp(-H))),do.points=FALSE,ylim=c(0.7,1),     ylab="S(t)",xlab="t (days)",main="",lwd=2)lines(stepfun(te,c(1,exp(-H+se))),do.points=FALSE)lines(stepfun(te,c(1,exp(-H-se))),do.points=FALSE)rug(pbcseq$day[pbcseq$id==25]) ## measurement times

GAM censored Poisson family

Description

Family for use withgam orbam, implementing regression for censoredPoisson data. Observations may be left, interval or right censored or uncensored.

Usage

cpois(link="log")

Arguments

Details

If the family is used with a vector response, then it is assumed that there is no censoring, and a regular Poisson regression results. If there is censoring then the response should be supplied as a two column matrix. The first column is always numeric. Entries in the second column are as follows.

• If an entry is identical to the corresponding first column entry, then it is an uncensored observation.

• If an entry is numeric and different to the first column entry then there is interval censoring. The first column entry is the lower interval limit and the second column entry is the upper interval limit (both should be non-integer).y is only known to be between these limits.

• If the second column entry is-Inf then the observation is left censored at the value of the entry (non-integer) in the first column. It is only known thaty is below the first column value.

• If the second column entry isInf then the observation is right censored at the value (non-integer) of the entry in the first column. It is only known thaty is above the first column value.

Any mixture of censored and uncensored data is allowed, but be aware that data consisting only of right and/or left censored data contain very little information. It is strongly recommended to use non-integer values for censoring limits, to avoid any possibility of ambiguity. For example ify is known to be 3 or above, set the lower censoring limit to 2.5, or any other value in(2,3).

Value

An object of classextended.family.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Examples

library(mgcv)set.seed(6); n <- 2000dat <- gamSim(1,n=n,dist="poisson",scale=.1) ## simulate Poi data## illustration that cpois an poisson give same results if there## is no censoring...b0 <- gam(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+             s(x3,bs="cr"),family=poisson,data=dat,method="REML")plot(b0,pages=1,scheme=2)b1 <- gam(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+            s(x3,bs="cr"),family=cpois,data=dat) plot(b1,pages=1,scheme=2)b0;b1## Now censor some observations...dat1 <- datm <- 300y <- matrix(dat$y,n,ncol=2) ## response matrixii <- sample(n,3*m) ## censor thesefor (i in 1:m) { ## right, left, interval...  j <- ii[i]; if (y[j,1] > 4.5) y[j,] <- c(4.5,Inf)  j <- ii[m+i]; if (y[j,1] < 2.5) y[j,] <- c(2.5,-Inf)  j <- ii[2*m+i];if (y[j,1] > 1.5 & y[j,1]< 5.5) y[j,] <- c(1.5,5.5)}n - sum(y[,1]==y[,2]) ## number of censored obsdat1$y <- y## now fit model with censoring...b2 <- gam(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+          s(x3,bs="cr"),family=cpois,data=dat1) plot(b2,pages=1,scheme=2);b2

Obtaining derivative w.r.t. linear predictor

Description

INTERNAL function. Distribution families provide derivatives of the deviance and link w.r.t.mu = inv_link(eta).This routine converts these to the required derivatives of the deviance w.r.t. eta, the linear predictor.

Usage

dDeta(y, mu, wt, theta, fam, deriv = 0)

Arguments

Value

A list of derivatives.

Author(s)

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

Stable evaluation of difference between normal c.d.f.s

Description

Evaluates the difference between twoN(0,1) cumulative distribution functions avoiding cancellation error.

Usage

dpnorm(x0,x1,log.p=FALSE)

Arguments

Details

Equivalent topnorm(x1)-pnorm(x0), but stable whenx0 andx1 values are very close, or in the upper tail of the standard normal.

Author(s)

Simon N. Woodsimon.wood@r-project.org

Examples

require(mgcv)x <- seq(-10,10,length=10000)eps <- 1e-10y0 <- pnorm(x+eps)-pnorm(x) ## cancellation proney1 <- dpnorm(x,x+eps)       ## stable## illustrate stable computation in black, and## cancellation prone in red...par(mfrow=c(1,2),mar=c(4,4,1,1))plot(log(y1),log(y0),type="l")lines(log(y1[x>0]),log(y0[x>0]),col=2)plot(x,log(y1),type="l")lines(x,log(y0),col=2)

Stable evaluation of difference between Poisson c.d.f.s

Description

Evaluates the difference between two Poisson cumulative distribution functions avoiding cancellation error.

Usage

dppois(y0,y1,mu,log.p=TRUE)

Arguments

Details

Equivalent tolog(ppois(y1,mu)-ppois(y0,mu)), but stable wheny0 andy1 are highly improbable and in the same tail of the distribution (avoids cancellation tolog(0)).

Author(s)

Simon N. Woodsimon.wood@r-project.org

Examples

require(mgcv)n <- 10mu <- rep(2,n)y0 <- c(0:5,10,23,0,2)y1 <- c(1:6,11,1000,100,40)dppois(y0,y1,mu)log(ppois(y1,mu)-ppois(y0,mu))

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 elementsTRUE if the corresponding node is too far from data andFALSE otherwise. Basically a service routine forvis.gam andplot.gam.

Usage

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

Arguments

Details

Linear scalings of the axes are first determined so that the grid defined by the nodes ing1 andg2 lies exactly in the unit square (i.e. on [0,1] by [0,1]). These scalings are applied tog1,g2,d1 andd2. The minimum Euclidean distance from each node to a datum is then determined and if it is greater thandist the corresponding entry in the returned array is set toTRUE (otherwise toFALSE). The distance calculations are performed in compiled code for speed without storage overheads.

Value

A logical array withTRUE indicating a node in the grid defined byg1,g2 that is ‘too far’ from any datum.

Author(s)

Simon N. Woodsimon.wood@r-project.org

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 forgamm. Extracts the estimated covariance matrix of the data from anlme 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=NULL,start.level=1)extract.lme.cov2(b,data=NULL,start.level=1)

Arguments

.

Details

The random effects, correlation structure and variance structure usedfor 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!).

Theextract.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 thecovariance matrix is diagonal, then only the leading diagonal is returned; ifit can be written as a block diagonal matrix (under some permutation of theoriginal data) then a list of matrices defining the non-zero blocks isreturned along with an index indicating which row of the original data eachrow/column of the block diagonal matrix relates to. The block sizes are defined bythe coarsest level of grouping in the random effect structure.

gamm usesextract.lme.cov2.

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

Value

Forextract.lme.cov an estimated covariance matrix.

Forextract.lme.cov2 a list containing the estimated covariance matrixand an indexing array. The covariance matrix is stored as the elements on theleading diagonal, a list of the matrices defining a block diagonal matrix, ora full matrix if the previous two options are not possible.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Forlme 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 usinglme see:

Wood, S.N. (2006) Low rank scale invariant tensor product smooths forGeneralized Additive Mixed Models. Biometrics 62(4):1025-1036

or

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

See Also

gamm,formXtViX

Examples

## see also ?formXtViX for use of extract.lme.cov2require(mgcv)library(nlme)data(Rail)b <- lme(travel~1,Rail,~1|Rail)extract.lme.cov(b)extract.lme.cov2(b)

Factor smooth interactions in GAMs

Description

The interaction of one or more factors with a smooth effect, produces a separate smooth for each factor level. These smooths can have different smoothing parameters, or all have the same smoothing parameter. There are several vays to set them up.

Factorby variables.

If theby variables for a smooth (specified usings,te,ti ort2) is a factor, then a separate smooth is produced for each factor level. If the factor is ordered, then no smooth is produced for its first level: this is useful for setting up models which have a reference level smooth and then difference to reference smooths for each factor level except the first (which is the reference). Giving the smooth anid forces the same smoothing parameter to be used for all levels of the factor. For examples(x,by=fac,id=1) would produce a separate smooth ofx for each level offac, with each smooth having the same smoothing parameter. Seegam.models for more.

Sum to zero smooth interactions

bs="sz" These factor smooth interactions are specified usings(...,bs="sz"). There may be several factors supplied, and a smooth is produced for each combination of factor levels. The smooths are constructed to exclude the ‘main effect’ smooth, or the effects of individual smooths produced for lower order combinations of factor levels. For example, with a single factor, the smooths for the different factor levels are so constrained that the sum over all factor levels of equivalent spline coefficients are all zero. This allows the meaningful and identifiable construction of models with a main effect smooth plus smooths for the difference between each factor level and the main effect. Such a construction is often more natural than theby variable with ordered factors construction. Seesmooth.construct.sz.smooth.spec.

Random wiggly curves

bs="fs" This approach produces a smooth curve for each level of a single factor, treating the curves as entirely random. This means that in principle a model can be constructed with a main effect plus factor level smooth deviations from that effect. However the model is not forced to make the main effect do as much of the work as possible, in the way that the"sz" approach does. This approach can be very efficient withgamm as it exploits the nested estimation available inlme. Seesmooth.construct.fs.smooth.spec.

Author(s)

Simon N. Woodsimon.wood@r-project.org with input from Matteo Fasiolo.

See Also

smooth.construct.fs.smooth.spec,smooth.construct.sz.smooth.spec

Examples

library(mgcv)set.seed(0)## simulate data...f0 <- function(x) 2 * sin(pi * x)f1 <- function(x,a=2,b=-1) exp(a * x)+bf2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 *             (10 * x)^3 * (1 - x)^10n <- 500;nf <- 25fac <- sample(1:nf,n,replace=TRUE)x0 <- runif(n);x1 <- runif(n);x2 <- runif(n)a <- rnorm(nf)*.2 + 2;b <- rnorm(nf)*.5f <- f0(x0) + f1(x1,a[fac],b[fac]) + f2(x2)fac <- factor(fac)y <- f + rnorm(n)*2## so response depends on global smooths of x0 and ## x2, and a smooth of x1 for each level of fac.## fit model...bm <- gamm(y~s(x0)+ s(x1,fac,bs="fs",k=5)+s(x2,k=20))plot(bm$gam,pages=1)summary(bm$gam)bd <- bam(y~s(x0)+ s(x1) + s(x1,fac,bs="sz",k=5)+s(x2,k=20),discrete=TRUE)plot(bd,pages=1)summary(bd)## Could also use...## b <- gam(y~s(x0)+ s(x1,fac,bs="fs",k=5)+s(x2,k=20),method="ML")## ... but its slower (increasingly so with increasing nf)## b <- gam(y~s(x0)+ t2(x1,fac,bs=c("tp","re"),k=5,full=TRUE)+##        s(x2,k=20),method="ML"))## ... is exactly equivalent.

Distribution families in mgcv

Description

As well as the standard families (of classfamily) documented infamily (see alsoglm) which can be used with functionsgam,bam andgamm,mgcv also supplies some extra families, most of which are currently only usable withgam, although some can also be used withbam. These are described here.

Details

The following families (classfamily) are in the exponential family given the value of a single parameter. They are usable with all modelling functions.

• Tweedie An exponential family distribution for which the variance of the response is given by the mean response to the powerp.p is in (1,2) and must be supplied. Alternatively, seetw to estimatep (gam/bam only).

• negbin The negative binomial. Alternatively seenb to estimate thetheta parameter of the negative binomial (gam/bam only).

The following families (classextended.family) are for regression type models dependent on a single linear predictor, and with a log likelihoodwhich is a sum of independent terms, each corresponding to a single response observation. Usable withgam, with smoothing parameter estimation by"NCV","REML" or"ML" (the latter does not integrate the unpenalized and parameteric effects out of the marginal likelihood optimized for the smoothing parameters). Also usable withbam.

• betar for proportions data on (0,1) when the binomial is not appropriate.

• cnorm censored normal distribution, for log normal accelerated failure time models, Tobit regression and rounded data, for example.

• bcg Box-Cox transformed censored Gaussian.

• clog censored logistic distribution, for accelerated failure time models.

• cpois censored Poisson distribution.

• nb for negative binomial data when thetheta parameter is to be estimated.

• ocat for ordered categorical data.

• scat scaled t for heavy tailed data that would otherwise be modelled as Gaussian.

• tw for Tweedie distributed data, when the power parameter relating the variance to the mean is to be estimated.

• ziP for zero inflated Poisson data, when the zero inflation rate depends simply on the Poisson mean.

The above families of classfamily andextended.family can be combined to model data where different response observations come from different distributions. For example, when modelling the combination of presence-absence and abundance data,binomial andnb families might be used.

• gfam creates a 'grouped family' (or 'family group') from a list of families. The response is supplied as a two column matrix, the first containing the response observations, and the second an index of the family to which each observation relates.

The following families (classgeneral.family) implement more general model classes. Usable only withgam and only with REML or NCV smoothing parameter estimation.

• cox.ph the Cox Proportional Hazards model for survival data (no NCV).

• gammals a gamma location-scale model, where the mean and standard deviation are modelled with separate linear predictors.

• gaulss a Gaussian location-scale model where the mean and the standard deviation are both modelled using smooth linear predictors.

• gevlss a generalized extreme value (GEV) model where the location, scale and shape parameters are each modelled using a linear predictor.

• gumbls a Gumbel location-scale model (2 linear predictors).

• multinom: multinomial logistic regression, for unordered categorical responses.

• mvn: multivariate normal additive models (no NCV).

• shash Sinh-arcsinh location scale and shape model family (4 linear predicors).

• twlss Tweedie location scale and variance power model family (3 linear predicors). Can only be fitted using EFS method.

• ziplss a ‘two-stage’ zero inflated Poisson model, in which 'potential-presence' is modelled with one linear predictor, and Poisson mean abundancegiven potential presence is modelled with a second linear predictor.

Author(s)

Simon N. Wood (s.wood@r-project.org), Natalya Pya, Matteo Fasiolo and Chris Shen.

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Modify families for use in GAM fitting and checking

Description

Generalized Additive Model fitting by ‘outer’ iteration,requires extra derivatives of the variance and link functions to be added to family objects. The first 3 functions add what is needed. Model checking canbe aided by adding quantile and random deviate generating functions to the family. The final two functions do this.

Usage

fix.family.link(fam)fix.family.var(fam)fix.family.ls(fam)fix.family.qf(fam)fix.family.rd(fam)

Arguments

Details

Consider the first 3 function first.

Outer iteration GAM estimation requires derivatives of the GCV, UBRE/gAIC,GACV, REML or ML score, which are obtained by finding the derivatives of the modelcoefficients w.r.t. the log smoothing parameters, using the implicit function theorem. The expressions for the derivatives require the second and third derivatives of the link w.r.t. the mean (and the 4th derivatives if Fisher scoring is not used). Also required are the firstand second derivatives of the variance function w.r.t. the mean (plus the third derivative if Fisher scoring is not used). Finally REML or ML estimation of smoothing parametersrequires the log saturated likelihood and its first two derivatives w.r.t. the scale parameter.These functions add functions evaluating these quantities to a family.

If the family already has functionsdvar,d2var,d3var,d2link,d3link,d4link and for RE/MLls, then these functions simply return the family unmodified: this allows non-standard linksto be used withgam when using outer iteration (performanceiteration operates with unmodified families). Note that if you only need Fisher scoring thend4link andd3var can be dummy, as they are ignored. Similalryls is only needed for RE/ML.

Thedvar function is a function of a mean vector,mu, and returnsa vector of corresponding first derivatives of the family variancefunction. Thed2link function is also a function of a vector of meanvalues,mu: it returns a vector of second derivatives of the link,evaluated atmu. Higher derivatives are defined similarly.

If modifying your own family, note that you can often get away with supplyingonly advar andd2var, function if your family only requires links that occur inone of the standard families.

The second two functions are useful for investigating the distribution of residuals and are used byqq.gam. If possible the functions add quantile (qf) or random deviate (rd) generating functions to the family. If a family already hasqf orrd functions then it is left unmodified.qf functions are only available for some families, and for quasi families neither type of function is available.

Value

A family object with extra component functionsdvar,d2var,d2link,d3link,d4link,ls, and possiblyqf andrd, depending on which functions are called.fix.family.var also adds a variablescale set tonegative to indicate that family has a free scale parameter.

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

gam.fit3,qq.gam

Detect linear dependencies of one matrix on another

Description

Identifies columns of a matrixX2 which are linearlydependent on columns of a matrixX1. Primarily of use in setting up identifiability constraints for nested GAMs.

Usage

fixDependence(X1,X2,tol=.Machine$double.eps^.5,rank.def=0,strict=FALSE)

Arguments

Details

The algorithm uses a simple approach based on QR decomposition: seeWood (2017, section 5.6.3) for details.

Value

A vector of the columns ofX2 which are linearly dependent oncolumns ofX1 (or which need to be deleted to acheive independence and full rank ifstrict==FALSE).NULL if the two matrices are independent.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

Examples

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

Form component of GAMM covariance matrix

Description

This is a service routine forgamm. Given,V, an estimated covariance matrix obtained usingextract.lme.cov2 thisroutine forms a matrix square root of X^TV^{-1}X as efficiently as possible, giventhe structure ofV (usually sparse).

Usage

formXtViX(V,X)

Arguments

Details

The covariance matrix returned byextract.lme.cov2 maybe in a packed and re-ordered format, since it is usually sparse. Hence aspecial service routine is required to form the required products involvingthis matrix.

Value

A matrix, R such thatcrossprod(R) gives X^TV^{-1}X.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Forlme 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 for estimation usinglme see:

Wood, S.N. (2006) Low rank scale invariant tensor product smooths forGeneralized Additive Mixed Models. Biometrics 62(4):1025-1036

See Also

gamm,extract.lme.cov2

Examples

require(mgcv)library(nlme)data(ergoStool)b <- lme(effort ~ Type, data=ergoStool, random=~1|Subject)V1 <- extract.lme.cov(b, ergoStool)V2 <- extract.lme.cov2(b, ergoStool)X <- model.matrix(b, data=ergoStool)crossprod(formXtViX(V2, X))t(X)

GAM formula

Description

Description ofgam formula (see Details), and how to extract it from a fittedgam object.

Usage

## S3 method for class 'gam'formula(x,...)

Arguments

Details

gam will accept a formula or, with some families, a list of formulae. Othermgcv modelling functions will not accept a list. The list form provides a mechanism for specifying several linear predictors, and allows these to share terms: see below.

The formulae supplied togam are exactly like those supplied toglm except that smooth terms,s,te,ti andt2 canbe added to the right hand side (and. is not supported ingam formulae).

Smooth terms are specified by expressions of the form:
s(x1,x2,...,k=12,fx=FALSE,bs="tp",by=z,id=1)
wherex1,x2, etc. are the covariates which the smoothis a function of, andk is the dimension of the basis used torepresent the smooth term. Ifk is notspecified then basis specific defaults are used. Note that these defaults areessentially arbitrary, and it is important to check that they are not so small that they cause oversmoothing (too large just slows down computation). Sometimes the modelling context suggests sensible values fork, but if notinformal checking is easy: seechoose.k andgam.check.

fx is used to indicate whether or not this term should be unpenalized, and therefore have a fixed number of degrees of freedom set byk (almost alwaysk-1).bs indicates the basis to use for the smooth:the built in options are described insmooth.terms, and user defined smooths can be added (seeuser.defined.smooth). Ifbs is not supplied then the default"tp" (tprs) basis is used.by can be used to specify a variable by whichthe smooth should be multiplied. For examplegam(y~s(x,by=z))would specify a model E(y) = f(x)z wheref(\cdot) is a smooth function. Theby option is particularly useful for models inwhich different functions of the same variable are required foreach level of a factor and for ‘varying coefficient models’: seegam.models.id is used to give smooths identities: smooths with the same identity havethe same basis, penalty and smoothing parameter (but different coefficients, so they are different functions).

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 covariatesx andz 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.te provides further information.ti terms are a variant designed to be used as interaction terms when the main effects (and any lower order interactions) are present.t2 produces tensor productsmooths that are the natural low rank analogue of smoothing spline anova models.

s,te,ti andt2 terms accept ansp argument of supplied smoothing parameters: positive values are taken as fixed values to be used, negative to indicate that the parameter should be estimated. Ifsp is supplied then it over-rides whatever is in thesp argument togam, if it is not supplied then it defaults to all negative, but does not over-ride thesp argument togam.

Formulae can involve nested or “overlapping” terms such as
y~s(x)+s(z)+s(x,z) ory~s(x,z)+s(z,v)
but nested models should really be set up usingti terms:seegam.side for further details and examples.

Smooth terms in agam formula will accept matrix arguments as covariates (and correspondingby variable), in which case a ‘summation convention’ is invoked. Consider the example ofs(X,Z,by=L) whereX,ZandL are n by m matrices. LetF be the n by m matrix that results from evaluating the smooth at the values inX andZ. Then the contribution to the linear predictor from the term will berowSums(F*L) (note the element-wise multiplication). This convention allows the linear predictor of the GAMto depend on (a discrete approximation to) any linear functional of a smooth: seelinear.functional.terms for more information and examples (including functional linear models/signal regression).

Note thatgam allows any term in the model formula to be penalized (possibly by multiple penalties), via theparaPen argument. Seegam.models for details and example code.

When several formulae are provided in a list, then they can be used to specify multiple linear predictors for families for which this makes sense (e.g.mvn). The first formula in the list must include a response variable, but later formulae need not (depending on the requirements of the family). Let the linear predictors be indexed, 1 to d where d is the number of linear predictors, and the indexing is in the order in which the formulae appear in the list. It is possible to supply extra formulae specifying that several linear predictors should share some terms. To do this a formula is supplied in which the response is replaced by numbers specifying the indices of the linear predictors which will shre the terms specified on the r.h.s. For example1+3~s(x)+z-1 specifies that linear predictors 1 and 3 will share the termss(x) andz (but we don't want an extra intercept, as this would usually be unidentifiable). Note that it is possible that a linear predictor only includes shared terms: it must still have its own formula, but the r.h.s. would simply be-1 (e.g.y ~ -1 or~ -1). Seemultinom for an example.

Value

Returns the model formula,x$formula. Provided so thatanova methodsprint an appropriate description of the model.

WARNING

Agam formula should not refer to variables using e.g.dat[["x"]].

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

gam

FELSPLINE test function

Description

Implements a finite area test function based on one proposed by Tim Ramsay (2002).

Usage

fs.test(x,y,r0=.1,r=.5,l=3,b=1,exclude=TRUE)fs.boundary(r0=.1,r=.5,l=3,n.theta=20)

Arguments

Details

The function details are not given in the source article: but this is prettyclose. The function is modified from Ramsay (2002), in that it bulges, ratherthan being flat: this makes a better test of the smoother.

Value

fs.test returns function evaluations, orNAs for pointsoutside the boundary.fs.boundary returns a list ofx,y pointsto be jointed up in order to define/draw the boundary.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Tim Ramsay (2002) "Spline smoothing over difficult regions" J.R.Statist. Soc. B 64(2):307-319

Examples

require(mgcv)## plot the function, and its boundary...fsb <- fs.boundary()m<-300;n<-150 xm <- seq(-1,4,length=m);yn<-seq(-1,1,length=n)xx <- rep(xm,n);yy<-rep(yn,rep(m,n))tru <- matrix(fs.test(xx,yy),m,n) ## truthimage(xm,yn,tru,col=heat.colors(100),xlab="x",ylab="y")lines(fsb$x,fsb$y,lwd=3)contour(xm,yn,tru,levels=seq(-5,5,by=.25),add=TRUE)

Generalized additive models with integrated smoothness estimation

Description

Fits a generalized additive model (GAM) todata, the term ‘GAM’ being taken to include any quadratically penalized GLM and a variety of other models estimated by a quadratically penalised likelihood type approach (seefamily.mgcv). The degree of smoothness of model terms is estimated as part offitting.gam can also fit any GLM subject to multiple quadratic penalties (including estimation of degree of penalization). Confidence/credible intervals are readilyavailable for any quantity predicted using a fitted model.

Smooth terms are represented using penalized regression splines (or similar smoothers)with smoothing parameters selected by GCV/UBRE/AIC/REML/NCV or by regression splines withfixed 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), and users can add smooths. Linear functionals of smooths can also be included in models.For an overview of the smooths available seesmooth.terms. For more on specifying models seegam.models,random.effects andlinear.functional.terms. For more on model selection seegam.selection. Do readgam.check andchoose.k.

See packagegam, for GAMs via the original Hastie and Tibshirani approach (see details for differences to this implementation).

For very large datasets seebam, for mixed GAM seegamm andrandom.effects.

Usage

gam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,    na.action,offset=NULL,method="GCV.Cp",    optimizer=c("outer","newton"),control=list(),scale=0,    select=FALSE,knots=NULL,sp=NULL,min.sp=NULL,H=NULL,gamma=1,    fit=TRUE,paraPen=NULL,G=NULL,in.out,drop.unused.levels=TRUE,    drop.intercept=NULL,nei=NULL,discrete=FALSE,...)

Arguments

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)\} = \alpha + f_1(x_{1i})+f_2(x_{2i})

where the (independent) response variablesy_i \sim {\rm Poi }, andf_1 andf_2 are smooth functions of covariatesx_1 andx_2. The log is an example of a link function. Note that to be identifiable the modelrequires constraints on the smooth functions. By default these are imposed automatically and require that the function sums to zero over the observed covariate values (the presence of a metricby variable is the only case which usually suppresses this).

If absolutely any smooth functions were allowed in model fitting then maximum likelihood estimation of such models would invariably result in complex over-fitting estimates off_1 andf_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 trade-off 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.

Themgcv implementation ofgam 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 inmgcv 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

whereD is the deviance,n the number of data,sthe scale parameter andDoF the effective degrees of freedom of the model. Notice that UBRE is effectivelyjust AIC rescaled, but is only used whens is known.

Alternatives are GACV,NCV or a Laplace approximation to REML. Thereis some evidence that the latter may actually be the most effective choice. The main computational challenge solved by themgcv package is to optimize the smoothness selection criteria efficiently and reliably.

Broadlygam 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. The linear identifiability constraints are also obtained at this point. The model is fit usinggam.fit3 or variants, which are modifications ofglm.fit. The GAM penalized likelihood maximization problem is solved by Penalized Iteratively Re-weighted Least Squares (P-IRLS) (see e.g. Wood 2000). Smoothing parameter selection is possible in one of three ways. (i)‘Performance iteration’ uses the fact that at each P-IRLS step a working penalized linear modelis estimated, and the smoothing parameter estimation can be performed for each such working model.Eventually, in most cases, both model parameter estimates and smoothing parameter estimates converge. This option is available inbam andgamm.(ii) Alternatively the P-IRLS scheme is iterated to convergence for each trial set of smoothing parameters,and GCV, UBRE or REML scores are only evaluated on convergence - optimization is then ‘outer’ to the P-IRLSloop: in this case the P-IRLS iteration has to be differentiated, tofacilitate optimization, andgam.fit3 or one of its variants is used. (iii) The extended Fellner-Schall algorithm of Wood and Fasiolo (2017) alternates estimation of model coefficients with simple updates of smoothing parameters, eventually approximately maximizing the marginal likelihood of the model (REML).gam uses the second method, outer iteration, by default.

Several alternative basis-penalty types are built in for representing modelsmooths, but alternatives can easily be added (seesmooth.terms for an overview andsmooth.construct for how to add smooth classes). The choice of the basis dimension (k in thes,te,ti andt2 terms) is something that should be considered carefully (the exact value is not critical, but it is important not to make it restrictively small, nor very large and computationally costly). The basis should be chosen to be larger than is believed to be necessary to approximate the smooth function concerned. The effective degrees of freedom for the smooth will then be controlled by the smoothing penalty on the term, and (usually) selected automatically (with an upper limit set byk-1 or occasionallyk). Of course thek should not be made too large, or computation will be slow (or inextreme cases there will be more coefficients to estimate than there are data).

Note thatgam assumes a very inclusive definition of what counts as a GAM: basically any penalized GLM can be used: to this endgam allows the non smooth model components to be penalized via argumentparaPen and allows the linear predictor to depend on general linear functionals of smooths, via the summation convention mechanism described inlinear.functional.terms.link{family.mgcv} details what is available beyond GLMs and the exponential family.

Details of the default underlying fitting methods are given in Wood (2011, 2004) and Wood, Pya and Saefken (2016). Some alternative methods are discussed in Wood (2000, 2017).

gam() is not a clone of Trevor Hastie's original (as supplied in S-PLUS or packagegam). The majordifferences are (i) that by default estimation of thedegree of smoothness of model terms is part of model fitting, (ii) aBayesian approach to variance estimation is employed that makes for easierconfidence interval calculation (with good coverage probabilities), (iii) that the modelcan depend on any (bounded) linear functional of smooth terms, (iv) the parametric part of the model can be penalized, (v) simple random effects can be incorporated, and (vi) the facilities for incorporating smooths of more than one variable aredifferent: specifically there are nolo smooths, but instead (a)sterms can have more than one argument, implying an isotropic smooth and (b)te,ti ort2 smooths areprovided as an effective means for modelling smooth interactions of anynumber of variables via scale invariant tensor product smooths. Splines on the sphere, Duchon splines and Gaussian Markov Random Fields are also available. (vii) Models beyond the exponential family are available. See packagegam, for GAMs via the original Hastie and Tibshirani approach.

Value

Iffit=FALSE the function returns a listG of items needed tofit a GAM, but doesn't actually fit it.

Otherwise the function returns an object of class"gam" as described ingamObject.

WARNINGS

The default basis dimensions used for smooth terms are essentially arbitrary, and it should be checked that they are not too small. Seechoose.k andgam.check.

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

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 binomialfamilies. The problem is that with e.g. log or logit links, mean value zero corresponds toan infinite range on the linear predictor scale.

Author(s)

Simon N. Woodsimon.wood@r-project.org

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

References

Key References on this implementation:

Wood, S.N. (2025) Generalized Additive Models. Annual Review ofStatistics and Its Applications 12:497-526doi:10.1146/annurev-statistics-112723-034249

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models (with discussion).Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36doi:10.1111/j.1467-9868.2010.00749.x

Wood, S.N. (2004) Stable and efficient multiple smoothing parameter estimation forgeneralized additive models. J. Amer. Statist. Ass. 99:673-686. [Defaultmethod for additive case by GCV (but no longer for generalized)]

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

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

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.doi:10.1201/9781315370279

Wood, S.N. and M. Fasiolo (2017) A generalized Fellner-Schall method for smoothingparameter optimization with application to Tweedie location, scale and shape models.Biometrics 73 (4), 1071-1081doi:10.1111/biom.12666

Wood S.N., F. Scheipl and J.J. Faraway (2013) Straightforward intermediate rank tensor product smoothingin mixed models. Statistics and Computing 23: 341-360.doi:10.1007/s11222-012-9314-z

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals for Generalized AdditiveModel Components. Scandinavian Journal of Statistics, 39(1), 53-74.doi:10.1111/j.1467-9469.2011.00760.x

Key Reference on GAMs and related models:

Wood, S. N. (2020) Inference and computation with generalizedadditive models and their extensions. Test 29(2): 307-339.doi:10.1007/s11749-020-00711-5

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

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

Wahba (1990) Spline Models of Observational Data. SIAM

Wood, S.N. (2000) Modelling and Smoothing Parameter Estimationwith Multiple Quadratic Penalties. J.R.Statist.Soc.B 62(2):413-428 [The originalmgcv paper, but no longer the default methods.]

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 viathe Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

Gu (2002) Smoothing Spline ANOVA Models, Springer.

McCullagh and Nelder (1989) Generalized Linear Models 2nd ed. Chapman & Hall.

O'Sullivan, Yandall and Raynor (1986) Automatic smoothing of regressionfunctions 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

https://webhomes.maths.ed.ac.uk/~swood34/

See Also

mgcv-package,gamObject,gam.models,smooth.terms,linear.functional.terms,s,tepredict.gam,plot.gam,summary.gam,gam.side,gam.selection,gam.controlgam.check,linear.functional.termsnegbin,magic,vis.gam

Examples

## see also examples in ?gam.models (e.g. 'by' variables, ## random effects and tricks for large binary datasets)library(mgcv)set.seed(2) ## simulate some data... dat <- gamSim(1,n=400,dist="normal",scale=2)b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)summary(b)plot(b,pages=1,residuals=TRUE)  ## show partial residualsplot(b,pages=1,seWithMean=TRUE) ## `with intercept' CIs## plot derivaives of smooths w.r.t. covariates - the## direct analogue of linear model regression coefficients## - also selecting a more colourful plotting scheme...plot(b,pages=1,scale=0,deriv=TRUE,scheme=2)## run some basic model checks, including checking## smoothing basis dimensions...gam.check(b)## same fit in two parts .....G <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),fit=FALSE,data=dat)b <- gam(G=G)print(b)## 2 part fit enabling manipulation of smoothing parameters...G <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),fit=FALSE,data=dat,sp=b$sp)G$lsp0 <- log(b$sp*10) ## provide log of required sp vecgam(G=G) ## it's smoother## change the smoothness selection method to REMLb0 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,method="REML")## use alternative plotting scheme, and way intervals include## smoothing parameter uncertainty...plot(b0,pages=1,scheme=1,unconditional=TRUE) ## Would a smooth interaction of x0 and x1 be better?## Use tensor product smooth of x0 and x1, basis ## dimension 49 (see ?te for details, also ?t2).bt <- gam(y~te(x0,x1,k=7)+s(x2)+s(x3),data=dat,          method="REML")plot(bt,pages=1) plot(bt,pages=1,scheme=2) ## alternative visualizationAIC(b0,bt) ## interaction worse than additive## Alternative: test for interaction with a smooth ANOVA ## decomposition (this time between x2 and x1)bt <- gam(y~s(x0)+s(x1)+s(x2)+s(x3)+ti(x1,x2,k=6),            data=dat,method="REML")summary(bt)## If it is believed that x0 and x1 are naturally on ## the same scale, and should be treated isotropically ## then could try...bs <- gam(y~s(x0,x1,k=40)+s(x2)+s(x3),data=dat,          method="REML")plot(bs,pages=1)AIC(b0,bt,bs) ## additive still better. ## Now do automatic terms selection as wellb1 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,       method="REML",select=TRUE)plot(b1,pages=1)## 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),data=dat)plot(bp,pages=1,scheme=1)## alternatively...bp <- gam(y~s(x0,sp=.01)+s(x1)+s(x2)+s(x3),data=dat)# 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),data=dat) print(b);print(bp)# same with REMLbp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),        min.sp=c(0.1,0.1,0,10),data=dat,method="REML") print(b0);print(bp)## now a GAM with 3df regression spline term & 2 penalized termsb0 <- gam(y~s(x0,k=4,fx=TRUE,bs="tp")+s(x1,k=12)+s(x2,k=15),data=dat)plot(b0,pages=1)## now simulate poisson data...set.seed(6)dat <- gamSim(1,n=2000,dist="poisson",scale=.1)## use "cr" basis to save time, with 2000 data...b2<-gam(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+        s(x3,bs="cr"),family=poisson,data=dat,method="REML")plot(b2,pages=1)## drop x3, but initialize sp's from previous fit, to ## save more time...b2a<-gam(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr"),         family=poisson,data=dat,method="REML",         in.out=list(sp=b2$sp[1:3],scale=1))par(mfrow=c(2,2))plot(b2a)par(mfrow=c(1,1))## similar example using GACV...dat <- gamSim(1,n=400,dist="poisson",scale=.25)b4<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,        data=dat,method="GACV.Cp",scale=-1)plot(b4,pages=1)## repeat using REML as in Wood 2011...b5<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,        data=dat,method="REML")plot(b5,pages=1) ## a binary example (see ?gam.models for large dataset version)...dat <- gamSim(1,n=400,dist="binary",scale=.33)lr.fit <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=binomial,              data=dat,method="REML")## plot model components with truth overlaid in redop <- par(mfrow=c(2,2))fn <- c("f0","f1","f2","f3");xn <- c("x0","x1","x2","x3")for (k in 1:4) {  plot(lr.fit,residuals=TRUE,select=k)  ff <- dat[[fn[k]]];xx <- dat[[xn[k]]]  ind <- sort.int(xx,index.return=TRUE)$ix  lines(xx[ind],(ff-mean(ff))[ind]*.33,col=2)}par(op)anova(lr.fit)lr.fit1 <- gam(y~s(x0)+s(x1)+s(x2),family=binomial,               data=dat,method="REML")lr.fit2 <- gam(y~s(x1)+s(x2),family=binomial,               data=dat,method="REML")AIC(lr.fit,lr.fit1,lr.fit2)## For a Gamma example, see ?summary.gam...## For inverse Gaussian, see ?rig## now 2D smoothing...eg <- gamSim(2,n=500,scale=.1)attach(eg)op <- par(mfrow=c(2,2),mar=c(4,4,1,1))contour(truth$x,truth$z,truth$f) ## contour truthb4 <- gam(y~s(x,z),data=data) ## fit modelfit1 <- matrix(predict.gam(b4,pr,se=FALSE),40,40)contour(truth$x,truth$z,fit1)   ## contour fitpersp(truth$x,truth$z,truth$f)    ## persp truthvis.gam(b4)                     ## persp fitdetach(eg)par(op)#################################################### largish dataset example with user defined knots##################################################par(mfrow=c(2,2))n <- 5000eg <- gamSim(2,n=n,scale=.5)attach(eg)ind<-sample(1:n,200,replace=FALSE)b5<-gam(y~s(x,z,k=40),data=data,        knots=list(x=data$x[ind],z=data$z[ind]))## various visualizationsvis.gam(b5,theta=30,phi=30)plot(b5)plot(b5,scheme=1,theta=50,phi=20)plot(b5,scheme=2)par(mfrow=c(1,1))## and a pure "knot based" spline of the same datab6<-gam(y~s(x,z,k=64),data=data,knots=list(x= rep((1:8-0.5)/8,8),        z=rep((1:8-0.5)/8,rep(8,8))))vis.gam(b6,color="heat",theta=30,phi=30)## varying the default large dataset behaviour via `xt'b7 <- gam(y~s(x,z,k=40,xt=list(max.knots=500,seed=2)),data=data)vis.gam(b7,theta=30,phi=30)detach(eg)

Some diagnostics for a fitted gam model

Description

Takes a fittedgam object produced bygam() and produces some diagnostic informationabout the fitting procedure and results. The default is to produce 4 residualplots, some information about the convergence of the smoothness selection optimization, and to run diagnostic tests of whether the basis dimension choises are adequate. Care should be taken in interpreting the results when applied togam objects returned bygamm.

Usage

gam.check(b, old.style=FALSE,          type=c("deviance","pearson","response"),          k.sample=5000,k.rep=200,          rep=0, level=.9, rl.col=2, rep.col="gray80", ...)

Arguments

Details

Checking a fittedgam is like checking a fittedglm, with two main differences. Firstly, the basis dimensions used for smooth terms need to be checked, to ensure that they are not so small that they force oversmoothing: the defaults are arbitrary.choose.k provides more detail, but the diagnostic tests described below and reported by this function may also help. Secondly, fitting may not always be as robust to violation of the distributional assumptions as would be the case for a regular GLM, so slightly more care may be needed here. In particular, the thoery of quasi-likelihood implies that if the mean variance relationship is OK for a GLM, then other departures from the assumed distribution are not problematic: GAMs can sometimes be more sensitive. For example, un-modelled overdispersion will typically lead to overfit, as the smoothness selection criterion tries to reduce the scale parameter to the one specified. Similarly, it is not clear how sensitive REML and ML smoothness selection will be to deviations from the assumed response dsistribution. For these reasons this routine uses an enhanced residual QQ plot.

This function plots 4 standard diagnostic plots, some smoothing parameter estimationconvergence information and the results of tests which may indicate if the smoothing basis dimensionfor a term is too low.

Usually the 4 plots are various residual plots. For the default optimization methods the convergence information is summarized in a readable way, but for other optimization methods, whatever is returned by way ofconvergence diagnostics is simply printed.

The test of whether the basis dimension for a smooth is adequate (Wood, 2017, section 5.9) is based on computing an estimate of the residual variance based on differencing residuals that are near neighbours according to the (numeric) covariates of the smooth. This estimate divided by the residual variance is thek-index reported. The further below 1 this is, the more likely it is that there is missed pattern left in the residuals. Thep-value is computed by simulation: the residuals are randomly re-shuffledk.rep times to obtain the null distribution of the differencing variance estimator, if there is no pattern in the residuals. For models fitted to more thank.sample data, the tests are based ofk.sample randomly sampled data. Low p-values may indicate that the basis dimension,k, has been set too low, especially if the reportededf is close to k', the maximum possible EDF for the term. Note the disconcerting fact that if the test statistic itself is based on random resampling and the null is true, then the associated p-values will of course vary widely from one replicate to the next. Currently smooths of factor variables are not supported and will give anNA p-value.

Doubling a suspectk and re-fitting is sensible: if the reportededf increases substantially then you may have been missing something in the first fit. Of course p-values can be low for reasons other than a too lowk. Seechoose.k for fuller discussion.

The QQ plot produced is usually created by a call toqq.gam, and plots deviance residuals against approximate theoretical quantilies of the deviance residual distribution, according to the fitted model. If this looks odd then investigate further usingqq.gam. Note that residuals for models fitted to binary data contain very little information useful for model checking (it is necessary to find some way of aggregating them first), so the QQ plot is unlikely to be useful in this case.

Take care when interpreting results from applying this function to a model fitted usinggamm. In this case the returnedgam object is based on the working model used for estimation, and will treat all the random effects as part of the error. This means that the residuals extracted from thegam object are not standardized for the family used or for the random effects or correlation structure. Usually it is necessary to produce your own residual checks based on consideration of the model structure you have used.

Value

A vector of reference quantiles for the residual distribution, if these can be computed.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

N.H. Augustin, E-A Sauleaub, S.N. Wood (2012) On quantile quantile plots for generalized linear models.Computational Statistics & Data Analysis. 56(8), 2404-3409.

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

See Also

choose.k,gam,magic

Examples

library(mgcv)set.seed(0)dat <- gamSim(1,n=200)b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)plot(b,pages=1)gam.check(b,pch=19,cex=.3)

Setting GAM fitting defaults

Description

This is an internal function of packagemgcv which allows control of the numerical options for fitting a GAM. Typically users will want to modify the defaults if model fitting fails toconverge, or if the warnings are generated which suggest a loss of numerical stability during fitting. To change the defaultchoise of fitting method, seegam argumentsmethod andoptimizer.

Usage

gam.control(nthreads=1,ncv.threads=1,irls.reg=0.0,epsilon = 1e-07,            maxit = 200,mgcv.tol=1e-7,mgcv.half=15, trace = FALSE,            rank.tol=.Machine$double.eps^0.5,nlm=list(),    optim=list(),newton=list(),    idLinksBases=TRUE,scalePenalty=TRUE,efs.lspmax=15,    efs.tol=.1,keepData=FALSE,scale.est="fletcher",    edge.correct=FALSE)

Arguments

Details

Outer iteration usingnewton is controlled by the listnewtonwith the following elements:conv.tol (default1e-6) is the relative convergence tolerance;maxNstep is the maximumlength allowed for an element of the Newton search direction (default 5);maxSstep is the maximum length allowed for an element of the steepestdescent direction (only used if Newton fails - default 2);maxHalf isthe maximum number of step halvings to permit before giving up (default 30).

If outer iteration usingnlm is used for fitting, then the control listnlm stores control arguments for calls to routinenlm. The list has the following named elements: (i)ndigit isthe number of significant digits in the GCV/UBRE score - by default this isworked out fromepsilon; (ii)gradtol is the tolerance used tojudge convergence of the gradient of the GCV/UBRE score to zero - by defaultset to10*epsilon; (iii)stepmax is the maximum allowable logsmoothing parameter step - defaults to 2; (iv)steptol is the minimumallowable step length - defaults to 1e-4; (v)iterlim is the maximumnumber of optimization steps allowed - defaults to 200; (vi)check.analyticals indicates whether the built in exact derivativecalculations should be checked numerically - defaults toFALSE. Any ofthese which are not supplied and named in the list are set to their defaultvalues.

Outer iteration usingoptim is controlled using listoptim, which currently has one element:factr which takesdefault value 1e7.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

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

See Also

gam,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 fitting methods used bygam opt for certainty of convergence over speed offit.bam opts for speed.

gam uses a nested iteration method (seegam.outer), inwhich each trial set of smoothing parameters proposed by an outer Newton algorithmrequire an inner Newton algorithm (penalized iteratively re-weighted least squares, PIRLS)to find the corresponding best fit model coefficients. Implicit differentiation is used tofind the derivatives of the coefficients with respect to log smoothing parameters, so that thederivatives of the smoothness selection criterion can be obtained, as required by the outer iteration.This approach is less expensive than it at first appears, since excellent starting values for the inneriteration are available as soon as the smoothing parameters start to converge. See Wood (2011) and Wood, Pya and Saefken (2016).

bam uses an alternative approach similar to ‘performance iteration’ or ‘PQL’. A single PIRLS iteration is run to find the model coefficients. At each step this requires the estimation of a working penalized linear model. Smoothing parameter selection is applied directly to this working model at each step (as if it were a Gaussian additive model). This approach is more straightforward to code and in principle less costly than the nested approach. However it is not guaranteed to converge, since the smoothness selection criterion is changing at each iteration. It is sometimes possible for the algorithm to cycle around a small set of smoothing parameter, coefficient combinations without ever converging.bam includes some checks to limit this behaviour, and the further checks in the algorithm used bybam(...,discrete=TRUE) actually guarantee convergence in some cases, but in general guarantees are not possible. See Wood, Goude and Shaw (2015) and Wood et al. (2017).

gam when used with ‘general’ families (such asmultinom orcox.ph) can also use a potentially faster scheme based on the extended Fellner-Schall method (Wood and Fasiolo, 2017). This also operates with a single iteration and is not guaranteed to converge, theoretically.

There are three things that you can try to speed up GAM fitting. (i) if you have large numbers of smoothing parameters in the generalized case, then try the"bfgs" method option ingam argumentoptimizer: this can be faster than the default. (ii) Try usingbam(iii) For large datasets it may be worth changingthe smoothing basis to usebs="cr" (sees for details)for 1-d smooths, and to usete smooths in place ofs smooths for smooths of more than one variable. This is becausethe default thin plate regression spline basis"tp" is costly to set upfor large datasets.

If you have convergence problems, it's worth noting that a GAM is just a (penalized)GLM and the IRLS scheme used to estimate GLMs is not guaranteed toconverge. Hence non convergence of a GAM may relate to a lack of stability inthe basic IRLS scheme. Therefore it is worth trying to establish whether the IRLS iterationsare capable of converging. To do this fit the problematic GAM with all smoothterms specified withfx=TRUE so that the smoothing parameters are allfixed 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 itslikely that your model is just too flexible for the IRLS process itself. Having triedincreasingmaxit ingam.control, there are several otherpossibilities for stabilizing the iteration. It is possible to try (i) setting lower bounds on thesmoothing parameters using themin.sp argument ofgam: this mayor may not change the model being fitted; (ii)reducing the flexibility of the model by reducing the basis dimensionsk in the specification ofs andte model terms: thisobviously changes the model being fitted somewhat.

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

Please report convergence problems, especially if you there is no obvious pathology in the data/model thatsuggests convergence should fail.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Key References on this implementation:

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models (with discussion).Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

Wood, S.N., Goude, Y. & Shaw S. (2015) Generalized additive models for large datasets. Journal of the Royal Statistical Society, Series C 64(1): 139-155.

Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association.

Wood, S.N. and M. Fasiolo (2017) A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models, Biometrics.

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

GAM P-IRLS estimation with GCV/UBRE smoothness estimation (deprecated)

Description

This is a deprecated internal function of packagemgcv. It is a modificationof the functionglm.fit, designed to be called fromgam when perfomance iteration is selected (not the default). The majormodification is that rather than solving a weighted least squares problem at each IRLS step, a weighted, penalized least squares problemis solved at each IRLS step with smoothing parameters associated with each penalty chosen by GCV or UBRE,using routinemagic. For further information on usage see code forgam. Some regularization of the IRLS weights is also permitted as a way of addressing identifiability related problems (seegam.control). Negative binomial parameter estimation issupported.

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

Usage

 gam.fit(G, start = NULL, etastart = NULL,         mustart = NULL, family = gaussian(),         control = gam.control(),gamma=1,        fixedSteps=(control$maxit+1),...)

Arguments

Value

A list of fit information.

Author(s)

Simon N. Woodsimon.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 viathe Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

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

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

See Also

gam.fit3,gam,magic

P-IRLS GAM estimation with GCV, UBRE/AIC or RE/ML derivative calculation

Description

Estimation of GAM smoothing parameters is most stable ifoptimization of the UBRE/AIC, GCV, GACV, REML or ML score is outer to the penalized iterativelyre-weighted least squares scheme used to estimate the model given smoothing parameters.

This routine estimates a GAM (any quadratically penalized GLM) given log smoothing paramaters, and evaluates derivatives of the smoothness selection scores of the model with respect to thelog smoothing parameters. Calculation of exact derivatives is generally fasterthan approximating them by finite differencing, as well as generally improvingthe reliability of GCV/UBRE/AIC/REML score minimization.

The approach is to run the P-IRLSto convergence, and only then to iterate for first and second derivatives.

Not normally called directly, but rather service routines forgam.

Usage

gam.fit3(x, y, sp, Eb ,UrS=list(),          weights = rep(1, nobs), start = NULL, etastart = NULL,          mustart = NULL, offset = rep(0, nobs), U1 = diag(ncol(x)),          Mp = -1, family = gaussian(), control = gam.control(),          intercept = TRUE,deriv=2,gamma=1,scale=1,         printWarn=TRUE,scoreType="REML",null.coef=rep(0,ncol(x)),         pearson.extra=0,dev.extra=0,n.true=-1,Sl=NULL,nei=NULL,...)

Arguments

Details

This routine is basicallyglm.fit with somemodifications to allow (i) for quadratic penalties on the log likelihood;(ii) derivatives of the model coefficients with respect tolog smoothing parameters to be obtained by use of the implicit function theorem and (iii) derivatives of the GAM GCV, UBRE/AIC, REML or ML scores to beevaluated at convergence.

In addition the routines apply step halving to any step that increases thepenalized deviance substantially.

The most costly parts of the calculations are performed by calls to compiled Ccode (which in turn calls LAPACK routines) in place of the compiled code thatwould usually perform least squares estimation on the working model in theIRLS iteration.

Estimation of smoothing parameters by optimizing GCV scores obtained atconvergence 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 modificationof the families as described infix.family.link.

Author(s)

Simon N. Woodsimon.wood@r-project.org

The routine has been modified fromglm.fit in R 2.0.1, writtenby the R core (seeglm.fit for further credits).

References

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

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

See Also

gam,magic

Post-processing output of gam.fit5

Description

INTERNAL function for post-processing the output ofgam.fit5.

Usage

gam.fit5.post.proc(object, Sl, L, lsp0, S, off, gamma)

Arguments

Value

A list containing:

• R: unpivoted Choleski of estimated expected hessian of log-likelihood.

• Vb: the Bayesian covariance matrix of the model parameters.

• Ve: "frequentist" alternative toVb.

• Vc: corrected covariance matrix.

• F: matrix of effective degrees of freedom (EDF).

• edf:diag(F).

• edf2:diag(2F-FF).

Author(s)

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

Simple posterior simulation with gam fits

Description

GAM coefficients can be simulated directly from the Gaussian approximation to the posterior for the coefficients, or using a simple Metropolis Hastings sampler. See alsoginla.

Usage

gam.mh(b,ns=10000,burn=1000,t.df=40,rw.scale=.25,thin=1)

Arguments

Details

Posterior simulation is particularly useful for making inferences about non-linear functions of the model coefficients. Simulate random draws from the posterior, compute the function for each draw, and you have a draw from the posterior for the function. In many cases the Gaussian approximation to the posterior of the model coefficients is accurate, and samples generated from it can be treated as samples from the posterior for the coefficients. See example code below. This approach is computationally very efficient.

In other cases the Gaussian approximation can become poor. A typical example is in a spatial model with a log or logit link when there is a large area of observations containing only zeroes. In this case the linear predictor is poorly identified and the Gaussian approximation can become useless (an example is provided below). In that case it can sometimes be useful to simulate from the posterior using a Metropolis Hastings sampler. A simple approach alternates fixed proposals, based on the Gaussian approximation to the posterior, with random walk proposals, based on a shrunken version of the approximate posterior covariane matrix.gam.mh implements this. The fixed proposal often promotes rapid mixing, while the random walk component ensures that the chain does not become stuck in regions for which the fixed Gaussian proposal density is much lower than the posterior density.

The function reports the acceptance rate of the two types of step. If the random walk acceptance probability is higher than a quarter thenrw.step should probably be increased. Similarly if the acceptance rate is too low, it should be decreased. The random walk steps can be turned off altogether (see above), but it is important to check the chains for stuck sections if this is done.

Value

A list containing the retained simulated coefficients in matrixbs and two entries for the acceptance probabilities.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N. (2015) Core Statistics, Cambridge

Examples

library(mgcv)set.seed(3);n <- 400############################################## First example: simulated Tweedie model...############################################dat <- gamSim(1,n=n,dist="poisson",scale=.2)dat$y <- rTweedie(exp(dat$f),p=1.3,phi=.5) ## Tweedie responseb <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=tw(),          data=dat,method="REML")## simulate directly from Gaussian approximate posterior...br <- rmvn(1000,coef(b),vcov(b))## Alternatively use MH sampling...br <- gam.mh(b,thin=2,ns=2000,rw.scale=.15)$bs## If 'coda' installed, can check effective sample size## require(coda);effectiveSize(as.mcmc(br))## Now compare simulation results and Gaussian approximation for## smooth term confidence intervals...x <- seq(0,1,length=100)pd <- data.frame(x0=x,x1=x,x2=x,x3=x)X <- predict(b,newdata=pd,type="lpmatrix")par(mfrow=c(2,2))for(i in 1:4) {  plot(b,select=i,scale=0,scheme=1)  ii <- b$smooth[[i]]$first.para:b$smooth[[i]]$last.para  ff <- X[,ii]%*%t(br[,ii]) ## posterior curve sample  fq <- apply(ff,1,quantile,probs=c(.025,.16,.84,.975))  lines(x,fq[1,],col=2,lty=2);lines(x,fq[4,],col=2,lty=2)  lines(x,fq[2,],col=2);lines(x,fq[3,],col=2)}################################################################# Second example, where Gaussian approximation is a failure...###############################################################y <- c(rep(0, 89), 1, 0, 1, 0, 0, 1, rep(0, 13), 1, 0, 0, 1,        rep(0, 10), 1, 0, 0, 1, 1, 0, 1, rep(0,4), 1, rep(0,3),         1, rep(0, 3), 1, rep(0, 10), 1, rep(0, 4), 1, 0, 1, 0, 0,        rep(1, 4), 0, rep(1, 5), rep(0, 4), 1, 1, rep(0, 46))set.seed(3);x <- sort(c(0:10*5,rnorm(length(y)-11)*20+100))b <- gam(y ~ s(x, k = 15),method = 'REML', family = binomial)br <- gam.mh(b,thin=2,ns=2000,rw.scale=.4)$bsX <- model.matrix(b)par(mfrow=c(1,1))plot(x, y, col = rgb(0,0,0,0.25), ylim = c(0,1))ff <- X%*%t(br) ## posterior curve samplelinv <- b$family$linkinv## Get intervals for the curve on the response scale...fq <- linv(apply(ff,1,quantile,probs=c(.025,.16,.5,.84,.975)))lines(x,fq[1,],col=2,lty=2);lines(x,fq[5,],col=2,lty=2)lines(x,fq[2,],col=2);lines(x,fq[4,],col=2)lines(x,fq[3,],col=4)## Compare to the Gaussian posterior approximationfv <- predict(b,se=TRUE)lines(x,linv(fv$fit))lines(x,linv(fv$fit-2*fv$se.fit),lty=3)lines(x,linv(fv$fit+2*fv$se.fit),lty=3)## ... Notice the useless 95% CI (black dotted) based on the## Gaussian approximation!

Specifying generalized additive models

Description

This page is intended to provide some more information onhow to specify GAMs. A GAM is a GLM in which the linear predictor depends, in part, on a sum of smooth functions of predictors and (possibly) linear functionals of smooth functions of (possibly dummy) predictors.

Specifically lety_i denote an independent random variable with mean\mu_i and an exponential family distribution, or failing that a known mean variance relationship suitable for use of quasi-likelihood methods. Then the the linear predictor of a GAM has a structure something like

g(\mu_i) = {\bf X}_i{\beta} + f_1(x_{1i},x_{2i}) + f_2(x_{3i}) + L_i f_3(x_4) + \ldots

whereg is a known smooth monotonic ‘link’ function,{\bf X}_i\beta is the parametric part of the linear predictor, thex_j are predictor variables,thef_j are smooth functions andL_i is some linear functional off_3. There may of course be multiple linear functional terms, or none.

The key idea here is that thedependence of the response on the predictors can be represented as aparametric sub-model plus the sum of some (functionals of) smooth functions of one ormore of the predictor variables. Thus the model is quite flexiblerelative to strictly parametric linear or generalized linear models,but still has much more structure than the completely general modelthat says that the response is just some smooth function of all thecovariates.

Note one important point. In order for the model to be identifiablethe smooth functions usually have to be constrained to have zero mean (usuallytaken over the set of covariate values). The constraint is needed if the term involving the smooth includes a constant function in its span.gam always applies such constraints unless there is aby variable present, in which case an assessment is made of whether the constraint is needed or not (see below).

The following sections discuss specifying model structures forgam. Specification of the distribution and link function is done using thefamily argument togam and works in the same way as forglm. This page therefore concentrates on the model formula forgam.

Models with simple smooth terms

Consider the example model.

g(\mu_i) = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + f_1(x_{3i}) + f_2(x_{4i},x_{5i})

where the response variablesy_i has expectation\mu_iandg is a link function.

Thegam formula for this would be
y ~ x1 + x2 + s(x3) + s(x4,x5).
This would use the default basis for the smooths (a thin plateregression spline basis for each), with automatic selection of theeffective degrees of freedom for both smooths. The dimension of thesmoothing basis is given a default value as well (the dimension of thebasis sets an upper limit on the maximum possible degrees offreedom for the basis - the limit is typically one less than basisdimension). Full details of how to control smooths are given ins andte, and further discussion of basisdimension choice can be found inchoose.k. For the moment suppose that we would like to changethe basis of the first smooth to a cubic regression spline basis witha dimension of 20, while fixing the second term at 25 degrees offreedom. The appropriate formula would be:
y ~ x1 + x2 + s(x3,bs="cr",k=20) + s(x4,x5,k=26,fx=TRUE).

The above assumes thatx_{4} andx_5 are naturally on similar scales (e.g. they might be co-ordinates), so that isotropic smoothing is appropriate. If this assumption is false then tensor product smoothing might be better (seete).
y ~ x1 + x2 + s(x3) + te(x4,x5)
would generate a tensor product smooth ofx_{4} andx_5. By default this smooth would have basis dimension 25 and use cubic regression spline marginals. Varying the defaults is easy. For example
y ~ x1 + x2 + s(x3) + te(x4,x5,bs=c("cr","ps"),k=c(6,7))
specifies that the tensor product should use a rank 6 cubic regression spline marginaland a rank 7 P-spline marginal to create a smooth with basis dimension 42.

Nested terms/functional ANOVA

Sometimes it is interesting to specify smooth models with a main effects + interaction structure such as

E(y_i) = f_1(x_i) + f_2(z_i) + f_3(x_i,z_i)

or

E(y_i)=f_1(x_i) + f_2(z_i) + f_3(v_i)+ f_4(x_i,z_i) + f_5(z_i,v_i) + f_6(z_i,v_i) + f_7(x_i,z_i,v_i)

for example. Such models should be set up usingti terms in the model formula. For example:
y ~ ti(x) + ti(z) + ti(x,z), or
y ~ ti(x) + ti(z) + ti(v) + ti(x,z) + ti(x,v) + ti(z,v)+ti(x,z,v).
Theti terms produce interactions with the component main effects excluded appropriately. (There is in fact no need to useti terms for the main effects here,s terms could also be used.)

gam allows nesting (or ‘overlap’) ofte ands smooths, and automatically generates side conditions to make such models identifiable, but the resulting models are much less stable and interpretable than those constructed usingti terms.

‘by’ variables

by variables are the means for constructing ‘varying-coefficient models’ (geographic regression models) and for letting smooths ‘interact’ with factors or parametric terms. They are also the key to specifying general linear functionals of smooths.

Thes andte terms used to specify smooths accept an argumentby, which is a numeric or factor variable of the same dimension as the covariates of the smooth. If aby variable is numeric, then itsi^{th} element multiples thei^{th}row of the model matrix corresponding to the smooth term concerned.

Factor smooth interactions (see alsofactor.smooth.interaction).If aby variable is afactor then it generates an indicator vector for each level of the factor, unless it is anordered factor.In the non-ordered case, the model matrix for the smooth term is then replicated for each factor level,and each copy has its rows multiplied by the corresponding rows of itsindicator variable. The smoothness penalties are also duplicated for eachfactor level. In short a different smooth is generatedfor each factor level (theid argument tos andte can be used to force all such smooths to have the same smoothing parameter).orderedby variables are handled in the same way, except that no smooth is generated for the first level of the ordered factor (seeb3 example below). This is useful for setting up identifiable models when the same smooth occurs more than once in a model, with different factorby variables.

As an example, consider the model

E(y_i) = \beta_0+ f(x_i)z_i

wheref is a smooth function, andz_i is a numeric variable.The appropriate formula is:
y ~ s(x,by=z)
- theby argument ensures that the smooth function gets multiplied bycovariatez. Note that when using factor by variables, centering constraints are applied to the smooths,which usually means that the by variable should be included as a parametric term, as well.

The example code below also illustrates the use of factorby variables.

by variables may be supplied as numeric matrices as part of specifying general linear functional terms.

If aby variable is present and numeric (rather than a factor) then the corresponding smooth is only subjected to an identifiability constraint if (i) theby variable is a constant vector, or, (ii) for a matrixby variable,L, ifL%*%rep(1,ncol(L)) is constant or (iii) if a user defined smooth constructor supplies an identifiability constraint explicitly, and that constraint has an attibute"always.apply".

Linking smooths with ‘id’

It is sometimes desirable to insist that different smooth terms have the same degree of smoothness. This can be done by using theid argument tos orte terms. Smooths which share anid will have the same smoothing parameter. Really this only makes sense if the smooths use the same basis functions, and the default behaviour is to force this to happen: all smooths sharing anid have the same basis functions as the first smooth occurring with thatid. Note that if you want exactly the same function for each smooth, then this is best achieved by making use of the summation convention covered under ‘linear functional terms’.

As an example suppose thatE(y_i)\equiv\mu_i and

g(\mu_i) = f_1(x_{1i}) + f_2(x_{2i},x_{3i}) + f_3(x_{4i})

but thatf_1 andf_3 should have the same smoothing parameters (andx_2andx_3 are on different scales). Then thegam formula
y ~ s(x1,id=1) + te(x_2,x3) + s(x4,id=1)
would achieve the desired result.id can be numbers or character strings. Giving anid to a term with a factorby variable causes the smooths at each level of the factor to have the same smoothing parameter.

Smooth termids are not supported bygamm.

Linear functional terms

General linear functional terms have a long history in the spline literature including in the penalized GLM context (see e.g. Wahba 1990). Such terms encompass varying coefficient models/ geographic regression, functional GLMs (i.e. GLMs with functional predictors), GLASS models, etc, and allow smoothing with respect to aggregated covariate values, for example.

Such terms are implemented inmgcv using a simple ‘summation convention’ for smooth terms: If the covariates of a smooth are supplied as matrices, then summation of the evaluated smooth over the columns of the matrices is implied. Each covariate matrix and anyby variable matrix must be of the same dimension. Consider, for example the term
s(X,Z,by=L)
whereX,Z andL aren \times p matrices. Letf denote the thin plate regression spline specified. The resulting contibution to thei^{\rm th} element of the linear predictor is

\sum_{j=1}^p L_{ij}f(X_{ij},Z_{ij})

If noL is supplied then all its elements are taken as 1. In R code terms, letF denote then \times p matrix obtained by evaluating the smooth at the values inX andZ. Then the contribution of the term to the linear predictor isrowSums(L*F) (note that it's element by element multiplication here!).

The summation convention applies tote terms as well ass terms. More details and examples are provided inlinear.functional.terms.

Random effects

Random effects can be added togam models usings(...,bs="re") terms (seesmooth.construct.re.smooth.spec), or theparaPen argument togam covered below. Seegam.vcomp,random.effects andsmooth.construct.re.smooth.spec for further details. An alternative is to use the approach ofgamm.

Penalizing the parametric terms

In case the ability to add smooth classes, smooth identities,by variables and the summation convention are still not sufficient to implement exactly the penalized GLM that you require,gam also allows you to penalize the parametric terms in the model formula. This is mostly useful in allowing one or more matrix terms to be included in the formula, along with a sequence of quadratic penalty matrices for each.

Suppose that you have set up a model matrix\bf X, and want to penalize the corresponding coefficients,\betawith two penalties\beta^T {\bf S}_1 \beta and\beta^T {\bf S}_2 \beta. Then something like the following would be appropriate:
gam(y ~ X - 1,paraPen=list(X=list(S1,S2)))
TheparaPen argument should be a list with elements having names corresponding to the terms being penalized. Each element ofparaPen is itself a list, with optional elementsL,rank andsp: all other elements must be penalty matrices. If present,rank is a vector giving the rank of each penalty matrix (if absent this is determined numerically).L is a matrix that maps underlying log smoothing parameters to the log smoothing parameters that actually multiply the individual quadratic penalties: taken as the identity if not supplied.sp is a vector of (underlying) smoothing parameter values: positive values are taken as fixed, negative to signal that the smoothing parameter should be estimated. Taken as all negative if not supplied.

An obvious application ofparaPen is to incorporate random effects, and an example of this is provided below. In thiscase the supplied penalty matrices will be (generalized) inverse covariance matrices for the random effects — i.e. precision matrices. The final estimate of the covariance matrix corresponding to one of these penalties is given by the (generalized) inverse of the penalty matrix multiplied by the estimated scale parameter and divided by the estimated smoothing parameter for the penalty. For example, if you use an identity matrix to penalize some coefficients that are to be viewed as i.i.d. Gaussian random effects, then their estimated variance will be the estimated scale parameter divided by the estimate of the smoothing parameter, for this penalty. See the ‘rail’ example below.

P-values for penalized parametric terms should be treated with caution. If you must have them, then use the optionfreq=TRUE inanova.gam andsummary.gam, which will tend to give reasonable results for random effects implemented this way, but not for terms with a rank defficient penalty (or penalties with a wide eigen-spectrum).

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wahba (1990) Spline Models of Observational Data SIAM.

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

Examples

require(mgcv)set.seed(10)## simulate date from y = f(x2)*x1 + errordat <- gamSim(3,n=400)b<-gam(y ~ s(x2,by=x1),data=dat)plot(b,pages=1)summary(b)## Factor `by' variable example (with a spurious covariate x0)## simulate data...dat <- gamSim(4)## fit model...b <- gam(y ~ fac+s(x2,by=fac)+s(x0),data=dat)plot(b,pages=1)summary(b)## note that the preceding fit is the same as....b1<-gam(y ~ s(x2,by=as.numeric(fac==1))+s(x2,by=as.numeric(fac==2))+            s(x2,by=as.numeric(fac==3))+s(x0)-1,data=dat)## ... the `-1' is because the intercept is confounded with the ## *uncentred* smooths here.plot(b1,pages=1)summary(b1)## repeat forcing all s(x2) terms to have the same smoothing param## (not a very good idea for these data!)b2 <- gam(y ~ fac+s(x2,by=fac,id=1)+s(x0),data=dat)plot(b2,pages=1)summary(b2)## now repeat with a single reference level smooth, and ## two `difference' smooths...dat$fac <- ordered(dat$fac)b3 <- gam(y ~ fac+s(x2)+s(x2,by=fac)+s(x0),data=dat,method="REML")plot(b3,pages=1)summary(b3)rm(dat)## An example of a simple random effects term implemented via ## penalization of the parametric part of the model...dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth## Now add some random effects to the simulation. Response is ## grouped into one of 20 groups by `fac' and each groups has a## random effect added....fac <- as.factor(sample(1:20,400,replace=TRUE))dat$X <- model.matrix(~fac-1)b <- rnorm(20)*.5dat$y <- dat$y + dat$X%*%b## now fit appropriate random effect model...PP <- list(X=list(rank=20,diag(20)))rm <- gam(y~ X+s(x0)+s(x1)+s(x2)+s(x3),data=dat,paraPen=PP)plot(rm,pages=1)## Get estimated random effects standard deviation...sig.b <- sqrt(rm$sig2/rm$sp[1]);sig.b ## a much simpler approach uses "re" terms...rm1 <- gam(y ~ s(fac,bs="re")+s(x0)+s(x1)+s(x2)+s(x3),data=dat,method="ML")gam.vcomp(rm1)## Simple comparison with lme, using Rail data. ## See ?random.effects for a simpler methodrequire(nlme)b0 <- lme(travel~1,data=Rail,~1|Rail,method="ML") Z <- model.matrix(~Rail-1,data=Rail,     contrasts.arg=list(Rail="contr.treatment"))b <- gam(travel~Z,data=Rail,paraPen=list(Z=list(diag(6))),method="ML")b0 (b$reml.scale/b$sp)^.5 ## `gam' ML estimate of Rail sdb$reml.scale^.5         ## `gam' ML estimate of residual sdb0 <- lme(travel~1,data=Rail,~1|Rail,method="REML") Z <- model.matrix(~Rail-1,data=Rail,     contrasts.arg=list(Rail="contr.treatment"))b <- gam(travel~Z,data=Rail,paraPen=list(Z=list(diag(6))),method="REML")b0 (b$reml.scale/b$sp)^.5 ## `gam' REML estimate of Rail sdb$reml.scale^.5         ## `gam' REML estimate of residual sd################################################################## Approximate large dataset logistic regression for rare events## based on subsampling the zeroes, and adding an offset to## approximately allow for this.## Doing the same thing, but upweighting the sampled zeroes## leads to problems with smoothness selection, and CIs.################################################################n <- 50000  ## simulate n data dat <- gamSim(1,n=n,dist="binary",scale=.33)p <- binomial()$linkinv(dat$f-6) ## make 1's raredat$y <- rbinom(p,1,p)      ## re-simulate rare response## Now sample all the 1's but only proportion S of the 0'sS <- 0.02                   ## sampling fraction of zeroesdat <- dat[dat$y==1 | runif(n) < S,] ## sampling## Create offset based on total sampling fractiondat$s <- rep(log(nrow(dat)/n),nrow(dat))lr.fit <- gam(y~s(x0,bs="cr")+s(x1,bs="cr")+s(x2,bs="cr")+s(x3,bs="cr")+              offset(s),family=binomial,data=dat,method="REML")## plot model components with truth overlaid in redop <- par(mfrow=c(2,2))fn <- c("f0","f1","f2","f3");xn <- c("x0","x1","x2","x3")for (k in 1:4) {       plot(lr.fit,select=k,scale=0)       ff <- dat[[fn[k]]];xx <- dat[[xn[k]]]       ind <- sort.int(xx,index.return=TRUE)$ix       lines(xx[ind],(ff-mean(ff))[ind]*.33,col=2)}par(op)rm(dat)## A Gamma example, by modify `gamSim' output... dat <- gamSim(1,n=400,dist="normal",scale=1)dat$f <- dat$f/4 ## true linear predictor Ey <- exp(dat$f);scale <- .5 ## mean and GLM scale parameter## Note that `shape' and `scale' in `rgamma' are almost## opposite terminology to that used with GLM/GAM...dat$y <- rgamma(Ey*0,shape=1/scale,scale=Ey*scale)bg <- gam(y~ s(x0)+ s(x1)+s(x2)+s(x3),family=Gamma(link=log),          data=dat,method="REML")plot(bg,pages=1,scheme=1)

Minimize GCV or UBRE score of a GAM using ‘outer’ iteration

Description

Estimation of GAM smoothing parameters is most stable ifoptimization of the smoothness selection score (GCV, GACV, UBRE/AIC, REML, ML etc) is outer to the penalized iterativelyre-weighted least squares scheme used to estimate the model given smoothing parameters.

This routine optimizes a smoothness selection score in this way. Basically the score is evaluated for each trial set of smoothing parameters byestimating the GAM for those smoothing parameters. The score is minimizedw.r.t. the parameters numerically, usingnewton (default),bfgs,optim ornlm. Exact(first and second) derivatives of the score can be used by fitting withgam.fit3. Thisimproves efficiency and reliability relative to relying on finitedifference derivatives.

Not normally called directly, but rather a service routine forgam.

Usage

gam.outer(lsp,fscale,family,control,method,optimizer,          criterion,scale,gamma,G,start=NULL,nei=NULL,...)

Arguments

Details

See Wood (2008) for full details on ‘outer iteration’.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

See Also

gam.fit3,gam,magic

Finding stable orthogonal re-parameterization of the square root penalty.

Description

INTERNAL function for finding an orthogonal re-parameterization which avoids "dominant machine zero leakage" betweencomponents of the square root penalty.

Usage

gam.reparam(rS, lsp, deriv)

Arguments

Value

A list containing

• S: the total penalty matrix similarity transformed for stability.

• rS: the component square roots, transformed in the same way.

• Qs: the orthogonal transformation matrixS = t(Qs)%*%S0%*%Qs, whereS0 is theuntransformed total penalty implied bysp andrS on input.

• det: log|S|.

• det1: dlog|S|/dlog(sp) ifderiv >0.

• det2: hessian of log|S| wrt log(sp) ifderiv>1.

Author(s)

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

Scale parameter estimation in GAMs

Description

Scale parameter estimation ingam depends on the type offamily. For extended families then the RE/ML estimate is used. For conventional exponential families, estimated by the default outer iteration, the scale estimator can be controlled using argumentscale.est ingam.control. The options are"fletcher" (default),"pearson" or"deviance". The Pearson estimator is the (weighted) sum of squares of the pearson residuals, divided by the effective residual degrees of freedom. The Fletcher (2012) estimator is an improved version of the Pearson estimator.The deviance estimator simply substitutes deviance residuals for Pearson residuals.

Usually the Pearson estimator is recommended for GLMs, since it is asymptotically unbiased. However, it can also be unstable at finite sample sizes, if a few Pearson residuals are very large. For example, a very low Poisson mean with a non zero count can give a huge Pearson residual, even though the deviance residual is much more modest. The Fletcher (2012) estimator is designed to reduce these problems.

For performance iteration the Pearson estimator is always used.

gamm uses the estimate of the scale parameter from the underlying call tolme.bam uses the REML estimator if the method is"fREML". Otherwise the estimator is a Pearson estimator.

Author(s)

Simon N. Woodsimon.wood@r-project.org with help from Mark Bravington and David Peel

References

Fletcher, David J. (2012) Estimating overdispersion when fitting a generalized linear model to sparse data. Biometrika 99(1), 230-237.

See Also

gam.control

Generalized Additive Model Selection

Description

This page is intended to provide some more information onhow to select GAMs. In particular, it gives a brief overview of smoothness selection, and then discusses how this can be extended to select inclusion/exclusion of terms. Hypothesis testing approaches to the latter problem are also discussed.

Smoothness selection criteria

Given a model structure specified by a gam model formula,gam() attempts to find the appropriate smoothness for each applicable model term using prediction error criteria or likelihood based methods. The prediction error criteria used are Generalized (Approximate) Cross Validation (GCV or GACV) when the scale parameter is unknown or an Un-Biased Risk Estimator (UBRE) when it is known. UBRE is essentiallyscaled AIC (Generalized case) or Mallows' Cp (additive model case). GCV and UBRE are covered in Craven and Wahba (1979) and Wahba (1990). Alternatively REML of maximum likelihood (ML) may be used for smoothness selection, by viewing the smooth components as random effects (in this case the variance component for each smooth random effect will be given by the scale parameter divided by the smoothing parameter — for smooths with multiple penalties, there will be multiple variance components). Themethod argument togam selects the smoothness selection criterion.

Automatic smoothness selection is unlikely to be successful with few data, particularly withmultiple terms to be selected. In addition GCV and UBRE/AIC score can occasionallydisplay local minima that can trap the minimisation algorithms. GCV/UBRE/AICscores become constant with changing smoothing parameters at very low or veryhigh smoothing parameters, and on occasion these ‘flat’ regions can beseparated from regions of lower score by a small ‘lip’. This seems to be themost common form of local minimum, but is usually avoidable by avoidingextreme smoothing parameters as starting values in optimization, and byavoiding big jumps in smoothing parameters while optimizing. Never the less, if you aresuspicious of smoothing parameter estimates, try changing fit method (seegam argumentsmethod andoptimizer) and see if the estimates change, or try changingsome or all of the smoothing parameters ‘manually’ (argumentsp ofgam,orsp arguments tos orte).

REML and ML are less prone to local minima than the other criteria, and may therefore be preferable.

Automatic term selection

Unmodified smoothness selection by GCV, AIC, REML etc. will not usually remove a smooth from a model. This is because most smoothing penalties view some spaceof (non-zero) functions as ‘completely smooth’ and once a term is penalized heavily enough that it is in this space, further penalization does not change it.

However it is straightforward to modify smooths so that under heavy penalization they are penalized to the zero function and thereby ‘selected out’ of the model. There are two approaches.

The first approach is to modify the smoothing penalty with an additional shrinkage term. Smooth classescs.smooth andtprs.smooth (specified by"cs" and"ts" respectively) have smoothness penalties which include a smallshrinkage component, so that for large enough smoothing parameters the smooth becomes identically zero. This allows automatic smoothing parameter selectionmethods to effectively remove the term from the model altogether. Theshrinkage component of the penalty is set at a level that usually makes negligablecontribution to the penalization of the model, only becoming effective whenthe term is effectively ‘completely smooth’ according to the conventional penalty.

The second approach leaves the original smoothing penalty unchanged, but constructs an additional penalty for each smooth, which penalizes only functions in the null space of the original penalty (the ‘completely smooth’ functions). Hence, if all the smoothing parameters for a term tend to infinity, the term will be selected out of the model. This latter approach is more expensive computationally, but has the advantage that it can be applied automatically to any smooth term. Theselect argument togam turns on this method.

In fact, as implemented, both approaches operate by eigen-decomposiong the original penalty matrix. A new penalty is created on the null space: it is the matrix with the same eigenvectors as theoriginal penalty, but with the originally postive egienvalues set to zero, and the originally zero eigenvalues set to something positive. The first approach just addes a multiple of this penalty to the original penalty, where the multiple is chosen so that the new penalty can not dominate the original. The second approach treats the new penalty as an extra penalty, with its own smoothing parameter.

Of course, as with all model selection methods, some care must be take to ensure that the automatic selection is sensible, and a decision about the effective degrees of freedom at which to declare a term ‘negligible’ has to be made.

Interactive term selection

In general the most logically consistent method to use for deciding whichterms to include in the model is to compare GCV/UBRE/ML scores for models with and without the term (REML scores should not be used to compare models with different fixed effects structures). When UBRE is the smoothness selection method this willgive the same result as comparing byAIC (the AIC in this caseuses the model EDF in place of the usual model DF). Similarly, comparison viaGCV score and via AIC seldom yields different answers. Note that the negativebinomial with estimatedtheta parameter is a special case: the GCVscore is not informative, because of thetheta estimationscheme used. 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 replacement by parametric terms are smooth terms with estimated degrees of freedom close to their minimum possible.

Candidates for removal can also be identified by reference to the approximate p-values provided bysummary.gam, and by looking at theextent to which the confidence band for an estimated term includes the zero function. It is perfectly possible to perform backwards selection using p-valuesin the usual way: that is by sequentially dropping the single term with the highest non-significant p-value from the model and re-fitting, until all terms are significant. This suffers from the same problems as stepwise procedures for anyGLM/LM, with the additional caveat that the p-values are only approximate. If adopting this approach, it is probably best to use ML smoothness selection.

Note that GCV and UBRE are not appropriate for comparing models using different families: in that case AIC should be used.

Caveats/platitudes

Formal model selection methods are only appropriate for selecting between reasonable models.If formal model selection is attempted starting from a model that simply doesn't fit the data, then it is unlikely to provide meaningful results.

The more thought is given to appropriate model structure up front, the more successful modelselection is likely to be. Simply starting with a hugely flexible model with ‘everything in’ and hoping that automatic selection will find the right structure is not often successful.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Marra, G. and S.N. Wood (2011) Practical variable selection for generalized additive models.Computational Statistics and Data Analysis 55,2372-2387.

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. (2003) Thin plate regression splines. J.R.Statist.Soc.B 65(1):95-114

Wood, S.N. (2008) Fast stable direct fitting and smoothness selection forgeneralized additive models. J.R.Statist. Soc. B 70(3):495-518

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

See Also

gam,step.gam

Examples

## an example of automatic model selection via null space penalizationlibrary(mgcv)set.seed(3);n<-200dat <- gamSim(1,n=n,scale=.15,dist="poisson") ## simulate datadat$x4 <- runif(n, 0, 1);dat$x5 <- runif(n, 0, 1) ## spuriousb<-gam(y~s(x0)+s(x1)+s(x2)+s(x3)+s(x4)+s(x5),data=dat,         family=poisson,select=TRUE,method="REML")summary(b)plot(b,pages=1)

Identifiability side conditions for a GAM

Description

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

The method identifies nested and repeated variables by their names, butnumerically evaluates which constraints need to be imposed. Constraints are alwaysapplied to smooths of more variables in preference to smooths of fewervariables. The numerical approach allows appropriate constraints to beapplied to models constructed using any smooths, including user defined smooths.

Usage

gam.side(sm,Xp,tol=.Machine$double.eps^.5,with.pen=TRUE)

Arguments

Details

Models such asy~s(x)+s(z)+s(x,z) can be estimated bygam, but require identifiability constraints to be applied, tomake them identifiable. This routine does this, effectively setting redundant parametersto zero. When the redundancy is between smooths of lower and higher numbersof variables, the constraint is always applied to the smooth of the highernumber of variables.

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

Xp is used to check whether there is a constant term in the model (or columns that can be linearly combined to give a constant). This is because centred smooths can appear independent, when they would be dependent if there is a constant in the model, so dependence testing needs to take account of this.

Value

A list of smooths, with model matrices and penalty matrices adjustedto automatically impose the required constraints. Any smooth that has beenmodified will have an attribute"del.index", listing the columns of itsmodel matrix that were deleted. This index is used in the creation ofprediction matrices for the term.

WARNINGS

Much better statistical stability will be obtained by using models likey~s(x)+s(z)+ti(x,z) ory~ti(x)+ti(z)+ti(x,z) rather thany~s(x)+s(z)+s(x,z), since the former are designed not to require further constraint.

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

ti,gam.models

Examples

## The first two examples here iluustrate models that cause## gam.side to impose constraints, but both are a bad way ## of estimating such models. The 3rd example is the right## way.... set.seed(7)require(mgcv)dat <- gamSim(n=400,scale=2) ## simulate data## estimate model with redundant smooth interaction (bad idea).b<-gam(y~s(x0)+s(x1)+s(x0,x1)+s(x2),data=dat)plot(b,pages=1)## Simulate data with real interation...dat <- gamSim(2,n=500,scale=.1)old.par<-par(mfrow=c(2,2))## a fully nested tensor product example (bad idea)b <- gam(y~s(x,bs="cr",k=6)+s(z,bs="cr",k=6)+te(x,z,k=6),       data=dat$data)plot(b)old.par<-par(mfrow=c(2,2))## A fully nested tensor product example, done properly,## so that gam.side is not needed to ensure identifiability.## ti terms are designed to produce interaction smooths## suitable for adding to main effects (we could also have## used s(x) and s(z) without a problem, but not s(z,x) ## or te(z,x)).b <- gam(y ~ ti(x,k=6) + ti(z,k=6) + ti(x,z,k=6),       data=dat$data)plot(b)par(old.par)rm(dat)

Report gam smoothness estimates as variance components

Description

GAMs can be viewed as mixed models, where the smoothing parameters are related to variance components. This routine extracts the estimated variance components associated with each smooth term, and if possible returns confidence intervals on the standard deviation scale.

Usage

gam.vcomp(x,rescale=TRUE,conf.lev=.95)## S3 method for class 'gam.vcomp'print(x,...)

Arguments

Details

The (pseudo) inverse of the penalty matrix penalizing a term is proportional to the covariance matrix of the term's coefficients, when these are viewed as random. For single penalty smooths, it is possible to compute the variance component for the smooth (which multiplies the inverse penalty matrix to obtain the covariance matrix of the smooth's coefficients). Thisvariance component is given by the scale parameter divided by the smoothing parameter.

This routine computes such variance components, forgam models, and associated confidence intervals, if smoothing parameter estimation was likelihood based. Note that variance components are also returned for tensor product smooths, but that their interpretation is not so straightforward.

The routine is particularly useful for model fitted bygam in which random effects have been incorporated.

Value

Either a vector of variance components for each smooth term (as standarddeviations), or a list. A list will be returned if the smoothness selectionmethod used was REML or ML, or where a model has more smoothing parametersthan actually estimated.

If a list, it will contain some or all of the elements described depending onthe method of smoothness selection used, and whether some smoothing parameterswere not estimated.

Depending on the smoothness selection method, elementvc may contain amatrix or a vector. If REML or ML smoothness selection was used,vc will be a matrix whose first column gives standard deviations foreach term, while the subsequent columns give lower and upper confidencebounds, on the same scale. For models fitted with other smoothness selectionmethods, componentvc will be a vector of standard deviations.

Theall element is a vector of variance components forall the smoothing parameters (estimated + fixed or replicated).

Additionally, for REML or ML smoothness selection, the numerical rank of thecovariance matrix and the Hessian matrix are returns as elementsrankandrank.hess respectively.

Author(s)

Simon N. Woodsimon.wood@r-project.org modified by Gavin Simpson

References

Wood, S.N. (2008) Fast stable direct fitting and smoothnessselection for generalized additive models. Journal of the RoyalStatistical Society (B) 70(3):495-518

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

See Also

smooth.construct.re.smooth.spec

Examples

   set.seed(3)   require(mgcv)  ## simulate some data, consisting of a smooth truth + random effects  dat <- gamSim(1,n=400,dist="normal",scale=2)  a <- factor(sample(1:10,400,replace=TRUE))  b <- factor(sample(1:7,400,replace=TRUE))  Xa <- model.matrix(~a-1)    ## random main effects  Xb <-  model.matrix(~b-1)  Xab <- model.matrix(~a:b-1) ## random interaction  dat$y <- dat$y + Xa%*%rnorm(10)*.5 +            Xb%*%rnorm(7)*.3 + Xab%*%rnorm(70)*.7  dat$a <- a;dat$b <- b  ## Fit the model using "re" terms, and smoother linkage      mod <- gam(y~s(a,bs="re")+s(b,bs="re")+s(a,b,bs="re")+s(x0,id=1)+s(x1,id=1)+               s(x2,k=15)+s(x3),data=dat,method="ML")  gam.vcomp(mod)

Objective functions for GAM smoothing parameter estimation

Description

Estimation of GAM smoothing parameters is most stable ifoptimization of the UBRE/AIC or GCV score is outer to the penalized iterativelyre-weighted least squares scheme used to estimate the model given smoothing parameters. These functions evaluate the GCV/UBRE/AIC score of a GAM model, givensmoothing parameters, in a manner suitable for use byoptim ornlm.Not normally called directly, but rather service routines forgam.outer.

Usage

gam2objective(lsp,args,...)gam2derivative(lsp,args,...)

Arguments

Details

gam2objective andgam2derivative are functions suitablefor calling byoptim, to evaluate the GCV/UBRE/AIC score and itsderivatives w.r.t. log smoothing parameters.

gam4objective is an equivalent togam2objective, suitable foroptimization bynlm - derivatives of the GCV/UBRE/AIC function arecalculated and returned as attributes.

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

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

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

Wood, S.N. (2008) Fast stable direct fitting and smoothness selection for generalizedadditive models. J.R.Statist.Soc.B 70(3):495-518

See Also

gam.fit3,gam,magic

Fitted gam object

Description

A fitted GAM object returned by functiongam and of class"gam" inheriting from classes"glm" and"lm". Methodfunctionsanova,logLik,influence,plot,predict,print,residuals andsummary exist forthis 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, sothat the elements relating to the QR decomposition of the model matrix areabsent.

Value

Agam object has the following elements:

WARNINGS

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

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

A Key Reference on this implementation:

Wood, S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman& Hall/ CRC, Boca Raton, Florida

Key Reference on GAMs generally:

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

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

See Also

gam

Simulate example data for GAMs

Description

Function used to simulate data sets to illustrate the use ofgam andgamm. Mostly used in help files to keep down the length of the example code sections.

Usage

gamSim(eg=1,n=400,dist="normal",scale=2,verbose=TRUE)

Arguments

Details

See the source code for exactly what is simulated in each case.

1. Gu and Wahba 4 univariate term example.

2. A smooth function of 2 variables.

3. Example with continuous by variable.

4. Example with factor by variable.

5. An additive example plus a factor variable.

6. Additive + random effect.

7. As 1 but with correlated covariates.

Value

Depends oneg, but usually a dataframe, which may also contain some information on the underlying truth. Sometimes a list with more items, including a data frame for model fitting.See source code or helpfile examples where the function is used for further information.

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

gam,gamm

Examples

## see ?gam

Transform derivatives wrt mu to derivatives wrt linear predictor

Description

Mainly intended for internal use in specifying location scale models.Letg(mu) = lp, wherelp is the linear predictor, andg is the linkfunction. Assume that we have calculated the derivatives of the log-likelihood wrtmu.This function uses the chain rule to calculate the derivatives of the log-likelihood wrtlp. Seetrind.generator for array packing conventions.

Usage

gamlss.etamu(l1, l2, l3 = NULL, l4 = NULL, ig1, g2, g3 = NULL,  g4 = NULL, i2, i3 = NULL, i4 = NULL, deriv = 0)

Arguments

Value

A list where the arraysl1,l2,l3,l4 contain the derivatives (upto order four) of the log-likelihood wrt the linear predictor.

Author(s)

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

See Also

trind.generator

Calculating derivatives of log-likelihood wrt regression coefficients

Description

Mainly intended for internal use with location scale model families.Given the derivatives of the log-likelihood wrt the linear predictor, this function obtainsthe derivatives and Hessian wrt the regression coefficients and derivatives ofthe Hessian w.r.t. the smoothing parameters. For input derivative array packing conventions seetrind.generator.

Usage

gamlss.gH(X, jj, l1, l2, i2, l3 = 0, i3 = 0, l4 = 0, i4 = 0, d1b = 0,  d2b = 0, deriv = 0, fh = NULL, D = NULL,sandwich=FALSE)

Arguments

Value

A list containinglb - the grad vector w.r.t. coefs;lbb - the Hessian matrix w.r.t. coefs;d1H - either a list of the derivatives of the Hessian w.r.t. the smoothing parameters, or a single matrix whose columns are the leading diagonals of these dervative matrices;trHid2H - the trace of the inverse Hessian multiplied by the second derivative of the Hessian w.r.t. all combinations of smoothing parameters.

Author(s)

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

See Also

trind.generator

Generalized Additive Mixed Models

Description

Fits the specified generalized additive mixed model (GAMM) todata, by a call tolme in the normal errors identity link case, or by a call togammPQL (a modification ofglmmPQL from theMASS library) otherwise. In the latter case estimates are only approximately MLEs. The routine is typically slower thangam, and not quite as numerically robust.

To uselme4 in place ofnlme as the underlying fitting engine, seegamm4 from packagegamm4.

Smooths are specified as in a call togam 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 available forlme 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 interestis 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 ofgam andlme.

Usage

gamm(formula,random=NULL,correlation=NULL,family=gaussian(),data=list(),weights=NULL,subset=NULL,na.action,knots=NULL,control=list(niterEM=0,optimMethod="L-BFGS-B",returnObject=TRUE),niterPQL=20,verbosePQL=TRUE,method="ML",drop.unused.levels=TRUE,mustart=NULL, etastart=NULL,...)

Arguments

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 modelas 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 aslme. 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 inglmmPQL 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 thenlme library (using correlation structures beyond the strictly additive case amounts to using a GEE approach to fitting).

Some details of how GAMs are represented as mixed models and estimated usinglme orgammPQL ingamm can be found in Wood (2004 ,2006a,b). In additiongamm 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 ingam, although adaptivesmoothing bases are not available. Prediction from the returnedgam object is straightforward usingpredict.gam, but this will set the random effects to zero. If you want to predict with random effects set to their predicted values then you can adapt the prediction code given in the examples below.

In the event oflme convergence failures, considermodifyingoptions(mgcv.vc.logrange): reducing it helps to removeindefiniteness in the likelihood, if that is the problem, but too large areduction can force over or undersmoothing. SeenotExp2 for moreinformation on this option. Failing that, you can try increasing theniterEM option incontrol: this will perturb the starting valuesused in fitting, but usually to values with lower likelihood! Note that thisversion ofgamm works best with R 2.2.0 or above andnlme, 3.1-62 and above,since these use an improved optimizer.

Value

Returns a list with two items:

WARNINGS

gamm has a somewhat different argument list togam,gam arguments such asgamma supplied togamm willjust be ignored.

gamm performs poorly with binary data, since it uses PQL. It is better to usegam withs(...,bs="re") terms, orgamm4.

gamm assumes that you know what you are doing! For example, unlikeglmmPQL fromMASS it will return the completelme objectfrom the working model at convergence of the PQL iteration, including the 'log likelihood', even though this is not the likelihood of the fitted GAMM.

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

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

gamm is not as numerically stable asgam: anlme callwill occasionally fail. See details section for suggestions, or try the ‘gamm4’ package.

gamm is usually much slower thangam, and on some platforms you may need toincrease the memory available to R in order to use it with large data sets(seememory.limit).

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

Note that thegam object part of the returned object is not complete inthe sense of having all the elements defined ingamObject anddoes not inherit fromglm: hence e.g. multi-modelanova calls will not work.It is also based on the working model when PQL is used.

The parameterization used for the smoothing parameters ingamm, boundsthem above and below by an effective infinity and effective zero. SeenotExp2 for details of how to change this.

Linked smoothing parameters and adaptive smoothing are not supported.

Author(s)

Simon N. Woodsimon.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 Statisticswith S. Fourth edition. Springer.

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

Wood, S.N. (2004) Stable and efficient multiple smoothing parameter estimation forgeneralized 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. (2006a) Low rank scale invariant tensor product smooths forgeneralized additive mixed models. Biometrics 62(4):1025-1036

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

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

See Also

magic for an alternative for correlated data,te,s,predict.gam,plot.gam,summary.gam,negbin,vis.gam,pdTens,gamm4 (https://cran.r-project.org/package=gamm4)

Examples

library(mgcv)## simple examples using gamm as alternative to gamset.seed(0) dat <- gamSim(1,n=200,scale=2)b <- gamm(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)plot(b$gam,pages=1)summary(b$lme) # details of underlying lme fitsummary(b$gam) # gam style summary of fitted modelanova(b$gam) gam.check(b$gam) # simple checking plotsb <- gamm(y~te(x0,x1)+s(x2)+s(x3),data=dat) op <- par(mfrow=c(2,2))plot(b$gam)par(op)rm(dat)## Add a factor to the linear predictor, to be modelled as randomdat <- gamSim(6,n=200,scale=.2,dist="poisson")b2 <- gamm(y~s(x0)+s(x1)+s(x2),family=poisson,           data=dat,random=list(fac=~1))plot(b2$gam,pages=1)fac <- dat$facrm(dat)vis.gam(b2$gam)## In the generalized case the 'gam' object is based on the working## model used in the PQL fitting. Residuals for this are not## that useful on their own as the following illustrates...gam.check(b2$gam) ## But more useful residuals are easy to produce on a model## by model basis. For example...fv <- exp(fitted(b2$lme)) ## predicted values (including re)rsd <- (b2$gam$y - fv)/sqrt(fv) ## Pearson residuals (Poisson case)op <- par(mfrow=c(1,2))qqnorm(rsd);plot(fv^.5,rsd)par(op)## now an example with autocorrelated errors....n <- 200;sig <- 2x <- 0:(n-1)/(n-1)f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10e <- rnorm(n,0,sig)for (i in 2:n) e[i] <- 0.6*e[i-1] + e[i]y <- f + eop <- par(mfrow=c(2,2))## Fit model with AR1 residualsb <- gamm(y~s(x,k=20),correlation=corAR1())plot(b$gam);lines(x,f-mean(f),col=2)## Raw residuals still show correlation, of course...acf(residuals(b$gam),main="raw residual ACF")## But standardized are now fine...acf(residuals(b$lme,type="normalized"),main="standardized residual ACF")## compare with model without AR component...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 + eb <- 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 <- 10n<-n.g*10*4## simulate smooth part...dat <- gamSim(1,n=n,scale=2)f <- dat$f## 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 groupse<-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)}dat$y <- f + ra + rb + edat$fa <- fa;dat$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)summary(b$gam)vis.gam(b$gam)## Prediction from gam object, optionally adding ## in random effects. ## Extract random effects and make names more convenient...refa <- ranef(b$lme,level=5)rownames(refa) <- substr(rownames(refa),start=9,stop=20)refb <- ranef(b$lme,level=6)rownames(refb) <- substr(rownames(refb),start=9,stop=20)## make a prediction, with random effects zero...p0 <- predict(b$gam,data.frame(x0=.3,x1=.6,x2=.98,x3=.77))## add in effect for fa = "2" and fb="2/4"...p <- p0 + refa["2",1] + refb["2/4",1] ## and a "spatial" example...library(nlme);set.seed(1);n <- 100dat <- gamSim(2,n=n,scale=0) ## standard exampleattach(dat)old.par<-par(mfrow=c(2,2))contour(truth$x,truth$z,truth$f)  ## true functionf <- data$f                       ## true expected response## Now simulate correlated errors...cstr <- corGaus(.1,form = ~x+z)  cstr <- Initialize(cstr,data.frame(x=data$x,z=data$z))V <- corMatrix(cstr) ## correlation matrix for dataCv <- chol(V)e <- t(Cv) %*% rnorm(n)*0.05 # correlated errors## next add correlated simulated errors to expected valuesdata$y <- f + e ## ... to produce responseb<- gamm(y~s(x,z,k=50),correlation=corGaus(.1,form=~x+z),         data=data)plot(b$gam) # gamm fit accounting for correlation# overfits when correlation ignored.....  b1 <- gamm(y~s(x,z,k=50),data=data);plot(b1$gam) b2 <- gam(y~s(x,z,k=50),data=data);plot(b2)par(old.par)

Gamma location-scale model family

Description

Thegammals family implements gamma location scale additive models in which the log of the mean,\mu, and the log of the scale parameter,\phi (see details) can depend on additive smooth predictors. The parameterization is of the usual GLM type where the variance of the response is given by\phi\mu^2. Useable only withgam, the linear predictors are specified via a list of formulae.

Usage

gammals(link=list("identity","log"),b=-7)

Arguments

Details

Used withgam to fit gamma location - scale models parameterized in terms of the log mean and the log scale parameter (the response variance is the squared mean multiplied by the scale parameter). Note thatidentity links mean that the linear predictors give the log mean and log scale directly. By default thelog link for the scale parameter simply forces the log scale parameter to have a lower limit given by argumentb: if\eta is the linear predictor for the log scale parameter,\phi, then\log \phi = b + \log(1+e^{\eta-b}).

gam is called with a list containing 2 formulae, the first specifies the response on the left hand side and the structure of the linear predictor for the log mean on the right hand side. The second is one sided, specifying the linear predictor for the log scale on the right hand side.

The fitted values for this family will be a two column matrix. The first column is the mean (on original, not log, scale), and the second column is the log scale. Predictions usingpredict.gam will also produce 2 column matrices fortype"link" and"response". The first column is on the original data scale whentype="response" and on the log mean scale of the linear predictor whentype="link". The second column whentype="response" is again the log scale parameter, but is on the linear predictor whentype="link".

The null deviance reported for this family computed by setting the fitted values to the mean response, but using the model estimated scale.

Value

An object inheriting from classgeneral.family.

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Examples

library(mgcv)## simulate some dataf0 <- 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)^10f3 <- function(x) 0 * xn <- 400;set.seed(9)x0 <- runif(n);x1 <- runif(n);x2 <- runif(n);x3 <- runif(n);mu <- exp((f0(x0)+f2(x2))/5)th <- exp(f1(x1)/2-2)y <- rgamma(n,shape=1/th,scale=mu*th)b1 <- gam(list(y~s(x0)+s(x2),~s(x1)+s(x3)),family=gammals)plot(b1,pages=1)summary(b1)gam.check(b1)plot(mu,fitted(b1)[,1]);abline(0,1,col=2)plot(log(th),fitted(b1)[,2]);abline(0,1,col=2)

Gaussian location-scale model family

Description

Thegaulss family implements Gaussian location scale additive models in which the mean and the logb of the standard deviation (see details) can depend on additive smooth predictors. Useable only withgam, the linear predictors are specified via a list of formulae.

Usage

gaulss(link=list("identity","logb"),b=0.01)

Arguments

Details

Used withgam to fit Gaussian location - scale models.gam is called with a list containing 2 formulae, the first specifies the response on the left hand side and the structure of the linear predictor for the mean on the right hand side. The second is one sided, specifying the linear predictor for the standard deviation on the right hand side.

Link functions"identity","inverse","log" and"sqrt" are available for the mean. For the standard deviation only the"logb" link is implemented:\eta = \log(\sigma - b) and\sigma = b + \exp(\eta). This link is designed to avoid singularities in the likelihood caused by the standard deviation tending to zero. Note that internally the family is parameterized in terms of the\tau=\sigma^{-1}, i.e. the standard deviation of the precision, so the link and inverse link are coded to reflect this, however the relationships between the linear predictor and the standard deviation are as given above.

The fitted values for this family will be a two column matrix. The first column is the mean, and the second column is the inverse of the standard deviation. Predictions usingpredict.gam will also produce 2 column matrices fortype"link" and"response". The second column whentype="response" is again on the reciprocal standard deviation scale (i.e. the square root precision scale). The second column whentype="link" is\log(\sigma - b). Alsoplot.gam will plot smooths relating to\sigma on the\log(\sigma - b) scale (so high values correspond to high standard deviation and low values to low standard deviation). Similarly the smoothing penalties are applied on the (log) standard deviation scale, not the log precision scale.

The null deviance reported for this family is the sum of squares of the difference between the response and the mean response divided by the standard deviation of the response according to the model. The deviance is the sum of squares of residuals divided by model standard deviations.

Value

An object inheriting from classgeneral.family.

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Examples

library(mgcv);library(MASS)b <- gam(list(accel~s(times,k=20,bs="ad"),~s(times)),            data=mcycle,family=gaulss())summary(b) plot(b,pages=1,scale=0)

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 thenNULL is returned. The routine is useful forprocessing gam formulae. If the variable is a matrix then it is coerced to a numeric vector, by default.

Usage

 get.var(txt,data,vecMat=TRUE)

Arguments

Value

The evaluated variable orNULL. May be coerced to a numeric vector if it's a matrix.

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

gam

Examples

require(mgcv)y <- 1:4;dat<-data.frame(x=5:10)get.var("x",dat)get.var("y",dat)get.var("x==6",dat)dat <- list(X=matrix(1:6,3,2))get.var("X",dat)

Generalized Extreme Value location-scale model family

Description

Thegevlss family implements Generalized Extreme Value location scale additive models in which the location, scale and shape parameters depend on additive smooth predictors. Usable only withgam, the linear predictors are specified via a list of formulae.

Usage

gevlss(link=list("identity","identity","logit"))

Arguments

Details

Used withgam to fit Generalized Extreme Value distribution location scale and shape models. The p.d.f. is

t(y)^{\xi+1}e^{-t(y)}/\sigma

wheret(y) = \left \{1+\xi(y-\mu)/\sigma \right\}^{-1/\xi} if\xi \ne 0andt(y) = \exp\{-(y-\mu)/\sigma\} if\xi=0.

gam is called with a list containing 3 formulae: the first specifies the response on the left hand side and the structure of the linear predictor for the location parameter,\mu, on the right hand side. The second is one sided, specifying the linear predictor for the log scale parameter,\rho=\log(\sigma), on the right hand side. The third is one sided specifying the linear predictor for the shape parameter,\xi.

Link functions"identity" and"log" are available for the location (\mu) parameter. There is no choice of link for the log scale parameter (\rho = \log(\sigma)). The shape parameter (\xi) defaults to a modified logit link restricting its range to (-1,.5), the upper limit is required to ensure finite variance, while the lower limit ensures consistency of the MLE (Smith, 1985).

The fitted values for this family will be a three column matrix. The first column is the location parameter, the second column is the log scale parameter, the third column is the shape parameter.

This family does not produce a null deviance. Note that the distribution for\xi=0 is approximated by setting\xi to a small number.

The derivative system code for this family is mostly auto-generated, and the family is still somewhat experimental.

The GEV distribution is rather challenging numerically, and for small datasets or poorly fitting models improved numerical robustness may be obtained by using the extended Fellner-Schall method of Wood and Fasiolo (2017) for smoothing parameter estimation. See examples.

Value

An object inheriting from classgeneral.family.

References

Smith, R.L. (1985) Maximum likelihood estimation in a class ofnonregular cases. Biometrika 72(1):67-90

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Wood, S.N. and M. Fasiolo (2017) A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models. Biometrics 73(4): 1071-1081.doi:10.1111/biom.12666

Examples

library(mgcv)Fi.gev <- function(z,mu,sigma,xi) {## GEV inverse cdf.  xi[abs(xi)<1e-8] <- 1e-8 ## approximate xi=0, by small xi  x <- mu + ((-log(z))^-xi-1)*sigma/xi}## simulate test data...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)^10set.seed(1)n <- 500x0 <- runif(n);x1 <- runif(n);x2 <- runif(n)mu <- f2(x2)rho <- f0(x0)xi <- (f1(x1)-4)/9y <- Fi.gev(runif(n),mu,exp(rho),xi)dat <- data.frame(y,x0,x1,x2);pairs(dat)## fit model....b <- gam(list(y~s(x2),~s(x0),~s(x1)),family=gevlss,data=dat)## same fit using the extended Fellner-Schall method which## can provide improved numerical robustness... b <- gam(list(y~s(x2),~s(x0),~s(x1)),family=gevlss,data=dat,         optimizer="efs")## plot and look at residuals...plot(b,pages=1,scale=0)summary(b)par(mfrow=c(2,2))mu <- fitted(b)[,1];rho <- fitted(b)[,2]xi <- fitted(b)[,3]## Get the predicted expected response... fv <- mu + exp(rho)*(gamma(1-xi)-1)/xirsd <- residuals(b)plot(fv,rsd);qqnorm(rsd)plot(fv,residuals(b,"pearson"))plot(fv,residuals(b,"response"))

Grouped families

Description

Family for use withgam orbam allowing a univariate response vector to be made up of variables from several different distributions. The response variable is supplied as a 2 column matrix, where the first column contains the response observations and the second column indexes the distribution (family) from which it comes.gfam takes a list of families as its single argument.

Useful for modelling data from different sources that are linked by a model sharing some components. Smooth model components that are not shared are usually handled withby variables (seegam.models).

Usage

gfam(fl)

Arguments

Details

Each component function ofgfam uses the families supplied in the listfl to obtain the required quantities for that family's subset of data, and combines the results appropriately. For example it provides the total deviance (twice negative log-likelihood) of the model, along with its derivatives, by computing the family specific deviance and derivatives from each family applied to its subset of data, and summing them. Other quantities are computed in the same way.

Regular exponential families do not compute the same quantities as extended families, sogfam converts what these families produce toextended.family form internally.

Scale parameters obviously have to be handled separately for each family, and treated as parameters to be estimated, just like otherextended.family non-location distribution parameters. Again this is handled internally. This requirement is part of the reason that anextended.family is always produced, even if all elements offl are standard exponential families. In consequence smoothing parameter estimation is always by REML or NCV.

Note that the null deviance is currently computed by assuming a single parameter model for each family, rather than just one parameter, which may slightly lower explained deviances. Note also that residual checking should probably be done by disaggregating the residuals by family. For this reason functions are not provided to facilitate residual checking withqq.gam.

Prediction on the response scale requires that a family index vector is supplied, with the name of the response, as part of the new prediction data. However, families such asocat which usually produce matrix predictions for prediction type"response", will not be able to do so when part ofgfam.

gfam relies on the methods in Wood, Pya and Saefken (2016).

Value

An object of classextended.family.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Examples

library(mgcv)## a mixed family simulator function to play with...sim.gfam <- function(dist,n=400) {## dist can be norm, pois, gamma, binom, nbinom, tw, ocat (R assumed 4)## links used are identitiy, log or logit.  dat <- gamSim(1,n=n,verbose=FALSE)  nf <- length(dist) ## how many families  fin <- c(1:nf,sample(1:nf,n-nf,replace=TRUE)) ## family index  dat[,6:10] <- dat[,6:10]/5 ## a scale that works for all links used  y <- dat$y;  for (i in 1:nf) {    ii <- which(fin==i) ## index of current family    ni <- length(ii);fi <- dat$f[ii]    if (dist[i]=="norm") {      y[ii] <- fi + rnorm(ni)*.5    } else if (dist[i]=="pois") {      y[ii] <- rpois(ni,exp(fi))    } else if (dist[i]=="gamma") {      scale <- .5      y[ii] <- rgamma(ni,shape=1/scale,scale=exp(fi)*scale)    } else if (dist[i]=="binom") {      y[ii] <- rbinom(ni,1,binomial()$linkinv(fi))    } else if (dist[i]=="nbinom") {      y[ii] <- rnbinom(ni,size=3,mu=exp(fi))    } else if (dist[i]=="tw") {      y[ii] <- rTweedie(exp(fi),p=1.5,phi=1.5)    } else if (dist[i]=="ocat") {      alpha <- c(-Inf,1,2,2.5,Inf)      R <- length(alpha)-1      yi <- fi      u <- runif(ni)      u <- yi + log(u/(1-u))       for (j in 1:R) {        yi[u > alpha[j]&u <= alpha[j+1]] <- j      }      y[ii] <- yi    }  }  dat$y <- cbind(y,fin)  dat} ## sim.gfam## some examplesdat <- sim.gfam(c("binom","tw","norm"))b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),         family=gfam(list(binomial,tw,gaussian)),data=dat)predict(b,data.frame(y=1:3,x0=c(.5,.5,.5),x1=c(.3,.2,.3),        x2=c(.2,.5,.8),x3=c(.1,.5,.9)),type="response",se=TRUE)summary(b)plot(b,pages=1)## set up model so that only the binomial observations depend## on x0...dat$id1 <- as.numeric(dat$y[,2]==1)b1 <- gam(y~s(x0,by=id1)+s(x1)+s(x2)+s(x3),         family=gfam(list(binomial,tw,gaussian)),data=dat)plot(b1,pages=1) ## note the CI width increase

GAM Integrated Nested Laplace Approximation Newton Enhanced

Description

Apply Integrated Nested Laplace Approximation (INLA, Rue et al. 2009) to models estimable bygam orbam, using the INLA variant described in Wood (2019). Produces marginal posterior densities for each coefficient, selected coefficients or linear transformations of the coefficient vector.

Usage

ginla(G,A=NULL,nk=16,nb=100,J=1,interactive=FALSE,integ=0,approx=0)

Arguments

Details

Let\beta,\theta andy denote the model coefficients, hyperparameters/smoothing parameters and response data, respectively. In principle, INLA employs Laplace approximations for\pi(\beta_i|\theta,y) and\pi(\theta|y) and then obtains the marginal posterior distribution\pi(\beta_i|y) by intergrating the approximations to\pi(\beta_i|\theta,y)\pi(\theta|y) w.r.t\theta (marginals for the hyperparameters can also be produced). In practice the Laplace approximation for\pi(\beta_i|\theta,y) is too expensive to compute for each\beta_i and must itself be approximated. To this end, there are two quantities that have to be computed: the posterior mode\beta^*|\beta_i and the determinant of the Hessian of the joint log density\log \pi(\beta,\theta,y) w.r.t.\beta at the mode. Rue et al. (2009) originally approximated the posterior conditional mode by the conditional mode implied by a simple Gaussian approximation to the posterior\pi(\beta|y). They then approximated the log determinant of the Hessian as a function of\beta_i using a first order Taylor expansion, which is cheap to compute for the sparse model representaiton that they use, but not when using the dense low rank basis expansions used bygam. They also offer a more expensive alternative approximation based on computing the log determiannt with respect only to those elements of\beta with sufficiently high correlation with\beta_i according to the simple Gaussian posterior approximation: efficiency again seems to rest on sparsity. Wood (2020) suggests computing the required posterior modes exactly, and basing the log determinant approximation on a BFGS update of the Hessian at the unconditional model. The latter is efficient with or without sparsity, whereas the former is a ‘for free’ improvement. Both steps are efficient because it is cheap to obtain the Cholesky factor ofH[-i,-i] from that ofH - seecholdrop. This is the approach taken by this routine.

Theapprox argument allows two further approximations to speed up computations. Forapprox==1 the exact posterior conditional modes are not used, but instead the conditional modes implied by the simple Gaussian posterior approximation. Forapprox==2 the same approximation is used for the modes and the Hessian is assumed constant. The latter is quite fast as no log joint density gradient evaluations are required.

Note that for many models the INLA estimates are very close to the usual Gaussian approximation to the posterior, theinteractive argument is useful for investigating this issue.

bam models are only supported with thedisrete=TRUE option. Thediscrete=FALSE approach would be too inefficient. AR1 models are not supported (related arguments are simply ignored).

Value

A list with elementsbeta anddensity, both of which are matrices. Each row relates to one coefficient (or linear coefficient combination) of interest. Both matrices havenb columns. Ifint!=0 then a further elementreml gives the integration weights used in the CCD integration, with the central point weight given first.

WARNINGS

This routine is still somewhat experimental, so details are liable to change. Also currently not all steps are optimally efficient.

The routine is written for relatively expert users.

ginla is not designed to deal with rank deficient models.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Rue, H, Martino, S. & Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series B. 71: 319-392.

Wood (2020) Simplified Integrated Laplace Approximation. Biometrika 107(1): 223-230. [Note: There is an error in the theorem proof - theoretical properties are weaker than claimed - under investigation]

Examples

  require(mgcv); require(MASS)  ## example using a scale location model for the motorcycle data. A simple  ## plotting routine is produced first...  plot.inla <- function(x,inla,k=1,levels=c(.025,.1,.5,.9,.975),               lcol = c(2,4,4,4,2),lwd = c(1,1,2,1,1),lty=c(1,1,1,1,1),       xlab="x",ylab="y",cex.lab=1.5) {    ## a simple effect plotter, when distributions of function values of    ## 1D smooths have been computed    require(splines)    p <- length(x)     betaq <- matrix(0,length(levels),p) ## storage for beta quantiles     for (i in 1:p) { ## work through x and betas      j <- i + k - 1      p <- cumsum(inla$density[j,])*(inla$beta[j,2]-inla$beta[j,1])      ## getting quantiles of function values...      betaq[,i] <- approx(p,y=inla$beta[j,],levels)$y    }    xg <- seq(min(x),max(x),length=200)    ylim <- range(betaq)    ylim <- 1.1*(ylim-mean(ylim))+mean(ylim)    for (j in 1:length(levels)) { ## plot the quantiles      din <- interpSpline(x,betaq[j,])      if (j==1) {        plot(xg,predict(din,xg)$y,ylim=ylim,type="l",col=lcol[j],             xlab=xlab,ylab=ylab,lwd=lwd[j],cex.lab=1.5,lty=lty[j])      } else lines(xg,predict(din,xg)$y,col=lcol[j],lwd=lwd[j],lty=lty[j])    }  } ## plot.inla  ## set up the model with a `gam' call...  G <- gam(list(accel~s(times,k=20,bs="ad"),~s(times)),            data=mcycle,family=gaulss(),fit=FALSE)  b <- gam(G=G,method="REML") ## regular GAM fit for comparison  ## Now use ginla to get posteriors of estimated effect values  ## at evenly spaced times. Create A matrix for this...    rat <- range(mcycle$times)  pd0 <- data.frame(times=seq(rat[1],rat[2],length=20))  X0 <- predict(b,newdata=pd0,type="lpmatrix")  X0[,21:30] <- 0  pd1 <- data.frame(times=seq(rat[1],rat[2],length=10))  X1 <- predict(b,newdata=pd1,type="lpmatrix")  X1[,1:20] <- 0  A <- rbind(X0,X1) ## A maps coefs to required function values  ## call ginla. Set integ to 1 for integrated version.  ## Set interactive = 1 or 2 to plot marginal posterior distributions  ## (red) and simple Gaussian approximation (black).   inla <- ginla(G,A,integ=0)  par(mfrow=c(1,2),mar=c(5,5,1,1))  fv <- predict(b,se=TRUE) ## usual Gaussian approximation, for comparison  ## plot inla mean smooth effect...  plot.inla(pd0$times,inla,k=1,xlab="time",ylab=expression(f[1](time)))   ## overlay simple Gaussian equivalent (in grey) ...  points(mcycle$times,mcycle$accel,col="grey")  lines(mcycle$times,fv$fit[,1],col="grey",lwd=2)  lines(mcycle$times,fv$fit[,1]+2*fv$se.fit[,1],lty=2,col="grey",lwd=2)  lines(mcycle$times,fv$fit[,1]-2*fv$se.fit[,1],lty=2,col="grey",lwd=2)  ## same for log sd smooth...  plot.inla(pd1$times,inla,k=21,xlab="time",ylab=expression(f[2](time)))  lines(mcycle$times,fv$fit[,2],col="grey",lwd=2)  lines(mcycle$times,fv$fit[,2]+2*fv$se.fit[,2],col="grey",lty=2,lwd=2)  lines(mcycle$times,fv$fit[,2]-2*fv$se.fit[,2],col="grey",lty=2,lwd=2)  ## ... notice some real differences for the log sd smooth, especially  ## at the lower and upper ends of the time interval.

Gumbel location-scale model family

Description

Thegumbls family implements Gumbel location scale additive models in which the location and scale parameters (see details) can depend on additive smooth predictors. Useable only withgam, the linear predictors are specified via a list of formulae.

Usage

gumbls(link=list("identity","log"),b=-7)

Arguments

Details

Letz = (y-\mu) e^{-\beta}, then the log Gumbel density isl = -\beta - z - e^{-z}. The expected value of a Gumbel r.v. is\mu + \gamma e^{\beta} where\gamma is Eulers constant (about 0.57721566). The corresponding variance is\pi^2 e^{2\beta}/6.

gumbls is used withgam to fit Gumbel location - scale models parameterized in terms of location parameter\mu and the log scale parameter\beta. Note thatidentity link for the scale parameter means that the corresponding linear predictor gives\beta directly. By default thelog link for the scale parameter simply forces the log scale parameter to have a lower limit given by argumentb: if\eta is the linear predictor for the log scale parameter,\beta, then\beta = b + \log(1+e^{\eta-b}).

gam is called with a list containing 2 formulae, the first specifies the response on the left hand side and the structure of the linear predictor for location parameter,\mu, on the right hand side. The second is one sided, specifying the linear predictor for the lg scale,\beta, on the right hand side.

The fitted values for this family will be a two column matrix. The first column is the mean, and the second column is the log scale parameter,\beta. Predictions usingpredict.gam will also produce 2 column matrices fortype"link" and"response". The first column is on the original data scale whentype="response" and on the log mean scale of the linear predictor whentype="link". The second column whentype="response" is again the log scale parameter, but is on the linear predictor whentype="link".

Value

An object inheriting from classgeneral.family.

References

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models.Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Examples

library(mgcv)## simulate some dataf0 <- 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)^10n <- 400;set.seed(9)x0 <- runif(n);x1 <- runif(n);x2 <- runif(n);x3 <- runif(n);mu <- f0(x0)+f1(x1)beta <- exp(f2(x2)/5)y <- mu - beta*log(-log(runif(n))) ## Gumbel quantile functionb <- gam(list(y~s(x0)+s(x1),~s(x2)+s(x3)),family=gumbls)plot(b,pages=1,scale=0)summary(b)gam.check(b)

Identifiability constraints

Description

Smooth terms are generally only identifiable up to an additive constant. In consequence sum-to-zero identifiability constraints are imposed on most smooth terms. The exceptions are terms withby variables which cause the smooth to be identifiable without constraint (that doesn't include factorby variables), and random effect terms. Alternatively smooths can be set up to pass through zero at a user specified point.

Details

By default each smooth term is subject to the sum-to-zero constraint

\sum_i f(x_i) = 0.

The constraint is imposed by reparameterization. The sum-to-zero constraint causes the term to be orthogonal to the intercept: alternative constraints lead to wider confidence bands for the constrained smooth terms.

No constraint is used for random effect terms, since the penalty (random effect covariance matrix) anyway ensures identifiability in this case. Also if aby variable means that the smooth is anyway identifiable, then no extra constraint is imposed. Constraints are imposed for factorby variables, so that the main effect of the factor must usually be explicitly added to the model (the example below is an exception).

Occasionally it is desirable to substitute the constraint that a particular smooth curve should pass through zero at a particular point: thepc argument tos,te,ti andt2 allows this: if specified then such constraints are always applied.

Author(s)

Simon N. Wood (s.wood@r-project.org)

Examples

## Example of three groups, each with a different smooth dependence on x## but each starting at the same value...require(mgcv)set.seed(53)n <- 100;x <- runif(3*n);z <- runif(3*n)fac <- factor(rep(c("a","b","c"),each=100))y <- c(sin(x[1:100]*4),exp(3*x[101:200])/10-.1,exp(-10*(x[201:300]-.5))/       (1+exp(-10*(x[201:300]-.5)))-0.9933071) + z*(1-z)*5 + rnorm(100)*.4## 'pc' used to constrain smooths to 0 at x=0...b <- gam(y~s(x,by=fac,pc=0)+s(z)) plot(b,pages=1)

Which of a set of points lie within a polygon defined region

Description

Tests whether each of a set of points lie within a region defined by one or more (possibly nested) polygons. Points count as ‘inside’ if they are interior to an odd number of polygons.

Usage

in.out(bnd,x)

Arguments

Details

The algorithm works by counting boundary crossings (using compiled C code).

Value

A logical vector of lengthnrow(x).TRUE if the corresponding row ofx is inside the boundary andFALSE otherwise.

Author(s)

Simon N. Woodsimon.wood@r-project.org

Examples

library(mgcv)data(columb.polys)bnd <- columb.polys[[2]]plot(bnd,type="n")polygon(bnd)x <- seq(7.9,8.7,length=20)y <- seq(13.7,14.3,length=20)gr <- as.matrix(expand.grid(x,y))inside <- in.out(bnd,gr)points(gr,col=as.numeric(inside)+1)

Are points inside boundary?

Description

Assesses whether points are inside a boundary. The boundary must enclose thedomain, but may include islands.

Usage

inSide(bnd,x,y,xname=NULL,yname=NULL)

Arguments

Details

Segments of boundary are separated byNAs, or are in separate list elements.The boundary co-ordinates are taken to define nodes which are joined by straight line segments inorder to create the boundary. Each segment is assumed todefine a closed loop, and the last point in a segment will be assumed to bejoined to the first. Loops must not intersect (no test is made forthis).

The method used is to count how many times a line, in the y-direction from apoint, crosses a boundary segment. An odd number of crossings defines aninterior point. Hence in geographic applications it would be usual to havean outer boundary loop, possibly with some inner ‘islands’ completelyenclosed in the outer loop.

The routine calls compiled C code and operates by an exhaustive search foreach point inx, y.

Value

The function returns a logical array of the same dimension asx andy.TRUE indicates that the correspondingx, y point liesinside the boundary.

Author(s)

Simon N. Woodsimon.wood@r-project.org

Examples

require(mgcv)m <- 300;n <- 150xm <- seq(-1,4,length=m);yn <- seq(-1,1,length=n)x <- rep(xm,n);y <- rep(yn,rep(m,n))er <- matrix(fs.test(x,y),m,n)bnd <- fs.boundary()in.bnd <- inSide(bnd,x,y)plot(x,y,col=as.numeric(in.bnd)+1,pch=".")lines(bnd$x,bnd$y,col=3)points(x,y,col=as.numeric(in.bnd)+1,pch=".")## check boundary details ...plot(x,y,col=as.numeric(in.bnd)+1,pch=".",ylim=c(-1,0),xlim=c(3,3.5))lines(bnd$x,bnd$y,col=3)points(x,y,col=as.numeric(in.bnd)+1,pch=".")## alternatively, give the names of x and yd <- data.frame(x, y); rm(x, y)in.bnd2 <- inSide(bnd, d$x, d$y, xname="x", yname="y")all.equal(in.bnd, in.bnd2)

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

Description

Extracts the leading diagonal of the influence matrix (hatmatrix) of a fittedgam object.

Usage

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

Arguments

Details

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

Value

An array (see above).

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

gam

Starting values for multiple smoothing parameter estimation

Description

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

Usage

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

Arguments

Details

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

Usually only called bymagic andgam.

Value

An array of initial smoothing parameter estimates.

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

magic,gam.outer,gam,

Interpret a GAM formula

Description

This is an internal function of packagemgcv. It is a service routine forgam 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, extra.special = NULL)

Arguments

Value

An object of classsplit.gam.formula with the following items:

Author(s)

Simon N. Woodsimon.wood@r-project.org

See Also

gamgamm

Just Another Gibbs Additive Modeller: JAGS support for mgcv.

Description

Facilities to auto-generate model specification code and associated data to simulate with GAMs in JAGS (or BUGS). This is useful for inference about models with complex random effects structure best coded in JAGS. It is a very innefficient approach to making inferences about standard GAMs. The idea is thatjagam generates template JAGS code, and associated data, for the smooth part of the model. This template is then directly edited to include other stochastic components. After simulation with the resulting model, facilities are provided for plotting and prediction with the model smooth components.

Usage

jagam(formula,family=gaussian,data=list(),file,weights=NULL,na.action,offset=NULL,knots=NULL,sp=NULL,drop.unused.levels=TRUE,control=gam.control(),centred=TRUE,sp.prior = "gamma",diagonalize=FALSE)sim2jam(sam,pregam,edf.type=2,burnin=0)

Arguments

Details

Smooths are easily incorportated into JAGS models using multivariate normal priors on the smooth coefficients. The smoothing parameters and smoothing penalty matrices directly specifiy the prior multivariate normal precision matrix. Normally a smoothing penalty does not correspond to a full rank precision matrix, implying an improper prior inappropriate for Gibbs sampling. To rectify this problem the null space penalties suggested in Marra and Wood (2011) are added to the usual penalties.

In an additive modelling context it is usual to centre the smooths, to avoid the identifiability issues associated with having an intercept for each smooth term (in addition to a global intercept). Under Gibbs sampling with JAGS it is technically possible to omit this centring, since we anyway force propriety on the priors, and this propiety implies formal model identifiability. However, in most situations this formal identifiability is rather artificial and does not imply statistically meaningfull identifiability. Rather it serves only to massively inflate confidence intervals, since the multiple intercept terms are not identifiable from the data, but only from the prior. By default then,jagam imposes standard GAM identifiability constraints on all smooths. Thecentred argument does allow you to turn this off, but it is not recommended. If you do setcentred=FALSE then chain convergence and mixing checks should be particularly stringent.

The final technical issue for model setup is the setting of initial conditions for the coefficients and smoothing parameters. The approach taken is to take the default initial smoothing parameter values used elsewhere bymgcv, and to take a single PIRLS fitting step with these smoothing parameters in order to obtain starting values for the smooth coefficients. In the setting of fully conjugate updating the initial values of the coefficients are not critical, and good results are obtained without supplying them. But in the usual setting in which slice sampling is required for at least some of the updates then very poor results can sometimes be obtained without initial values, as the sampler simply fails to find the region of the posterior mode.

Thesim2jam function takes the partialgam object (pregam) fromjagam along with simulation output in standardrjags form and creates a reduced version of agam object, suitable for plotting and prediction of the model's smooth components.sim2gam computes effective degrees of freedom for each smooth, but it should be noted that there are several possibilites for doing this in the context of a model with a complex random effects structure. The simplest approach (edf.type=0) is to compute the degrees of freedom that the smooth would have had if it had been part of an unweighted Gaussian additive model. One might choose to use this option if the model has been modified so that the response distribution and/or link are not those that were specified tojagam. The second option is (edf.type=1) uses the edf that would have been computed bygam had it produced these estimates - in the context in which the JAGS model modifications have all been about modifying the random effects structure, this is equivalent to simply setting all the random effects to zero for the effective degrees of freedom calculation. The default option (edf.type=2) is to base the EDF on the sample covariance matrix,Vp, of the model coefficients. If the simulation output (sim) includes amu field, then this will be used to form the weight matrixW inXWX = t(X)%*%W%*%X, where the EDF is computed fromrowSums(Vp*XWX)*scale. Ifmu is not supplied then it is estimated from the the model matrixX and the mean of the simulated coefficients, but the resultingW may not be strictly comaptible with theVp matrix in this case. In the situation in which the fitted model is very different in structure from the regression model of the template produced byjagam then the default option may make no sense, and indeed it may be best to use option 0.

Value

Forjagam a three item list containing

Forsim2jam an object of class"jam": a partial version of anmgcvgamObject, suitable for plotting and predicting.

WARNINGS

Gibb's sampling is a very slow inferential method for standard GAMs. It is only likely to be worthwhile when complex random effects structures are required above what is possible with direct GAMM methods.

Check that the parameters of the priors on the parameters are fit for your purpose.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood, S.N. (2016) Just Another Gibbs Additive Modeller: Interfacing JAGS and mgcv. Journal of Statistical Software 75(7):1-15 doi:10.18637/jss.v075.i07)

Marra, G. and S.N. Wood (2011) Practical variable selection for generalized additive models.Computational Statistics & Data Analysis 55(7): 2372-2387

Here is a key early reference to smoothing using BUGS (although the approach and smooths used are a bit different to jagam)

Crainiceanu, C. M. D Ruppert, & M.P. Wand (2005) Bayesian Analysis for Penalized Spline Regression Using WinBUGS Journal of Statistical Software 14.

See Also

gam,gamm,bam

Examples

## the following illustrates a typical workflow. To run the ## 'Not run' code you need rjags (and JAGS) to be installed.require(mgcv)  set.seed(2) ## simulate some data... n <- 400dat <- gamSim(1,n=n,dist="normal",scale=2)## regular gam fit for comparison...b0 <- gam(y~s(x0)+s(x1) + s(x2)+s(x3),data=dat,method="REML")## Set directory and file name for file containing jags code.## In real use you would *never* use tempdir() for this. It is## only done here to keep CRAN happy, and avoid any chance of## an accidental overwrite. Instead you would use## setwd() to set an appropriate working directory in which## to write the file, and just set the file name to what you## want to call it (e.g. "test.jags" here). jags.file <- paste(tempdir(),"/test.jags",sep="") ## Set up JAGS code and data. In this one might want to diagonalize## to use conjugate samplers. Usually call 'setwd' first, to set## directory in which model file ("test.jags") will be written.jd <- jagam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,file=jags.file,            sp.prior="gamma",diagonalize=TRUE)## In normal use the model in "test.jags" would now be edited to add ## the non-standard stochastic elements that require use of JAGS....## Not run: require(rjags)load.module("glm") ## improved samplers for GLMs often worth loadingjm <-jags.model(jags.file,data=jd$jags.data,inits=jd$jags.ini,n.chains=1)list.samplers(jm)sam <- jags.samples(jm,c("b","rho","scale"),n.iter=10000,thin=10)jam <- sim2jam(sam,jd$pregam)plot(jam,pages=1)jampd <- data.frame(x0=c(.5,.6),x1=c(.4,.2),x2=c(.8,.4),x3=c(.1,.1))fv <- predict(jam,newdata=pd)## and some minimal checking...require(coda)effectiveSize(as.mcmc.list(sam$b))## End(Not run)## a gamma example...set.seed(1); n <- 400dat <- gamSim(1,n=n,dist="normal",scale=2)scale <- .5; Ey <- exp(dat$f/2)dat$y <- rgamma(n,shape=1/scale,scale=Ey*scale)jd <- jagam(y~s(x0)+te(x1,x2)+s(x3),data=dat,family=Gamma(link=log),            file=jags.file,sp.prior="log.uniform")## In normal use the model in "test.jags" would now be edited to add ## the non-standard stochastic elements that require use of JAGS....## Not run: require(rjags)## following sets random seed, but note that under JAGS 3.4 many## models are still not fully repeatable (JAGS 4 should fix this)jd$jags.ini$.RNG.name <- "base::Mersenne-Twister" ## setting RNGjd$jags.ini$.RNG.seed <- 6 ## how to set RNG seedjm <-jags.model(jags.file,data=jd$jags.data,inits=jd$jags.ini,n.chains=1)list.samplers(jm)sam <- jags.samples(jm,c("b","rho","scale","mu"),n.iter=10000,thin=10)jam <- sim2jam(sam,jd$pregam)plot(jam,pages=1)jampd <- data.frame(x0=c(.5,.6),x1=c(.4,.2),x2=c(.8,.4),x3=c(.1,.1))fv <- predict(jam,newdata=pd)## End(Not run)

Checking smooth basis dimension

Description

Takes a fittedgam object produced bygam() and runs diagnostic tests of whether the basis dimension choises are adequate.

Usage

k.check(b, subsample=5000, n.rep=400)

Arguments

Details

The test of whether the basis dimension for a smooth is adequate (Wood, 2017, section 5.9) is based on computing an estimate of the residual variance based on differencing residuals that are near neighbours according to the (numeric) covariates of the smooth. This estimate divided by the residual variance is thek-index reported. The further below 1 this is, the more likely it is that there is missed pattern left in the residuals. Thep-value is computed by simulation: the residuals are randomly re-shuffledn.rep times to obtain the null distribution of the differencing variance estimator, if there is no pattern in the residuals. For models fitted to more thansubsample data, the tests are based ofsubsample randomly sampled data. Low p-values may indicate that the basis dimension,k, has been set too low, especially if the reportededf is close tok\', the maximum possible EDF for the term. Note the disconcerting fact that if the test statistic itself is based on random resampling and the null is true, then the associated p-values will of course vary widely from one replicate to the next. Currently smooths of factor variables are not supported and will give anNA p-value.

Doubling a suspectk and re-fitting is sensible: if the reportededf increases substantially then you may have been missing something in the first fit. Of course p-values can be low for reasons other than a too lowk. Seechoose.k for fuller discussion.

Value

A matrix contaning the output of the tests described above.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapmanand Hall/CRC Press.

See Also

choose.k,gam,gam.check

Examples

library(mgcv)set.seed(0)dat <- gamSim(1,n=200)b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)plot(b,pages=1)k.check(b)

Log Tweedie density evaluation

Description

A function to evaluate the log of the Tweedie density for variance powers between 1 and 2, inclusive.Also evaluates first and second derivatives of log density w.r.t. its scale parameter,phi, andp, or w.r.t.rho=log(phi) andtheta wherep = (a+b*exp(theta))/(1+exp(theta)).

Usage

ldTweedie(y,mu=y,p=1.5,phi=1,rho=NA,theta=NA,a=1.001,b=1.999,all.derivs=FALSE)

Arguments

Details

A Tweedie random variable with 1<p<2 is a sum ofN gamma random variables whereN has a Poisson distribution. The p=1 case is a generalization of a Poisson distribution and is a discrete distribution supported on integer multiples of the scale parameter. For 1<p<2 the distribution is supported on the positive reals with a point mass at zero. p=2 is a gamma distribution. As p gets very close to 1 the continuous distribution begins to converge on the discretely supported limit at p=1.

ldTweedie is based on the series evaluation method of Dunn and Smyth (2005). Without the restriction onp the calculation of Tweedie densities is less straightforward. If you really need this case then thetweedie package is the place to start.

Therho,theta parameterization is useful for optimization ofp andphi, in order to keeppbounded well away from 1 and 2, andphi positive. The derivatives nearp=1 tend to infinity.

Note that ifp andphi (ortheta andrho) both contain only a single unique value, then the underlyingcode is able to use buffering to avoid repeated calls to expensive log gamma, di-gamma and tri-gamma functions (mu can still be a vector of different values). This is much faster than is possible when these parameters are vectors with different values.

Value

A matrix with 6 columns, or 10 ifall.derivs=TRUE. The first is the log density ofy (log probability ifp=1). The second and third are the first and second derivatives of the log density w.r.t.phi. 4th and 5th columns are first and second derivative w.r.t.p, final column is second derivative w.r.t.phi andp.

Ifrho andtheta were supplied then derivatives are w.r.t. these. In this case, and ifall.derivs=TRUE then the 7th colmn is the derivative w.r.t.mu, the 8th is the 2nd derivative w.r.t.mu, the 9th is the mixed derivative w.r.t.theta andmu and the 10th is the mixed derivative w.r.t.rho andmu.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Dunn, P.K. and G.K. Smith (2005) Series evaluation of Tweedie exponential dispersion model densities. Statistics and Computing 15:267-280

Tweedie, M. C. K. (1984). An index which distinguishes betweensome important exponential families. Statistics: Applications andNew Directions. Proceedings of the Indian Statistical InstituteGolden Jubilee International Conference (Eds. J. K. Ghosh and J.Roy), pp. 579-604. Calcutta: Indian Statistical Institute.

Examples

  library(mgcv)  ## convergence to Poisson illustrated  ## notice how p>1.1 is OK  y <- seq(1e-10,10,length=1000)  p <- c(1.0001,1.001,1.01,1.1,1.2,1.5,1.8,2)  phi <- .5  fy <- exp(ldTweedie(y,mu=2,p=p[1],phi=phi)[,1])  plot(y,fy,type="l",ylim=c(0,3),main="Tweedie density as p changes")  for (i in 2:length(p)) {    fy <- exp(ldTweedie(y,mu=2,p=p[i],phi=phi)[,1])    lines(y,fy,col=i)  }

Getting log generalized determinant of penalty matrices

Description

INTERNAL function calculating the log generalized determinant of penalty matrix S stored blockwise in an Sl list(which is the output ofSl.setup).

Usage

ldetS(Sl, rho, fixed, np, root = FALSE,Stot=FALSE,repara = TRUE,      nt = 1,deriv=2,sparse=FALSE)

Arguments

Value

A list containing:

• ldetS: the log-determinant of S.

• ldetS1: the gradient of the log-determinant of S.

• ldetS2: the Hessian of the log-determinant of S.

• Sl: with modified rS terms, if needed and rho added to each block

• rp: a re-parameterization list.

• E: a total penalty square root such thatt(E)%*%E = S_tot (ifroot==TRUE).

• S: the total penalty matrix (ifStot==TRUE).

Author(s)

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

Linear functionals of a smooth in GAMs

Description

gam allows the response variable to depend on linear functionals of smooth terms. Specifically dependancies of the form

g(\mu_i) = \ldots + \sum_j L_{ij} f(x_{ij}) + \ldots

are allowed, where thex_{ij} are covariate values and theL_{ij} are fixed weights. i.e. the response can depend on the weighted sum of the same smooth evaluated at different covariate values. This allows, for example, for the response to depend on the derivatives or integrals of a smooth (approximated by finite differencing or quadrature, respectively). It also allows dependence on predictor functions (sometimes called ‘signal regression’).

The mechanism by which this is achieved is to supply matrices of covariate values to the model smooth terms specified bys orte terms in the model formula.Each column of the covariate matrix gives rise to a corresponding column of predictions from the smooth. Let the resulting matrix of evaluated smooth values be F (F will have the same dimension as the covariate matrices). In the absense of aby variable then these columns are simply summed and added to the linear predictor. i.e. the contribution of the term to the linear predictor isrowSums(F). If aby variable is present then it must be a matrix, L,say, of the same dimension as F (and the covariate matrices), and it contains the weightsL_{ij} in the summation given above. So in this case the contribution to the linear predictor isrowSums(L*F).

Note that if a{\bf L1} (i.e.rowSums(L)) is a constant vector, or there is noby variable then the smooth will automatically be centred in order to ensure identifiability. Otherwise it will not be. Note also that for centred smooths it can be worth replacing the constant term in the model withrowSums(L) in order to ensure that predictions are automatically on the right scale.

predict.gam can accept matrix predictors for prediction with such terms, in which case itsnewdata argument will need to be a list. However when predicting from the model it is not necessary to provide matrix covariate andby variable values. For example to simply examine the underlying smooth function one would use vectors of covariate values and vectorby variables, with theby variable and equivalent ofL1, above, set to vectors of ones.

The mechanism is usable with random effect smooths which take factor arguments, by using a trick to create a 2D array of factors. Simply create a factor vector containing the columns of the factor matrix stacked end to end (column major order). Then reset the dimensions of this vector to create the appropriate 2D array: the first dimension should be the number of response data and the second the number of columns of the required factor matrix. You can not usematrix ordata.matrix to set up the required matrix of factor levels. See example below.

Author(s)

Simon N. Woodsimon.wood@r-project.org

Examples

### matrix argument `linear operator' smoothinglibrary(mgcv)set.seed(0)################################# simple summation example...################################n<-400sig<-2x <- runif(n, 0, .9)f2 <- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10x1 <- x + .1f <- f2(x) + f2(x1)  ## response is sum of f at two adjacent x values y <- f + rnorm(n)*sigX <- matrix(c(x,x1),n,2) ## matrix covariate contains both x valuesb <- gam(y~s(X))         plot(b)  ## reconstruction of fplot(f,fitted(b))## example of prediction with summation convention...predict(b,list(X=X[1:3,]))## example of prediction that simply evaluates smooth (no summation)...predict(b,data.frame(X=c(.2,.3,.7))) ######################################################################## Simple random effect model example.## model: y[i] = f(x[i]) + b[k[i]] - b[j[i]] + e[i]## k[i] and j[i] index levels of i.i.d. random effects, b.######################################################################set.seed(7)n <- 200x <- runif(n) ## a continuous covariate## set up a `factor matrix'...fac <- factor(sample(letters,n*2,replace=TRUE))dim(fac) <- c(n,2)## simulate data from such a model...nb <- length(levels(fac))b <- rnorm(nb)y <- 20*(x-.3)^4 + b[fac[,1]] - b[fac[,2]] + rnorm(n)*.5L <- matrix(-1,n,2);L[,1] <- 1 ## the differencing 'by' variable mod <- gam(y ~ s(x) + s(fac,by=L,bs="re"),method="REML")gam.vcomp(mod)plot(mod,page=1)## example of prediction using matrices...dat <- list(L=L[1:20,],fac=fac[1:20,],x=x[1:20],y=y[1:20])predict(mod,newdata=dat)######################################################################## multivariate integral example. Function `test1' will be integrated# ## (by midpoint quadrature) over 100 equal area sub-squares covering # ## the unit square. Noise is added to the resulting simulated data.  ### `test1' is estimated from the resulting data using two alternative### smooths.                                                          #######################################################################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))  }## create quadrature (integration) grid, in useful orderig <- 5 ## integration grid within squaremx <- mz <- (1:ig-.5)/igix <- rep(mx,ig);iz <- rep(mz,rep(ig,ig))og <- 10 ## observarion gridmx <- mz <- (1:og-1)/ogox <- rep(mx,og);ox <- rep(ox,rep(ig^2,og^2))oz <- rep(mz,rep(og,og));oz <- rep(oz,rep(ig^2,og^2))x <- ox + ix/og;z <- oz + iz/og ## full grid, subsquare by subsquare## create matrix covariates...X <- matrix(x,og^2,ig^2,byrow=TRUE)Z <- matrix(z,og^2,ig^2,byrow=TRUE)## create simulated test data...dA <- 1/(og*ig)^2  ## quadrature square areaF <- test1(X,Z)    ## evaluate on gridf <- rowSums(F)*dA ## integrate by midpoint quadraturey <- f + rnorm(og^2)*5e-4 ## add noise## ... so each y is a noisy observation of the integral of `test1'## over a 0.1 by 0.1 sub-square from the unit square## Now fit model to simulated data...L <- X*0 + dA## ... let F be the matrix of the smooth evaluated at the x,z values## in matrices X and Z. rowSums(L*F) gives the model predicted## integrals of `test1' corresponding to the observed `y'L1 <- rowSums(L) ## smooths are centred --- need to add in L%*%1## fit models to reconstruct `test1'....b <- gam(y~s(X,Z,by=L)+L1-1)   ## (L1 and const are confounded here)b1 <- gam(y~te(X,Z,by=L)+L1-1) ## tensor product alternative## plot results...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)plot(b)vis.gam(b,view=c("X","Z"),cond=list(L1=1,L=1),plot.type="contour")vis.gam(b1,view=c("X","Z"),cond=list(L1=1,L=1),plot.type="contour")###################################### A "signal" regression example...#####################################rf <- function(x=seq(0,1,length=100)) {## generates random functions...  m <- ceiling(runif(1)*5) ## number of components  f <- x*0;  mu <- runif(m,min(x),max(x));sig <- (runif(m)+.5)*(max(x)-min(x))/10  for (i in 1:m) f <- f+ dnorm(x,mu[i],sig[i])  f}x <- seq(0,1,length=100) ## evaluation points## example functional predictors...par(mfrow=c(3,3));for (i in 1:9) plot(x,rf(x),type="l",xlab="x")## simulate 200 functions and store in rows of L...L <- matrix(NA,200,100) for (i in 1:200) L[i,] <- rf()  ## simulate the functional predictorsf2 <- function(x) { ## the coefficient function  (0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10)/10 }f <- f2(x) ## the true coefficient functiony <- L%*%f + rnorm(200)*20 ## simulated response data## Now fit the model E(y) = L%*%f(x) where f is a smooth function.## The summation convention is used to evaluate smooth at each value## in matrix X to get matrix F, say. Then rowSum(L*F) gives E(y).## create matrix of eval points for each function. Note that## `smoothCon' is smart and will recognize the duplication...X <- matrix(x,200,100,byrow=TRUE) b <- gam(y~s(X,by=L,k=20)) par(mfrow=c(1,1))plot(b,shade=TRUE);lines(x,f,col=2)

AIC and Log likelihood for a fitted GAM

Description

Function to extract the log-likelihood for a fittedgammodel (note that the models are usually fitted by penalized likelihood maximization). Used byAIC. See details for more information on AIC computation.

Usage

## S3 method for class 'gam'logLik(object,...)

Arguments

Details

Modification oflogLik.glm which corrects the degrees offreedom for use withgam objects.

The function is provided so thatAIC functions correctly withgam objects, and uses the appropriate degrees of freedom (accountingfor penalization). See e.g. Wood, Pya and Saefken (2016) for a derivation ofan appropriate AIC.

Forgaussian family models the MLE of the scale parameter is used. For other familieswith a scale parameter the estimated scale parameter is used. This is usually not exactly the MLE, and is not the simple deviance based estimator used withglm models. This is because the simple deviance based estimator can be badly biased in some cases, for example when a Tweedie distribution is employed with low count data.

There are two possibile AIC's that might be considered for use with GAMs. MarginalAIC is based on the marginal likelihood of the GAM, that is the likelihood based ontreating penalized (e.g. spline) coefficients as random and integrating them out. Thedegrees of freedom is then the number of smoothing/variance parameters + the numberof fixed effects. The problem with Marginal AIC is that marginal likelihoodunderestimates variance components/oversmooths, so that the approach favours simpler modelsexcessively (substituting REML does not work, because REML is not comparable between modelswith different unpenalized/fixed components). Conditional AIC uses the likelihood of allthe model coefficients, evaluated at the penalized MLE. The degrees of freedom to use thenis the effective degrees of freedom for the model. However, Greven and Kneib (2010) showthat the neglect of smoothing parameter uncertainty can lead to this conditional AIC beingexcessively likely to select larger models. Wood, Pya and Saefken (2016) propose a simplecorrection to the effective degrees of freedom to fix this problem.mgcv applies thiscorrection whenever possible: that is when usingML orREML smoothing parameterselection withgam orbam. The correctionis not computable when using the Extended Fellner Schall or BFGS optimizer (since the correction requiresan estimate of the covariance matrix of the log smoothing parameters).

Value

StandardlogLik object: seelogLik.

Author(s)

Simon N. Woodsimon.wood@r-project.org based directly onlogLik.glm

References

Greven, S., and Kneib, T. (2010), On the Behaviour of Marginal andConditional AIC in Linear Mixed Models, Biometrika, 97, 773-789.

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter andmodel selection for general smooth models (with discussion).Journal of the American Statistical Association 111, 1548-1575doi:10.1080/01621459.2016.1180986

Wood S.N. (2017) Generalized Additive Models: An Introduction with R(2nd edition). Chapman and Hall/CRC Press.doi:10.1201/9781315370279

See Also

AIC

Basic linear programming

Description

Solves, for{\bf x}, linear programming problems in the standard form

\min {\bf c}^T {\bf x} \text{ s.t. } {\bf Ax}={\bf b},~~{\bf b}\ge {\bf 0}

or in the form

\min {\bf c}^T {\bf x} \text{ s.t. } {\bf Ax}\ge {\bf b},~~{\bf Cx}= {\bf d}

A lightweight implementation of the simplex method, designed for problems of modest size, not large scale problems requiring sparse methods.

feasible finds starting values meeting the constraints, which is also useful for finding feasible initial values for quadratic programming.

Usage

lp(c,A,b,C=NULL,d=NULL,Bi=NULL,maxit=max(1000, nrow(A) * 10), phase1 = FALSE)feasible(A,b,C=NULL,d=NULL,maxit = max(1000, nrow(A) * 10))

Arguments

Details

This code was written to provide a lightweight method for finding feasible initial coefficients for shape constrained splines, but is suitable for solving general linear programming problems of modest size. Given that it uses dense matrix computations, it is not suitable for large scale problems where it is important to exploit sparcity. Neither is it suitable for the sort of linear programming problems arising from integer programs, where degeneracy may be a substantial problem.

The function uses the simplex method (see e.g. Chapter 13 of Nocedal and Wright, 2006). Computational efficiency is ensured by QR decomposing{\bf A}[,Bi] and then efficiently updating the decomposition, every time the active setBi changes, using Givens based updating schemes broadly similar to those given in Golub and van Loan (2013) 5.1.3 p240 and 6.5.3 p337. Given the QR factorization solution of systems involving{\bf A}[,Bi] is efficient.

Value

A vector. Either the solution of the problem (lp), or a feaible initial vector (feasible). Produces an error if there is no solution ormaxit is exceeded.

WARNINGS

Not designed for large scale problems requiring sparse methods, nor for problems where significant degeneracy is expected.

Author(s)

Simon N. Woodsimon.wood@r-project.org

References

Golub G.H. and C.F. van Loan (2013) Matrix Computations (4th edition) Johns Hopkins

Nocedal J. and S. Wright (2006) Numerical Optimization (2nd edition) Springer

See Also

pcls

Examples

library(mgcv)## very simple linear program...c <- c(-4,-2,0,0)A <- matrix(c(1,2,1,.5,1,0,0,1),2,4)b <- c(5,8)x <- lp(c,A,b,c(3,4));sum(c*x);x ## Bi givenx <- lp(c,A,b);sum(c*x);x ## Bi found automatically## equivalent formulation Ax>bA <- -matrix(c(1,2,1,.5),2,2)b <- -c(5,8)C <- matrix(0,0,2); d <- numeric(0)c <- c(-4,-2)x <- lp(c,A,b,C=C,d=d);sum(c*x);x## equivalent formulation Ax>b, Cx=dc <- c(-4,-2,0)A <- matrix(c(0,-2,0,-.5,1,0),2,3)C <- matrix(c(1,1,1),1,3); d <- 5b <- c(0,-8)x <- lp(c,A,b,C=C,d=d);sum(c*x);x

Size of list elements

Description

Produces a named array giving the size, in bytes, of the elements of a list.

Usage

ls.size(x)

Arguments

Value

A numeric vector giving the size in bytes of each element of the listx. The elements of the array have the same names as the elements of the list. Ifx is not a list then its size in bytes is returned, un-named.

Author(s)

Simon N. Woodsimon.wood@r-project.org

Examples

library(mgcv)b <- list(M=matrix(runif(100),10,10),quote="The world is ruled by idiots because only an idiot would want to rule the world.",fam=binomial())ls.size(b)

Stable Multiple Smoothing Parameter Estimation by GCV or UBRE

Description

Function to efficiently estimate smoothing parameters in generalizedridge 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 1).

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,L=NULL,lsp0=NULL,rank=NULL,H=NULL,C=NULL,      w=NULL,gamma=1,scale=1,gcv=TRUE,ridge.parameter=NULL,      control=list(tol=1e-6,step.half=25,rank.tol=      .Machine$double.eps^0.5),extra.rss=0,n.score=length(y),nthreads=1)

Arguments