Commit d3f21f57 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.3-12

parent b3bce4bc
Package: mgcv
Version: 1.3-11
Version: 1.3-12
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: GAMs with GCV smoothness estimation and GAMMs by REML/PQL
......@@ -12,4 +12,4 @@ Imports: graphics, stats
Suggests: nlme (>= 3.1-64), MASS (>= 7.2-2)
LazyLoad: yes
License: GPL version 2 or later
Packaged: Sat Dec 10 14:35:06 2005; simon
Packaged: Tue Dec 13 21:38:39 2005; simon
......@@ -11,8 +11,8 @@ formXtViX Form component of GAMM covariance matrix
full.score GCV/UBRE score for use within nlm
gam Generalized additive models with integrated
smoothness estimation
gam2objective Objective functions for GAM smoothing
parameter estimation
gam2objective Objective functions for GAM smoothing parameter
estimation
gam.check Some diagnostics for a fitted gam model
gam.control Setting GAM fitting defaults
gam.convergence GAM convergence and performance issues
......@@ -47,8 +47,8 @@ mgcv Multiple Smoothing Parameter Estimation by GCV
mgcv.control Setting mgcv defaults
mgcv-package GAMs with GCV smoothness estimation and GAMMs
by REML/PQL
mono.con Monotonicity constraints for a cubic
regression spline
mono.con Monotonicity constraints for a cubic regression
spline
mroot Smallest square root of matrix
new.name Obtain a name for a new variable that is not
already in use
......@@ -59,10 +59,10 @@ notExp2 Alternative to log parameterization for
null.space.dimension The basis of the space of un-penalized
functions for a TPRS
pcls Penalized Constrained Least Squares Fitting
pdIdnot Overflow proof pdMat class for multiples of
the identity matrix
pdTens Functions implementing a pdMat class for
tensor product smooths
pdIdnot Overflow proof pdMat class for multiples of the
identity matrix
pdTens Functions implementing a pdMat class for tensor
product smooths
place.knots Automatically place a set of knots evenly
through covariate values
plot.gam Default GAM plotting
......@@ -74,8 +74,7 @@ residuals.gam Generalized Additive Model residuals
s Defining smooths in GAM formulae
smoothCon Prediction/Construction wrapper functions for
GAM smooth terms
smooth.construct Constructor functions for smooth terms in a
GAM
smooth.construct Constructor functions for smooth terms in a GAM
step.gam Alternatives to step.gam
summary.gam Summary for a GAM fit
te Define tensor product smooths in GAM formulae
......@@ -83,6 +82,6 @@ tensor.prod.model.matrix
Utility functions for constructing tensor
product smooths
uniquecombs find the unique rows in a matrix
vcov.gam Extract parameter (estimator) covariance
matrix from GAM fit
vcov.gam Extract parameter (estimator) covariance matrix
from GAM fit
vis.gam Visualization of GAM objects
......@@ -866,7 +866,7 @@ new.name <- function(proposed,old.names)
gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=list(),weights=NULL,
subset=NULL,na.action,knots=NULL,control=lmeControl(niterEM=0,optimMethod="L-BFGS-B"),
niterPQL=20,verbosePQL=TRUE,...)
niterPQL=20,verbosePQL=TRUE,method="ML",...)
## NOTE: niterEM modified after changed notLog parameterization - old version
## needed niterEM=3. 10/8/05.
# Routine to fit a GAMM to some data. Fixed and smooth terms are defined in the formula, but the wiggly
......@@ -903,7 +903,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
"+",paste(allvars,collapse="+")))
else mf$formula<-gp$fake.formula
mf$correlation<-mf$random<-mf$family<-mf$control<-mf$scale<-mf$knots<-mf$sp<-mf$weights<-
mf$min.sp<-mf$H<-mf$gamma<-mf$fit<-mf$niterPQL<-mf$verbosePQL<-mf$G<-mf$...<-NULL
mf$min.sp<-mf$H<-mf$gamma<-mf$fit<-mf$niterPQL<-mf$verbosePQL<-mf$G<-mf$method<-mf$...<-NULL
mf$drop.unused.levels<-TRUE
mf[[1]]<-as.name("model.frame")
pmf <- mf
......@@ -924,8 +924,10 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
G<-gamm.setup(gp,pterms=pTerms,data=mf,knots=knots,parametric.only=FALSE,absorb.cons=FALSE)
if (is.null(random)&&G$m==0)
stop("gamm models must have at least 1 smooth with unknown smoothing paraemter or at least one other random effect")
n.sr <- length(G$random) # number of random smooths (i.e. s(...,fx=FALSE,...) terms)
if (is.null(random)&&n.sr==0)
stop("gamm models must have at least 1 smooth with unknown smoothing parameter or at least one other random effect")
g<-as.factor(G$y*0+1)
......@@ -942,7 +944,6 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
}
fixed.formula <- as.formula(fixed.formula)
n.sr <- length(G$random) # number of random smooths (i.e. s(...,fx=FALSE,...) terms)
group.name<-rep("",n.sr)
r.name <- names(G$random)
if (n.sr) for (i in 1:n.sr) # adding the constructed variables to the model frame avoiding name duplication
......@@ -955,7 +956,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
if (!is.null(random)) # add explicit random effects
{ r.m<-length(random)
r.names<-c(names(rand),names(random))
for (i in 1:r.m) rand[[G$m+i]]<-random[[i]]
for (i in 1:r.m) rand[[n.sr+i]]<-random[[i]] # this line had G$m not n.sr <1.3-12
names(rand)<-r.names
}
## need to modify the correlation structure formula, in order that any
......@@ -974,12 +975,12 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
### Actually do fitting ....
if (family$family=="gaussian"&&family$link=="identity"&&
length(offset.name)==0) reml.used <- TRUE else reml.used <- FALSE
length(offset.name)==0) lme.used <- TRUE else lme.used <- FALSE
if (reml.used)
if (lme.used)
{ ## following construction is a work-around for problem in nlme 3-1.52
eval(parse(text=paste("ret$lme<-lme(",deparse(fixed.formula),
",random=rand,data=strip.offset(mf),correlation=correlation,control=control,weights=weights,method=\"ML\")"
",random=rand,data=strip.offset(mf),correlation=correlation,control=control,weights=weights,method=method)"
,sep="" )))
##ret$lme<-lme(fixed.formula,random=rand,data=mf,correlation=correlation,control=control)
} else
......@@ -1032,9 +1033,11 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
k <- 1
if (G$m>0) for (i in 1:G$m) # var.param in reverse term order, but forward order within terms!!
{ n.sp <- length(object$smooth[[i]]$S) # number of s.p.s for this term
if (inherits(object$smooth[[i]],"tensor.smooth"))#&&n.sp>1)
object$sp[k:(k+n.sp-1)] <- notExp2(var.param[(n.v-n.sp+1):n.v])
else object$sp[k:(k+n.sp-1)] <- 1/notExp2(var.param[(n.v-n.sp+1):n.v]) ## ^2 <- reinstate to revert
if (n.sp>0) {
if (inherits(object$smooth[[i]],"tensor.smooth"))
object$sp[k:(k+n.sp-1)] <- notExp2(var.param[(n.v-n.sp+1):n.v])
else object$sp[k:(k+n.sp-1)] <- 1/notExp2(var.param[(n.v-n.sp+1):n.v]) ## ^2 <- reinstate to revert
}
k <- k + n.sp
n.v <- n.v - n.sp
}
......@@ -1105,10 +1108,10 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
object$sig2 <- ret$lme$sigma^2
if (reml.used) { object$fit.method <- "lme";object$method <- "REML"}
if (lme.used) { object$fit.method <- "lme";object$method <- method}
else { object$fit.method <- "glmmPQL";object$method <- "PQL"}
if (!reml.used) Vb<-Vb*length(G$y)/(length(G$y)-G$nsdf)
if (!lme.used||method=="ML") Vb<-Vb*length(G$y)/(length(G$y)-G$nsdf)
object$Vp <- Vb
object$Ve <- Vb%*%Z%*%Vb
......
......@@ -752,8 +752,9 @@ Predict.matrix.tensor.smooth<-function(object,data)
{ m<-length(object$margin)
X<-list()
for (i in 1:m) X[[i]]<-Predict.matrix(object$margin[[i]],data)
if (length(object$XP)>0)
for (i in 1:m) if (!is.null(object$XP[[i]])) X[[i]] <- X[[i]]%*%object$XP[[i]]
mxp <- length(object$XP)
if (mxp>0)
for (i in 1:mxp) if (!is.null(object$XP[[i]])) X[[i]] <- X[[i]]%*%object$XP[[i]]
T <- tensor.prod.model.matrix(X)
if (object$by!="NA") # deal with "by" variable
{ by <- get.var(object$by,data)
......@@ -2030,16 +2031,19 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
}
} ## ... done
# get data from which to predict.....
nd.is.mf <- FALSE # need to flag if supplied newdata is already a model frame
if (newdata.guaranteed==FALSE)
{ if (missing(newdata)) # then "fake" an object suitable for prediction
{ newdata<-object$model
new.data.ok <- FALSE
nd.is.mf <- TRUE
}
else # do an R ``standard'' evaluation to pick up data
{ new.data.ok <- TRUE
if (is.data.frame(newdata)&&!is.null(attr(newdata,"terms"))) # it's a model frame
{ if (sum(!(names(object$model)%in%names(newdata)))) stop(
"newdata is a model.frame: it should contain all required variables\n")
nd.is.mf <- TRUE
} else
{ ## Following is non-standard to allow convenient splitting into blocks
## below, and to allow checking that all variables are in newdata ...
......@@ -2057,7 +2061,8 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
}
}
} else {na.act <- NULL}
if (new.data.ok)
{ ## check factor levels are right ...
names(newdata)->nn # new data names
......@@ -2116,8 +2121,10 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
## implements safe prediction for parametric part as described in
## http://developer.r-project.org/model-fitting-functions.txt
if (new.data.ok)
{ mf <- model.frame(Terms,data,xlev=object$xlevels)
if (!is.null(cl <- attr(object$pterms,"dataClasses"))) .checkMFClasses(cl,mf)
{ if (nd.is.mf) mf <- model.frame(data,xlev=object$xlevels) else
{ mf <- model.frame(Terms,data,xlev=object$xlevels)
if (!is.null(cl <- attr(object$pterms,"dataClasses"))) .checkMFClasses(cl,mf)
}
Xp <- model.matrix(Terms,mf,contrasts=object$contrasts)
} else
{ Xp <- model.matrix(Terms,object$model)
......@@ -3092,16 +3099,19 @@ vis.gam <- function(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,color="he
dnm <- names(list(...))
x$model <- strip.offset(x$model)
#x$model <- strip.offset(x$model)
## ... remove "offset(" and ")" from offset column name
v.names <- row.names(attr(delete.response(x$terms),"factors"))
if (is.null(view)) # get default view if none supplied
{ v.names<-attr(attr(x$model,"terms"),"term.labels")
{ # v.names<-attr(attr(x$model,"terms"),"term.labels") # BUG... too many of these!!
if (length(v.names)<2) stop("Model doesn't seem to have enough terms to do anything useful")
view<-v.names[1:2]
}
if (!sum(view%in%names(x$model))) stop(
paste("view variables must be one of",names(x$model)))
paste(c("view variables must be one of",v.names),collapse=", "))
if (length(unique(x$model[,view[1]]))<=1||length(unique(x$model[,view[2]]))<=1)
stop(paste("View variables must contain more than one value. view = c(",view[1],",",view[2],").",sep=""))
......
1.3-12
* vis.gam could fail if the original model formula contained functions of
covariates, since vis.gam calls predict.gam with a newdata argument
based on the *model frame* of the model object. predict.gam now
recognises that this has happened and doesn't fail if newdata is a model
frame which contains, e.g. log(x) rather than x itself. offset handling
simplified as a result.
* prediction from te smooths could fail becuase of a bug in handling the
list of re-parameterization matrices for 1-D terms in
Predict.matrix.tensor.smooth. Fixed. (tensor product docs also updated)
* gamm did not handle s(...,fx=TRUE) terms properly, due to several
failures to count s(...,fx=FALSE) terms properly if there were fixed
terms present. Now fixed.
* In the gaussian additive mixed model case `gamm' now allows "ML" or
"REML" to be selected (and is slightly more self consistent in
handling the results of the two alternatives).
1.3-11
* added package doc file
......
......@@ -29,7 +29,8 @@ To use this function effectively it helps to be quite familiar with the use of
\usage{
gamm(formula,random=NULL,correlation=NULL,family=gaussian(),
data=list(),weights=NULL,subset=NULL,na.action,knots=NULL,
control=lmeControl(niterEM=0,optimMethod="L-BFGS-B"),niterPQL=20,verbosePQL=TRUE,...)
control=lmeControl(niterEM=0,optimMethod="L-BFGS-B"),niterPQL=20,
verbosePQL=TRUE,method="ML",...)
}
\arguments{
......@@ -104,6 +105,11 @@ version of R does not have the \code{nlminb} optimizer function.}
\item{verbosePQL}{Should PQL report its progress as it goes along?}
\item{method}{Which of \code{"ML"} or \code{"REML"} to use in the Gaussian
additive mixed model case when \code{lme} is called directly. Ignored in the
generalized case (or if the model
has an offset), in which case \code{glmmPQL} is used.}
\item{...}{further arguments for passing on e.g. to \code{lme}}
}
%- maybe also `usage' for other objects documented here.
......@@ -139,7 +145,7 @@ since these use an improved optimizer.
\item{gam}{an object of class \code{gam}, less information relating to
GCV/UBRE model selection. At present this contains enough information to use
\code{predict}, \code{summary} and \code{print} methods and \code{vis.gam},
but not to use e.g. the \code{anova} method function to comapre models.}
but not to use e.g. the \code{anova} method function to compare models.}
\item{lme}{the fitted model object returned by \code{lme} or \code{glmmPQL}. Note that the model formulae and grouping
structures may appear to be rather bizarre, because of the manner in which the GAMM is split up and the calls to
\code{lme} and \code{glmmPQL} are constructed.}
......
\name{smooth.construct}
\name{smooth.construct}
\alias{smooth.construct}
\alias{smooth.construct.tp.smooth.spec}
\alias{smooth.construct.cr.smooth.spec}
......@@ -55,7 +55,9 @@ tensor product should be unpenalized.}
\item{full.call}{The full version of the \code{s} term, with all defaults expanded explicitly.}
\item{label}{A suitable label for use with this term.}
\item{dim}{is the number of covariates of which this smooth is a function.}
\item{null.space.dim}{The dimension of the null space of the wiggliness penalty.}
\item{mp}{\code{TRUE} if multiple penalties are to be used.}
\item{np}{\code{TRUE} if 1-D marginal smooths are to be re-parameterized in terms of
function values.}
}}
\item{data}{a data frame in which the covariates and any \code{by} variable can be found.}
\item{knots}{an optional data frame specifying knot locations for each covariate. If it is null then the knot
......@@ -74,9 +76,14 @@ if the term is to be left un-penalized.}
\item{df}{the degrees of freedom associated with this term (at least when unpenalized).}
Usually the returned object will also include extra information required to define the basis, and used by
\code{\link{Predict.matrix}} methods to make predictions using the basis. See the \code{Details} section for the infomation included for the built in smooth classes.
\code{tensor.smooth} returned objects will additionally have each element of the \code{margin} list updated in the same way.
\code{\link{Predict.matrix}} methods to make predictions using the basis. See the \code{Details} section for the information included for the built in smooth classes.
\code{tensor.smooth} returned objects will additionally have each element of
the \code{margin} list updated in the same way. \code{tensor.smooths} also
have a list, \code{XP}, containing re-parameterization matrices for any 1-D marginal terms
re-parameterized in terms of function values. This list will have \code{NULL}
entries for marginal smooths that are not re-parameterized, and is only long
enough to reach the last re-parameterized marginal in the list.
}
......
......@@ -81,7 +81,8 @@ smooths of a tensor product smooth so that the parameters are function values
at a set of points spread evenly through the range of values of the covariate
of the smooth. This means that the penalty of the tensor product associated
with any particular covariate direction can be interpreted as the penalty of
the appropriate marginal smooth applied in that direction and averaged over the smooth.
the appropriate marginal smooth applied in that direction and averaged over
the smooth. Currently this is only done for marginals of a single variable.
The function does not evaluate the variable arguments.
......@@ -93,26 +94,35 @@ The returned object contains the following items:
\item{margin}{A list of \code{smooth.spec} objects of the type returned by \code{\link{s}},
defining the basis from which the tensor product smooth is constructed.}
\item{term}{An array of text strings giving the names of the covariates that
\item{term}{An array of text strings giving the names of the covariates that
the term is a function of.}
\item{by}{is the name of any \code{by} variable as text (\code{"NA"} for none).}
\item{fx}{ logical array with element for each penalty of the term
\item{fx}{ logical array with element for each penalty of the term
(tensor product smooths have multiple penalties). \code{TRUE} if the penalty is to
be ignored, \code{FALSE}, otherwise. }
\item{p.order}{The order of the t.p.r.s. penalty, or 0 for
auto-selection of the penalty order.}
\item{full.call}{Text for pasting into a string to be converted to a
\item{full.call}{Text for pasting into a string to be converted to a
gam formula, which has the values of function options given explicitly -
this is useful for constructing a fully expanded gam formula which can
be used without needing access to any variables that may have been used
to define k, fx, bs or m in the original call. i.e. this is text which
when parsed and evaluated generates a call to \code{s()} with all the
options spelled out explicitly.}
\item{label}{A suitable text label for this smooth term.}
\item{dim}{The dimension of the smoother - i.e. the number of
covariates that it is a function of.}
\item{mp}{\code{TRUE} is multiple penalties are to be used (default).}
\item{np}{\code{TRUE} to re-parameterize 1-D marginal smooths in terms of function
values (defualt).}
}
......
......@@ -10,12 +10,14 @@ vis.gam(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,color="heat",
\arguments{
\item{x}{a \code{gam} object, produced by \code{gam()}}
\item{view}{an array containing the names of the two predictor variables to be displayed on the
x and y dimensions of the plot. If omitted the first two suitable variables
will be used. Names must be names from \code{names(x$model)}.
\item{view}{an array containing the names of the two main effect terms to be displayed on the
x and y dimensions of the plot. If omitted the first two suitable terms
will be used. Names must be names from
\code{names(x$model)}. i.e. \code{terms} are used, rather than \code{variables}.
}
\item{cond}{a named list of the values to use for the other predictor variables (not in \code{view}). Variables omitted from
\item{cond}{a named list of the values to use for the other predictor terms
(not in \code{view}). Terms omitted from
this list will have their values set to their mean for continuous variables,
or first level for factors. Names must correspond to \code{names(x$model)}.
}
......@@ -48,10 +50,10 @@ labelling to the plots. }
}
\value{Simply produces a plot.}
\description{ Produces perspective or contour plot views of \code{gam} model predictions, fixing all but the values in \code{view} to the
values supplied in \code{cond}.
\description{ Produces perspective or contour plot views of \code{gam} model
predictions, fixing all but the values in \code{view} to the values supplied in \code{cond}.
}
\details{ The x and y limits are determined by the ranges of the variables supplied in \code{view}. If \code{se}<=0 then
\details{ The x and y limits are determined by the ranges of the terms named in \code{view}. If \code{se}<=0 then
a single (height colour coded, by default) surface is produced, otherwise three (by default see-through) meshes are produced at
mean and +/- \code{se} standard errors. Parts of the x-y plane too far from
data can be excluded by setting \code{too.far}
......
......@@ -240,10 +240,10 @@ msgid "family not recognized"
msgstr "family not recognized"
msgid ""
"gamm models must have at least 1 smooth with unknown smoothing paraemter or "
"gamm models must have at least 1 smooth with unknown smoothing parameter or "
"at least one other random effect"
msgstr ""
"gamm models must have at least 1 smooth with unknown smoothing paraemter or "
"gamm models must have at least 1 smooth with unknown smoothing parameter or "
"at least one other random effect"
msgid "At least three knots required in call to mono.con."
......
This diff is collapsed.
......@@ -191,7 +191,7 @@ msgstr ""
msgid "family not recognized"
msgstr ""
msgid "gamm models must have at least 1 smooth with unknown smoothing paraemter or at least one other random effect"
msgid "gamm models must have at least 1 smooth with unknown smoothing parameter or at least one other random effect"
msgstr ""
msgid "At least three knots required in call to mono.con."
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment