Commit 72251544 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.3-3

parent 1eaabe88
Package: mgcv
Version: 1.3-2
Version: 1.3-3
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-52), MASS (>= 7.2-2)
LazyLoad: yes
License: GPL version 2 or later
Packaged: Mon Jun 20 10:33:55 2005; simon
Packaged: Thu Jun 30 13:52:32 2005; simon
......@@ -37,7 +37,7 @@ export(anova.gam,coef.pdIdnot,coef.pdTens,corMatrix.pdIdnot,Dim.pdIdnot,
importFrom(graphics, plot)
importFrom(stats, predict, residuals, influence, formula, logLik, anova,
cooks.distance,vcov)
cooks.distance,vcov,coef)
S3method(anova, gam)
S3method(influence, gam)
......@@ -55,6 +55,7 @@ S3method(summary, pdTens)
S3method(summary, pdIdnot)
S3method(vcov,gam)
S3method(smooth.construct, cc.smooth.spec)
S3method(smooth.construct, cr.smooth.spec)
S3method(smooth.construct, tp.smooth.spec)
......
......@@ -4,7 +4,9 @@
gam.fit2 <- function (x, y, sp, S=list(),rS=list(),off, H=NULL, weights =
rep(1, nobs), start = NULL, etastart = NULL,
mustart = NULL, offset = rep(0, nobs), family = gaussian(),
control = gam.control(), intercept = TRUE,deriv=TRUE,gamma=1,scale=1,pearson=FALSE)
control = gam.control(), intercept =
TRUE,deriv=TRUE,gamma=1,scale=1,pearson=FALSE,
printWarn=TRUE)
## deriv, sp, S, H added to arg list.
## need to modify family before call.
{ x <- as.matrix(x)
......@@ -127,8 +129,8 @@ rep(1, nobs), start = NULL, etastart = NULL,
{ d2g <- family$d2link(mug)
dV <- family$dvar(mug)
eta1 <- (x%*%upe$beta1)[good,]
z1 <- ((yg-mug)*d2g*mevg)*eta1
w1 <- (-0.5*w^3/weg*(dV/mevg + 2*var.mug*d2g))*eta1
z1 <- as.vector((yg-mug)*d2g*mevg)*eta1
w1 <- as.vector(-0.5*w^3/weg*(dV/mevg + 2*var.mug*d2g))*eta1
}
ngoodobs <- as.integer(nobs - sum(!good))
## Here a Fortran call has been replaced by update.beta call
......@@ -283,8 +285,8 @@ rep(1, nobs), start = NULL, etastart = NULL,
d2g <- family$d2link(mug)
dV <- family$dvar(mug)
eta1 <- (x%*%upe$beta1)[good,]
temp <- (-dV/V^2*mu.eta.val*(yg-mug)^2-2/V*(yg-mug)*mu.eta.val)*eta1
temp <- as.matrix(temp*weights[good])
temp <- as.vector(-dV/V^2*mu.eta.val*(yg-mug)^2-2/V*(yg-mug)*mu.eta.val)*eta1
temp <- as.matrix(temp*as.vector(weights[good]))
alpha1 <- colSums(temp) # deriv of alpha w.r.t. s.p.s
} else { ## deviance based GCV/UBRE scores
temp <- weights[good]*(yg-mug)*mu.eta.val/V
......@@ -298,16 +300,16 @@ rep(1, nobs), start = NULL, etastart = NULL,
} else UBRE1<-GCV1<-NULL
# end of inserted code
if (!conv)
if (!conv&&printWarn)
warning("Algorithm did not converge")
if (boundary)
if (printWarn&&boundary)
warning("Algorithm stopped at boundary value")
eps <- 10 * .Machine$double.eps
if (family$family == "binomial") {
if (printWarn&&family$family == "binomial") {
if (any(mu > 1 - eps) || any(mu < eps))
warning("fitted probabilities numerically 0 or 1 occurred")
}
if (family$family == "poisson") {
if (printWarn&&family$family == "poisson") {
if (any(mu < eps))
warning("fitted rates numerically 0 occurred")
}
......@@ -353,13 +355,14 @@ gam2derivative <- function(lsp,args)
## For use as optim() objective gradient
{ b<-gam.fit2(x=args$X, y=args$y, sp=lsp, S=args$S,rS=args$rS,off=args$off, H=args$H,
offset = args$offset,family = args$family,weights=args$w,deriv=TRUE,
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson)
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson,
printWarn=FALSE)
if (args$scoreType == "GCV") ret <- b$GCV1 else ret <- b$UBRE1
ret
}
gam2objective <- function(lsp,args)
gam2objective <- function(lsp,args,printWarn=FALSE)
## Performs IRLS GAM fitting for smoothing parameters given in lsp
## and returns the GCV or UBRE score for the model.
## args is a list containing the arguments for gam.fit2
......@@ -367,7 +370,8 @@ gam2objective <- function(lsp,args)
{
b<-gam.fit2(x=args$X, y=args$y, sp=lsp, S=args$S,rS=args$rS,off=args$off, H=args$H,
offset = args$offset,family = args$family,weights=args$w,deriv=FALSE,
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson)
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson,
printWarn=printWarn)
if (args$scoreType == "GCV") ret <- b$GCV else ret <- b$UBRE
attr(ret,"full.fit") <- b
ret
......@@ -381,7 +385,8 @@ gam3objective <- function(lsp,args)
{
b<-gam.fit2(x=args$X, y=args$y, sp=lsp, S=args$S,rS=args$rS,off=args$off, H=args$H,
offset = args$offset,family = args$family,weights=args$w,deriv=TRUE,
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson)
control=args$control,gamma=args$gamma,scale=args$scale,pearson=args$pearson,
printWarn=FALSE)
if (args$scoreType == "GCV") ret <- b$GCV else ret <- b$UBRE
attr(ret,"full.fit") <- b
if (args$scoreType == "GCV") at <- b$GCV1 else at <- b$UBRE1
......
This diff is collapsed.
This diff is collapsed.
30/6/2005 1.3-3
* te() smooths were not always estimated correctly by gamm(): invariance
lost and different results to equivalent s() smooths. The problem seems
to lie in a sensitivity of lme() estimation to the absolute size of the
`S' attribute matrices of a pdTens class pdMat object: the problem did
not occur at the last revision of the pdTens class, and there are no
changes logged for nlme that could have caused it, so I guess it's down
to a change in something that lme calls in the base distribution.
To avoid the problem, smooth.construct.tensor.smooth.spec has been
modified to scale all marginal penalty matrices so that they have
largest singular value 1.
* Changes to GLMs in R 2.1.1 mean that if the response is an array, gam
could fail, due to failure of terms like w * X when w is and array
rather than a vector. Code modified accordingly.
* Outer iteration now suppresses some warnings, until the final fitted
model is obtained, in order to avoid printing warnings that actually
don't apply to the final fit.
* Version number reporting made (hopefully) more robust.
* pdconstruct.pdTens removed absolute lower limit on coef - replaced with
relative lower limit.
* moved tensor product constraint construction to BEFORE by variable
stuff in smooth.construct.tensor.smooth.spec.
......@@ -20,7 +20,7 @@ gam.fit2(x, y, sp, S=list(),rS=list(),off, H=NULL,
weights = rep(1, nobs), start = NULL, etastart = NULL,
mustart = NULL, offset = rep(0, nobs), family = gaussian(),
control = gam.control(), intercept = TRUE,deriv=TRUE,
gamma=1,scale=1,pearson=FALSE)
gamma=1,scale=1,pearson=FALSE,printWarn=TRUE)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......@@ -70,6 +70,10 @@ scores can be varied (usually increased) using this parameter.}
\item{pearson}{The GCV/UBRE score can be based either on the Pearson statistic
or the deviance. The latter is generally to be preferred, as it is less prone
to severe undersmoothing.}
\item{printWarn}{Set to \code{FALSE} to suppress some warnings. Useful in
order to ensure that some warnings are only printed if they apply to the final
fitted model, rather than an intermediate used in optimization.}
}
\details{ This routine is basically \code{\link{glm.fit}} with some
modifications to allow (i) for quadratic penalties on the log likelihood;
......
......@@ -12,7 +12,7 @@ smoothing parameters, in a manner sutable for use by \code{\link{optim}} or \cod
Not normally called directly, but rather service routines for \code{\link{gam.outer}}.
}
\usage{
gam2objective(lsp,args)
gam2objective(lsp,args,printWarn=FALSE)
gam2derivative(lsp,args)
gam3objective(lsp,args)
}
......@@ -21,6 +21,10 @@ gam3objective(lsp,args)
\item{lsp}{The log smoothing parameters.}
\item{args}{List of arguments required to call \code{\link{gam.fit2}}.}
\item{printWarn}{Should \code{\link{gam.fit2}} print some warnings? Used to
suppress warnings that are only of interest for the final fitted model, until
that model is reached.}
}
\details{ \code{gam2objective} and \code{gam2derivative} are functions suitable
for calling by \code{\link{optim}}, to evaluate the GCV/UBRE score and it's
......
......@@ -38,7 +38,12 @@ covariance matrix for these random effects.}
\item{data}{data frame in which to evaluate formula.}
}
\details{ This appears to be the minimum set of functions required to implement a new \code{pdMat} class.
\details{ If using this class directly note that it is worthwhile scaling the
\code{S} matrices to be of `moderate size', for example by dividing each
matrix by its largest singular value: this avoids problems with \code{lme}
defaults (\code{\ref{smooth.construct.tensor.smooth.spec}} does this automatically).
This appears to be the minimum set of functions required to implement a new \code{pdMat} class.
Note that while the \code{pdFactor} and \code{pdMatrix} functions return the inverse of the scaled random
effect covariance matrix or its factor, the \code{pdConstruct} function is
......
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