Commit fe6ab539 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-13

parent b47683b8
......@@ -2,10 +2,98 @@
*** denotes really big changes
Currently deprecated and liable to be removed:
- bam(...,sparse=TRUE) [1.8-5]
- single penalty tensor product smooths.
- p.type!=0 in summary.gam.
1.8-13
* Added help file ?one.se.rule on the `one standard error rule' for obtaining
smoother models.
* bam(...,discrete=TRUE) no longer complains about more coefficients than data.
* 's', 'te', 'ti' and 't2' modified to allow user to specify that the smooth
should pass through zero at a specified point. See ?identifiability.
* anova.gam modified to use more appropriate reference degrees of freedom
for multiple model call, where possible. Also fixed to allow multiple
formulae models and to use -2*logLik in place of `deviance' for
general.family models.
* offsets allowed with multinomial, ziplss and gaulss families.
* gevlss family implementing generalized extreme value location, sclae and
shape models.
* Faster code used in 'uniquecombs'. Speeds up discretization step in
'bam(...,discrete=TRUE)'. Could still be improved for multi-column case.
* modification to 'smoothCon' to allow resetting of smooth supplied
constraints - enables fix of bug in bam handling of 't2' terms, where
parameterization of penalty and model matrix did not previously match
properly.
* clarification of `exclude' argument to predict.gam in help files.
* modification to 'plot.gam' etc, so that 'ylim' is no longer shifted by
'shift'.
* ylim and ... handling improved for 'fs' plot method (thanks Dave Miller)
* gam.check now recognises RStudio and plots appropriately.
* bam(...,sparse=TRUE) removed - not efficient, because of unavoidability
of dense off diagonal terms in X'X or equivalent. Deprecated since 1.8-5.
* tweak to initial.sp/g to avoid infinite loop in s.p. initialization, in
rather unusual circumstances. Thanks to Mark Bravington.
* bam and gam have `drop.intercept' argument to force the parametric terms not
to include a constant in their span, even when there are factor variables.
* Fix in Vb.corr (2nd order edf correction) for fixed smoothing parameter case.
* added 'all.vars1' to enable terms like s(x$y) in model formulae.
* modification to gam.fit4 to ignore 'start' if it is immediately worse than
'null.coef'.
* cSplineDes can now accept a 'derivs' argument.
* added drop.intercept handling for multiple formulae (mod by Matteo Fasiolo).
* 'gam.side' fix to avoid duplication of smooth model matrices to be tested
against, when checking for numerical degeneracy. Problem could occasionally
cause a failure (especially with bam), when the total matrix to be tested
against ended upo with more columns than rows.
* 4 occurances of as.name("model.frame") changed to quote(stats::model.frame)
* fix in predict.bamd discrete prediction code to be a bit more relaxed about
use of as.factor, etc in formulae.
* fix in predict.gam handling of 'na.action' to avoid trying to get type of
na.action from name of na.action function - this was fragile to nesting
and could cause predict.bam to fail in the presence of NA's.
* fix of gam.outer so that general families (e.g. cox.ph) can have all
their smoothing parameters supplied without then ignoring the penalties!
* fix in multiple formulae handling of fixed smoothing parameters.
* Fix of bug in zlim handling in vis.gam perspective plot with standard
errors. Thanks Petra Kuhnert.
* probit link added to 'jagam' (suggested by Kenneth Knoblauch).
* 'Sl.' routines revised to allow operation with non-linearly parameterized
smoothers.
* bug fix in Hessian computation in gam.fit5 - leading diagonal of Hessian of
log|Hp| could be wrong where Hp is penalized Hessian.
* better use of crossprod in gamlss.gH
1.8-12
......
Package: mgcv
Version: 1.8-12
Version: 1.8-13
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness
......@@ -16,6 +16,6 @@ LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2016-03-02 18:21:25 UTC; sw283
Packaged: 2016-07-20 08:12:21 UTC; sw283
Repository: CRAN
Date/Publication: 2016-03-03 07:57:37
Date/Publication: 2016-07-21 19:18:06
This diff is collapsed.
......@@ -9,7 +9,7 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
gam2objective,
gamm, gam.check, gam.control,gam.fit3,
gam.fit, gam.outer,gam.vcomp, gamSim ,
gaulss,gam.side,get.var,
gaulss,gam.side,get.var,gevlss,
influence.gam,
in.out,inSide,interpret.gam,initial.sp,
jagam,
......@@ -44,7 +44,7 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
residuals.gam,rig,rTweedie,rmvn,
Rrank,s,scat,sdiag,"sdiag<-",
sim2jam,
slanczos,
slanczos,smooth2random,
smoothCon,smooth.construct,smooth.construct2,
smooth.construct.bs.smooth.spec,
smooth.construct.cc.smooth.spec,
......
This diff is collapsed.
......@@ -33,7 +33,9 @@ cox.ph <- function (link = "identity") {
y.order <- order(G$y,decreasing=TRUE)
G$family$data$y.order <- y.order
G$y <- G$y[y.order]
attrX <- attributes(G$X)
G$X <- G$X[y.order,,drop=FALSE]
attributes(G$X) <- attrX
G$w <- G$w[y.order]
})
......@@ -108,7 +110,7 @@ cox.ph <- function (link = "identity") {
rd <- qf <- NULL ## these functions currently undefined for Cox PH
ll <- function(y,X,coef,wt,family,deriv=0,d1b=0,d2b=0,Hp=NULL,rank=0,fh=NULL,D=NULL) {
ll <- function(y,X,coef,wt,family,offset=NULL,deriv=0,d1b=0,d2b=0,Hp=NULL,rank=0,fh=NULL,D=NULL) {
## function defining the cox model log lik.
## Calls C code "coxlpl"
## deriv codes: 0 - evaluate the log likelihood
......@@ -122,7 +124,7 @@ cox.ph <- function (link = "identity") {
## or its Choleski factor
## D is the diagonal pre-conditioning matrix used to obtain Hp
## if Hr is the raw Hp then Hp = D*t(D*Hr)
if (!is.null(offset)&&sum(offset!=0)) stop("cox.ph does not yet handle offsets")
##tr <- sort(unique(y),decreasing=TRUE)
tr <- unique(y)
r <- match(y,tr)
......
......@@ -795,7 +795,7 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) {
y1 <- y + (y == 0)
theta <- if (p == 1) log(y1/mu) else (y1^(1 - p) - mu^(1 - p))/(1 - p)
kappa <- if (p == 2) log(y1/mu) else (y^(2 - p) - mu^(2 - p))/(2 - p)
2 * (y * theta - kappa) * wt
pmax(2 * (y * theta - kappa) * wt,0)
}
Dd <- function(y, mu, theta, wt, level=0) {
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -427,7 +427,7 @@ smooth2random.t2.smooth <- function(object,vnames,type=1) {
## 2. rind: and index vector such that if br is the vector of
## random coefficients for the term, br[rind] is the coefs in
## order for this term. rinc - dummy here.
## 3. A matrix, U, that transforms coefs, in order [rand1, rand2,... fix]
## 3. A matrix, trans.D, that transforms coefs, in order [rand1, rand2,... fix]
## back to original parameterization. If null, then not needed.
## 4. A matrix Xf for the fixed effects, if any.
## 5. fixed TRUE/FALSE if its fixed or not. If fixed the other stuff is
......@@ -483,7 +483,7 @@ smooth2random.mgcv.smooth <- function(object,vnames,type=1) {
## Returns 1. a list of random effects, including grouping factors, and
## a fixed effects matrix. Grouping factors, model matrix and
## model matrix name attached as attributes, to each element.
## 2. rind: and index vector such that if br is the vector of
## 2. rind: an index vector such that if br is the vector of
## random coefficients for the term, br[rind] is the coefs in
## order for this term. rinc - dummy here.
## 3. A matrix, U, + vec D that transforms coefs, in order [rand1, rand2,... fix]
......@@ -546,7 +546,7 @@ smooth2random.tensor.smooth <- function(object,vnames,type=1) {
## Returns 1. a list of random effects, including grouping factors, and
## a fixed effects matrix. Grouping factors, model matrix and
## model matrix name attached as attributes, to each element.
## 2. rind: and index vector such that if br is the vector of
## 2. rind: an index vector such that if br is the vector of
## random coefficients for the term, br[rind] is the coefs in
## order for this term. rinc - dummy here.
## 3. A matrix, U, that transforms coefs, in order [rand1, rand2,... fix]
......@@ -1201,7 +1201,10 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
# random terms. correlation describes the correlation structure. This routine is basically an interface
# between the basis constructors provided in mgcv and the gammPQL routine used to estimate the model.
{ if (inherits(family,"extended.family")) warning("family are not designed for use with gamm!")
## lmeControl turns sigma=NULL into sigma=0, but if you supply sigma=0 rejects it,
## which will break the standard handling of the control list. Following line fixes.
## but actually Martin M has now fixed lmeControl...
##if (!is.null(control$sigma)&&control$sigma==0) control$sigma <- NULL
control <- do.call("lmeControl",control)
# check that random is a named list
if (!is.null(random))
......@@ -1248,7 +1251,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
mf$min.sp <- mf$H <- mf$gamma <- mf$fit <- mf$niterPQL <- mf$verbosePQL <- mf$G <- mf$method <- mf$... <- NULL
}
mf$drop.unused.levels <- drop.unused.levels
mf[[1]] <- as.name("model.frame")
mf[[1]] <- quote(stats::model.frame) ## as.name("model.frame")
pmf <- mf
gmf <- eval(mf, parent.frame()) # the model frame now contains all the data, for the gam part only
gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly
......@@ -1268,7 +1271,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
## summarize the *raw* input variables
## note can't use get_all_vars here -- buggy with matrices
vars <- all.vars(gp$fake.formula[-2]) ## drop response here
vars <- all.vars1(gp$fake.formula[-2]) ## drop response here
inp <- parse(text = paste("list(", paste(vars, collapse = ","),")"))
dl <- eval(inp, data, parent.frame())
names(dl) <- vars ## list of all variables needed
......
......@@ -12,8 +12,8 @@ write.jagslp <- function(resp,family,file,use.weights,offset=FALSE) {
## write the JAGS code for the linear predictor
## and response distribution.
iltab <- ## table of inverse link functions
c("eta[i]","exp(eta[i])","ilogit(eta[i])","1/eta[i]","eta[i]^2")
names(iltab) <- c("identity","log","logit","inverse","sqrt")
c("eta[i]","exp(eta[i])","ilogit(eta[i])","phi(eta[i])","1/eta[i]","eta[i]^2")
names(iltab) <- c("identity","log","logit","probit","inverse","sqrt")
if (!family$link%in%names(iltab)) stop("sorry link not yet handled")
## code linear predictor and expected response...
......@@ -114,7 +114,7 @@ sp.prior = "gamma",diagonalize=FALSE) {
mf$family <- mf$knots <- mf$sp <- mf$file <- mf$control <-
mf$centred <- mf$sp.prior <- mf$diagonalize <- NULL
mf$drop.unused.levels <- drop.unused.levels
mf[[1]]<-as.name("model.frame")
mf[[1]] <- quote(stats::model.frame) ##as.name("model.frame")
pmf <- mf
pmf$formula <- gp$pf
......@@ -165,7 +165,7 @@ sp.prior = "gamma",diagonalize=FALSE) {
## get initial values, for use by JAGS, and to guess suitable values for
## uninformative priors...
lambda <- initial.spg(G$X,G$y,G$w,family,G$S,G$off,G$L) ## initial sp values
lambda <- initial.spg(G$X,G$y,G$w,family,G$S,G$off,offset=G$offset,L=G$L) ## initial sp values
jags.ini <- list()
lam <- if (is.null(G$L)) lambda else G$L%*%lambda
jin <- jini(G,lam)
......
This diff is collapsed.
......@@ -96,11 +96,13 @@ mvn <- function(d=2) {
nbeta <- ncol(G$X)
ntheta <- ydim*(ydim+1)/2 ## number of cov matrix factor params
lpi <- attr(G$X,"lpi")
#offs <- attr(G$X,"offset")
XX <- crossprod(G$X)
G$X <- cbind(G$X,matrix(0,nrow(G$X),ntheta)) ## add dummy columns to G$X
#G$cmX <- c(G$cmX,rep(0,ntheta)) ## and corresponding column means
G$term.names <- c(G$term.names,paste("R",1:ntheta,sep="."))
attr(G$X,"lpi") <- lpi
#offs -> attr(G$X,"offset")
attr(G$X,"XX") <- XX
## pad out sqrt of balanced penalty matrix to account for extra params
attr(G$Sl,"E") <- cbind(attr(G$Sl,"E"),matrix(0,nbeta,ntheta))
......@@ -157,7 +159,7 @@ mvn <- function(d=2) {
##rd <- qf <- NULL ## these functions currently undefined for
ll <- function(y,X,coef,wt,family,deriv=0,d1b=NULL,d2b=NULL,Hp=NULL,rank=0,fh=NULL,D=NULL) {
ll <- function(y,X,coef,wt,family,offset=NULL,deriv=0,d1b=NULL,d2b=NULL,Hp=NULL,rank=0,fh=NULL,D=NULL) {
## function defining the Multivariate Normal model log lik.
## Calls C code "mvn_ll"
## deriv codes: 0 - evaluate the log likelihood
......@@ -171,6 +173,9 @@ mvn <- function(d=2) {
## or its Choleski factor
## D is the diagonal pre-conditioning matrix used to obtain Hp
## if Hr is the raw Hp then Hp = D*t(D*Hr)
if (!is.null(offset)) {
for (i in 1:length(offset)) if (sum(offset[[i]]!=0)) stop("mvn does not yet handle offsets")
}
lpi <- attr(X,"lpi") ## lpi[[k]] is index of model matrix columns for dim k
overlap <- attr(lpi,"overlap") ## do dimensions share terms?
drop <- attr(X,"drop")
......
......@@ -267,7 +267,7 @@ gam.check <- function(b, old.style=FALSE,
napredict(b$na.action, b$linear.predictors[,1]) else
napredict(b$na.action, b$linear.predictors)
## if (b$method%in%c("GCV","GACV","UBRE","REML","ML","P-ML","P-REML","mle.REML","mle.ML","PQL")) {
old.par<-par(mfrow=c(2,2))
if (is.null(.Platform$GUI) || .Platform$GUI != "RStudio") old.par <- par(mfrow=c(2,2))
if (old.style)
qqnorm(resid,...)
else
......@@ -326,7 +326,7 @@ gam.check <- function(b, old.style=FALSE,
cat("indicate that k is too low, especially if edf is close to k\'.\n\n")
printCoefmat(kchck,digits=3);
}
par(old.par)
if (is.null(.Platform$GUI) ||.Platform$GUI != "RStudio") par(old.par)
## } else plot(linpred,resid,xlab="linear predictor",ylab="residuals",...)
} ## end of gam.check
......@@ -675,7 +675,6 @@ plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=
fac <- rep(x$flev,rep(n,nf))
dat <- data.frame(fac,xx)
names(dat) <- c(x$fterm,x$base$term)
# X <- Predict.matrix.fs.interaction(x,dat)
X <- PredictMat(x,dat)
if (is.null(xlab)) xlabel <- x$base$term else xlabel <- xlab
if (is.null(ylab)) ylabel <- label else ylabel <- ylab
......@@ -683,7 +682,8 @@ plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=
main="",x=xx,n=n,nf=nf))
} else { ## produce the plot
ind <- 1:P$n
plot(P$x[ind],P$fit[ind],ylim=range(P$fit),xlab=P$xlab,ylab=P$ylab,type="l")
if(is.null(ylim)) ylim <- range(P$fit)
plot(P$x[ind],P$fit[ind],ylim=ylim,xlab=P$xlab,ylab=P$ylab,type="l",...)
if (P$nf>1) for (i in 2:P$nf) {
ind <- ind + P$n
if (scheme==0) lines(P$x,P$fit[ind],lty=i,col=i) else
......@@ -867,17 +867,17 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (min.r < ylimit[1]) ylimit[1] <- min.r
}
}
if (!is.null(ylim)) ylimit <- ylim
ylimit <- if (is.null(ylim)) ylimit <- ylimit + shift else ylim
## plot the smooth...
if (shade) {
plot(P$x,trans(P$fit+shift),type="n",xlab=P$xlab,ylim=trans(ylimit+shift),
plot(P$x,trans(P$fit+shift),type="n",xlab=P$xlab,ylim=trans(ylimit),
xlim=P$xlim,ylab=P$ylab,main=P$main,...)
polygon(c(P$x,P$x[n:1],P$x[1]),
trans(c(ul,ll[n:1],ul[1])+shift),col = shade.col,border = NA)
lines(P$x,trans(P$fit+shift),...)
} else { ## ordinary plot
plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,ylim=trans(ylimit+shift),xlim=P$xlim,
plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,ylim=trans(ylimit),xlim=P$xlim,
ylab=P$ylab,main=P$main,...)
if (is.null(list(...)[["lty"]])) {
lines(P$x,trans(ul+shift),lty=2,...)
......@@ -936,9 +936,10 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (scale==0&&is.null(ylim)) {
if (partial.resids) ylimit <- range(P$p.resid,na.rm=TRUE) else ylimit <-range(P$fit)
}
if (!is.null(ylim)) ylimit <- ylim
ylimit <- if (is.null(ylim)) ylimit <- ylimit + shift else ylim
plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,
ylab=P$ylab,ylim=trans(ylimit+shift),xlim=P$xlim,main=P$main,...)
ylab=P$ylab,ylim=trans(ylimit),xlim=P$xlim,main=P$main,...)
if (rug) {
if (jit) rug(jitter(as.numeric(P$raw)),...)
else rug(as.numeric(P$raw),...)
......@@ -1442,15 +1443,15 @@ vis.gam <- function(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,color="he
lo.col <- "green"
hi.col <- "red"
}
if (!is.null(zlim))
{ if (length(zlim)!=2||zlim[1]>=zlim[2]) stop("Something wrong with zlim")
if (!is.null(zlim)) {
if (length(zlim)!=2||zlim[1]>=zlim[2]) stop("Something wrong with zlim")
min.z<-zlim[1]
max.z<-zlim[2]
} else
{ z.max<-max(fv$fit+fv$se.fit*se,na.rm=TRUE)
z.min<-min(fv$fit-fv$se.fit*se,na.rm=TRUE)
} else {
max.z<-max(fv$fit+fv$se.fit*se,na.rm=TRUE)
min.z<-min(fv$fit-fv$se.fit*se,na.rm=TRUE)
zlim<-c(min.z,max.z)
}
zlim<-c(z.min,z.max)
z<-fv$fit-fv$se.fit*se;z<-matrix(z,n.grid,n.grid)
if (plot.type=="contour") warning("sorry no option for contouring with errors: try plot.gam")
......
This diff is collapsed.
......@@ -74,7 +74,10 @@ Tweedie, M. C. K. (1984). An index which distinguishes between
Golden 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 (2015), Smoothing parameter and model selection for general smooth models. \url{http://arxiv.org/abs/1511.03864}
Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and
model selection for general smooth models.
Journal of the American Statistical Association.
\url{http://arxiv.org/abs/1511.03864}
}
\seealso{\code{\link{ldTweedie}}, \code{\link{rTweedie}}}
......
......@@ -7,10 +7,8 @@
\code{gam} objects. For a single fitted \code{gam} object, Wald tests of
the significance of each parametric and smooth term are performed, so interpretation
is analogous to \code{\link{drop1}} rather than \code{anova.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: this latter approach
should not be use to test the significance of terms which can be penalized
to zero. See details.
rather than a sequential type I ANOVA). Otherwise the fitted models are compared using an analysis of deviance table: 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{
\method{anova}{gam}(object, ..., dispersion = NULL, test = NULL,
......@@ -32,7 +30,8 @@ p-values. See \code{\link{summary.gam}} for details.}
}
\details{ If more than one fitted model is provided than \code{anova.glm} is
used, with the difference in model degrees of freedom being taken as the difference
in effective degress of freedom. The p-values resulting from this are only approximate,
in effective degress of freedom (when possible this is a smoothing parameter uncertainty corrected version).
The p-values resulting from this 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
......@@ -82,7 +81,7 @@ Wood, S.N. (2013b) A simple test for random effects in regression models. Biomet
\author{ Simon N. Wood \email{simon.wood@r-project.org} with substantial
improvements by Henric Nilsson.}
\section{WARNING}{ If models 'a' and 'b' differ only in terms with no un-penalized components then
\section{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=TRUE
......
......@@ -21,8 +21,8 @@ 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,
sparse=FALSE,cluster=NULL,nthreads=NA,gc.level=1,use.chol=FALSE,
samfrac=1,drop.unused.levels=TRUE,G=NULL,fit=TRUE,...)
cluster=NULL,nthreads=NA,gc.level=1,use.chol=FALSE,
samfrac=1,drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,...)
}
%- maybe also `usage' for other objects documented here.
......@@ -129,10 +129,6 @@ there are no breaks in AR1 correlaion.}
\item{discrete}{with \code{method="fREML"} it is possible to discretize covariates for storage and efficiency reasons.
If \code{discrete} is \code{TRUE}, a number or a vector of numbers for each smoother term, then discretization happens. If numbers are supplied they give the number of discretization bins. Experimental at present. }
\item{sparse}{Deprecated. If all smooths are P-splines and all tensor products are of the form \code{te(...,bs="ps",np=FALSE)}
then in principle computation could be made faster using sparse matrix methods, and you could set this to \code{TRUE}.
In practice the speed up is disappointing, and the computation is less well conditioned than the default. See details.}
\item{cluster}{\code{bam} can compute the computationally dominant QR decomposition in parallel using \link[parallel]{parLapply}
from the \code{parallel} package, if it is supplied with a cluster on which to do this (a cluster here can be some cores of a
single machine). See details and example code.
......@@ -160,6 +156,9 @@ involving factor variables you might want to turn this off. Only do so if you kn
\item{fit}{if \code{FALSE} then the model is set up for fitting but not estimated, and an object is returned, suitable for passing as the \code{G} argument to \code{bam}.}
\item{drop.intercept}{Set to \code{TRUE} to force the model to really not have the a constant in the parametric model part,
even with factor variables present.}
\item{...}{further arguments for
passing on e.g. to \code{gam.fit} (such as \code{mustart}). }
......@@ -170,8 +169,7 @@ An object of class \code{"gam"} as described in \code{\link{gamObject}}.
}
\details{ When \code{discrete=FALSE}, \code{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 using \code{\link{predict.gam}}. For each block the factor R,
\details{ When \code{discrete=FALSE}, \code{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 using \code{\link{predict.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.
......@@ -192,21 +190,14 @@ If \code{cluster} is provided as a cluster set up using \code{\link[parallel]{ma
When \code{discrete=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. Parallel computation is controlled using the \code{nthreads} argument. For this method no cluster computation is used, and the \code{parallel} package is not required.
If the deprecated argument \code{sparse=TRUE} then QR updating is replaced by an alternative scheme, in which the model matrix is stored whole
as a sparse matrix. This only makes sense if all smooths are P-splines and all tensor products are of the
form \code{te(...,bs="ps",np=FALSE)}, but no check is made. The computations are then based on the Choleski decomposition of
the crossproduct of the sparse model matrix. Although this crossproduct is nearly dense, sparsity should make its
formation efficient, which is useful as it is the leading order term in the operations count. However there is no benefit
in using sparse methods to form the Choleski decomposition, given that the crossproduct is dense.
In practice the sparse matrix handling overheads mean that modest or no speed ups are produced
by this approach, while the computation is less stable than the default, and the memory footprint often higher
(but please let the author know if you find an example where the speedup is really worthwhile).
}
\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.
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. (2016) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association.
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}
......@@ -218,7 +209,7 @@ Wood, S.N., Goude, Y. & Shaw S. (2015) Generalized additive models for large dat
The routine will be slow if the default \code{"tp"} basis is used.
You must have more unique combinations of covariates than the model has total
Unless discrete=TRUE, you must have more unique combinations of covariates than the model has total
parameters. (Total parameters is sum of basis dimensions plus sum of non-spline
terms less the number of spline terms).
......@@ -226,7 +217,6 @@ This routine is less stable than `gam' for the same dataset.
The negbin family is only supported for the *known theta* case.
AIC computation does not currently take account of an AR1 model, if used.
}
......
......@@ -48,7 +48,9 @@ to the data (e.g. default knot locations).
}
\section{WARNINGS }{
AIC computation does not currently take account of AR model, if used.
}
\seealso{\code{\link{mgcv-package}}, \code{\link{bam}}
}
......
......@@ -2,7 +2,7 @@
\alias{bug.reports.mgcv}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Reporting mgcv bugs.}
\description{\code{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. \code{mgcv} does not have a large development budget, so it is a big help if bug reports follow the following guidlines.
\description{\code{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. \code{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 to \email{simon.wood@r-project.org} with a subject line including \code{mgcv} somewhere, containing
......
......@@ -6,7 +6,7 @@
}
\usage{
cSplineDes(x, knots, ord = 4)
cSplineDes(x, knots, ord = 4, derivs=0)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -16,6 +16,7 @@ cSplineDes(x, knots, ord = 4)
\item{ord}{ order of the basis. 4 is a cubic spline basis. Must be >1.}
\item{derivs}{ order of derivative of the spline to evaluate, between 0 and \code{ord}-1. Recycled to length of \code{x}. }
}
\details{ The routine is a wrapper that sets up a B-spline basis, where the basis functions wrap at the first and
......
......@@ -44,7 +44,10 @@ Hastie and Tibshirani (1990) Generalized Additive Models, Chapman and Hall.
Klein, J.P and Moeschberger, M.L. (2003) Survival Analysis: Techniques for
Censored and Truncated Data (2nd ed.) Springer.
Wood, S.N., N. Pya and B. Saefken (2015), Smoothing parameter and model selection for general smooth models. \url{http://arxiv.org/abs/1511.03864}
Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and
model selection for general smooth models.
Journal of the American Statistical Association.
\url{http://arxiv.org/abs/1511.03864}
}
......
......@@ -43,7 +43,10 @@ given potential presence is modelled with a second linear predictor.
\author{ Simon N. Wood (s.wood@r-project.org) & Natalya Pya
}
\references{
Wood, S.N., N. Pya and B. Saefken (2015), Smoothing parameter and model selection for general smooth models. \url{http://arxiv.org/abs/1511.03864}
Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and
model selection for general smooth models.
Journal of the American Statistical Association.
\url{http://arxiv.org/abs/1511.03864}
}
......
......@@ -93,7 +93,9 @@ indices of the linear predictors which will shre the terms specified on the r.h.
}
\section{WARNING}{
A code{gam} formula should not refer to variables using e.g. \code{dat[["x"]]}.
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}}
......
......@@ -17,7 +17,7 @@ with smoothing parameters selected by GCV/UBRE/AIC/REML or by regression splines
fixed degrees of freedom (mixtures of the two are permitted). Multi-dimensional smooths are
available using penalized thin plate regression splines (isotropic) or tensor product splines
(when an isotropic smooth is inappropriate), and users can add smooths.
Linear functionals of smooths can alos be included in models.
Linear functionals of smooths can also be included in models.
For an overview of the smooths available see \code{\link{smooth.terms}}.
For more on specifying models see \code{\link{gam.models}}, \code{\link{random.effects}} and \code{\link{linear.functional.terms}}. For more on model selection see \code{\link{gam.selection}}. Do read \code{\link{gam.check}} and \code{\link{choose.k}}.
......@@ -31,7 +31,8 @@ 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,...)
fit=TRUE,paraPen=NULL,G=NULL,in.out,drop.unused.levels=TRUE,
drop.intercept=NULL,...)
}
%- maybe also `usage' for other objects documented here.
......@@ -85,7 +86,8 @@ exponential family \code{"REML"} is the default, and the only other option is \c
parameter estimation criterion (given by \code{method}). \code{"perf"} for performance iteration. \code{"outer"}
for the more stable direct approach. \code{"outer"} can use several alternative optimizers, specified in the
second element of \code{optimizer}: \code{"newton"} (default), \code{"bfgs"}, \code{"optim"}, \code{"nlm"}
and \code{"nlm.fd"} (the latter is based entirely on finite differenced derivatives and is very slow). }
and \code{"nlm.fd"} (the latter is based entirely on finite differenced derivatives and is very slow). \code{"efs"}
for the extended Fellner Schall method - only currently available for general families.}
\item{scale}{ If this is positive then it is taken as the known scale parameter. Negative signals that the
scale parameter is unknown. 0 signals that the scale parameter is 1 for Poisson and binomial and unknown otherwise.
......@@ -150,6 +152,9 @@ be estimated by RE/ML.}
\item{drop.unused.levels}{by default unused levels are dropped from factors before fitting. For some smooths
involving factor variables you might want to turn this off. Only do so if you know what you are doing.}
\item{drop.intercept}{Set to \code{TRUE} to force the model to really not have the a constant in the parametric model part,
even with factor variables present. Can be vector when \code{formula} is a list.}
\item{...}{further arguments for
passing on e.g. to \code{gam.fit} (such as \code{mustart}). }
......@@ -386,6 +391,11 @@ 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 vec
gam(G=G) ## it's smoother
## change the smoothness selection method to REML
b0 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,method="REML")
## use alternative plotting scheme, and way intervals include
......
......@@ -21,7 +21,7 @@ Not normally called directly, but rather a service routine for \code{\link{gam}}
}
\usage{
gam.outer(lsp,fscale,family,control,method,optimizer,
criterion,scale,gamma,G,...)
criterion,scale,gamma,G,start=NULL,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -47,7 +47,7 @@ finite differencing is being used.}
\item{G}{List produced by \code{mgcv:::gam.setup}, containing most of what's
needed to actually fit a GAM.}
\item{start}{starting parameter values.}
\item{...}{other arguments, typically for passing on to \code{gam.fit3} (ultimately).}
}
\details{
......
......@@ -222,7 +222,7 @@ will occasionally fail. See details section for suggestions, or try the
\code{gamm} is usually much slower than \code{gam}, and on some platforms you may need to
increase the memory available to R in order to use it with large data sets
(see \code{\link{mem.limits}}).
(see \code{\link{memory.limit}}).
Note that the weights returned in the fitted GAM object are dummy, and not
those used by the PQL iteration: this makes partial residual plots look odd.
......
......@@ -32,7 +32,10 @@ The null deviance reported for this family is the sum of squares of the differen
}
\references{
Wood, S.N., N. Pya and B. Saefken (2015), Smoothing parameter and model selection for general smooth models. \url{http://arxiv.org/abs/1511.03864}
Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and
model selection for general smooth models.
Journal of the American Statistical Association.
\url{http://arxiv.org/abs/1511.03864}
}
......
\name{gevlss}
\alias{gevlss}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Generalized Extreme Value location-scale model family}
\description{The \code{gevlss} family implements Generalized Extreme Value location scale additive models in which the location, scale and shape parameters depend on additive smooth predictors. Usable
only with \code{\link{gam}}, the linear predictors are specified via a list of formulae.
}
\usage{
gevlss(link=list("identity","identity","logit"))
}
\arguments{
\item{link}{three item list specifying the link for the location scale and shape parameters. See details.}
}
\value{
An object inheriting from class \code{general.family}.
}
\details{Used with \code{\link{gam}} to fit Generalized Extreme Value location scale and shape models. \code{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 on the right hand side. The second is one sided, specifying the linear predictor for the log scale parameter on the right hand side. The third is one sided specifying the linear predictor for the shape parameter.
Link functions \code{"identity"} and \code{"log"} are available for the location (mu) parameter. There is no choice of link for the log scale parameter (\eqn{\rho = \log \sigma}{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 existence of an MLE.
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.