Commit 6cfd6069 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.5-1

parent 81e6a96c
Package: mgcv
Version: 1.5-0
Version: 1.5-1
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: GAMs with GCV/AIC/REML smoothness estimation and GAMMs by PQL
Description: Routines for GAMs and other generalized ridge regression
with multiple smoothing parameter selection by GCV, REML
or UBRE/AIC. Also GAMMs by REML or PQL. Includes a gam()
function.
Description: Routines for GAMs and other generalized ridge regression
with multiple smoothing parameter selection by GCV, REML or
UBRE/AIC. Also GAMMs by REML or PQL. Includes a gam() function.
Priority: recommended
Depends: R (>= 2.3.0)
Imports: graphics, stats
Suggests: nlme (>= 3.1-64), splines
LazyLoad: yes
License: GPL (>=2)
Packaged: Tue Mar 3 10:50:59 2009; simon
Packaged: Tue Mar 17 10:28:50 2009; simon
Repository: CRAN
Date/Publication: 2009-03-17 18:17:32
## R routines for gam fitting with calculation of derivatives w.r.t. sp.s
## (c) Simon Wood 2004-2008
## (c) Simon Wood 2004-2009
## This routine is for type 3 gam fitting. The basic idea is that a P-IRLS
## These routine is for type 3 gam fitting. The basic idea is that a P-IRLS
## is run to convergence, and only then is a scheme for evaluating the
## derivatives iterated to convergence. The advantage is that many key
## quantities are fixed at this stage, including the key decompositions
## In addition the R side work is simplified considerably.The routine
## evaluates first and second derivatives of the deviance and tr(A).
## derivatives via the implicit function theorem used.
gam.reparam <- function(rS,lsp,deriv)
......@@ -81,7 +78,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
control = gam.control(), intercept = TRUE,deriv=2,
gamma=1,scale=1,printWarn=TRUE,scoreType="REML",null.coef=rep(0,ncol(x)),...)
## Version with new reparameterization and truncation strategy
## Version with new reparameterization and truncation strategy.
## ISSUES:
## 1. More efficient use of re-parameterization strategy could be employed,
## as repara is O(nq^2). Could pass in a repara object, containing sp's
......@@ -509,7 +506,11 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
if (control$trace) cat("done!\n")
rV <- matrix(oo$rV,ncol(x),ncol(x))
rV <- matrix(oo$rV,ncol(x),ncol(x)) ## rV%*%t(rV)*scale gives covariance matrix
Kmat <- matrix(0,nrow(x),ncol(x))
Kmat[good,] <- oo$X ## rV%*%t(K)%*%(sqrt(wf)*X) = F; diag(F) is edf array
coef <- oo$beta;
trA <- oo$trA;
scale.est <- dev/(nobs-trA)
......@@ -637,7 +638,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
names(mu) <- ynames
names(eta) <- ynames
wt <- rep.int(0, nobs)
wt[good] <- w
if (fisher) wt[good] <- w else wt[good] <- wf ## note that Fisher weights are returned
names(wt) <- ynames
names(weights) <- ynames
names(y) <- ynames
......@@ -659,10 +660,20 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
boundary = boundary,D1=D1,D2=D2,P=P,P1=P1,P2=P2,trA=trA,trA1=trA1,trA2=trA2,
GCV=GCV,GCV1=GCV1,GCV2=GCV2,GACV=GACV,GACV1=GACV1,GACV2=GACV2,UBRE=UBRE,
UBRE1=UBRE1,UBRE2=UBRE2,REML=REML,REML1=REML1,REML2=REML2,rV=rV,
scale.est=scale.est,reml.scale= reml.scale,aic=aic.model,rank=oo$rank.est)
scale.est=scale.est,reml.scale= reml.scale,aic=aic.model,rank=oo$rank.est,K=Kmat)
} ## end gam.fit3
gam.fit3.post.proc <- function(X,object) {
## get edf array and covariance matrices after a gam fit.
Vb <- object$rV%*%t(object$rV)*object$scale.est
PKt <- object$rV%*%t(object$K)
edf <- rowSums(PKt*t(sqrt(object$weights)*X))
Ve <- PKt%*%t(PKt)*object$scale
hat <- rowSums(object$K*object$K)
list(Vb=Vb,Ve=Ve,edf=edf,hat=hat)
}
score.transect <- function(ii, x, y, sp, Eb,UrS=list(),
weights = rep(1, length(y)), start = NULL, etastart = NULL,
......
......@@ -1008,7 +1008,7 @@ gam.outer <- function(lsp,fscale,family,control,method,optimizer,criterion,scale
object$sp <- exp(lsp)
return(object)
}
## some preparations for the other methods
## some preparations for the other methods, which all use gam.fit3...
family <- fix.family.link(family)
family <- fix.family.var(family)
......@@ -1070,7 +1070,8 @@ gam.outer <- function(lsp,fscale,family,control,method,optimizer,criterion,scale
if (scale>0) object$scale <- scale else object$scale <- object$scale.est
mv<-magic.post.proc(G$X,object,w=object$weights)
## mv<-magic.post.proc(G$X,object,w=object$weights)
mv <- gam.fit3.post.proc(G$X,object)
object$Vp <- mv$Vb
object$hat<-mv$hat
object$Ve <- mv$Ve
......
This diff is collapsed.
** denotes quite substantial/important changes
*** denotes really big changes
ISSUES
------
[none known]
1.5-1
* The stability of the fitting methods had become substantially greater
than the stability of the edf calculations after fitting. So it was
possible to fit very poor models, and then have non-sensical effective
degrees of freedom reported. This has been fixed by using much more stable
expressions for edf calculation, when outer iteration smoothness
selection is empolyed. (Changes are to gam.fit3, gam.outer and a new
routine gam.fit3.post.proc).
* edfs in version 1.5-0 were calculated using newton, rather than fisher
weights, in the matrix F=(X'WX+S)^{-1}X'WX, the diagonal of which gives
the edf array. The problem with this is that it is possible for X'WX
not to be +ve definite, and then degrees of freedom can be non-sensical.
Fisher weights are always used now (although the original problem is
exceedingly hard to generate an example of).
* The summation convention code could be *very* memory intensive for cases
in which the matrix arguments of a smooth feature many repeated values.
Code now fixed to make much more efficient use of any repeated rows in
matrix arguments. This enables muc larger signal regression problems to
be tackled.
* Some help file fixes.
1.5-0
*** Efficient and general REML/ML smoothing selection implemented.
......
......@@ -361,7 +361,7 @@ plot(b,pages=1)
## ... and do automatic terms selection as well
b1<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,
method="REML",select=TRUE)
plot(b,pages=1)
plot(b1,pages=1)
## set the smoothing parameter for the first term, estimate rest ...
......
\name{gam.fit3}
\alias{gam.fit3}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{P-IRLS GAM estimation with GCV \& UBRE/AIC derivative calculation}
\title{P-IRLS GAM estimation with GCV \& UBRE/AIC or RE/ML derivative calculation}
\description{Estimation of GAM smoothing parameters is most stable if
optimization of the UBRE/AIC or GCV score is outer to the penalized iteratively
optimization of the UBRE/AIC, GCV, GACV, REML or ML score is outer to the penalized iteratively
re-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 GCV and UBRE/AIC scores
smoothing paramaters, and evaluates derivatives of the smoothness selection scores
of the model with respect to the
log smoothing parameters. Calculation of exact derivatives is generally faster
than approximating them by finite differencing, as well as generally improving
the reliability of GCV/UBRE/AIC score minimization.
the reliability of GCV/UBRE/AIC/REML score minimization.
The approach is to run the P-IRLS
to convergence, and only then to iterate for first and second derivatives.
......
......@@ -20,7 +20,8 @@ difference derivatives.
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,...)
gam.outer(lsp,fscale,family,control,method,optimizer,
criterion,scale,gamma,G,...)
}
%- maybe also `usage' for other objects documented here.
\arguments{
......
......@@ -4,7 +4,7 @@
\title{Generalized Additive Mixed Models}
\description{ Fits the specified generalized additive mixed model (GAMM) to
data, by a call to \code{lme} in the normal errors identity link case, or by
a call to \code{gammPGL} (a modification of \code{glmmPQL} from the \code{MASS} library) otherwise.
a call to \code{gammPQL} (a modification of \code{glmmPQL} from the \code{MASS} library) otherwise.
In the latter case estimates are only approximately MLEs. The routine is typically
slower than \code{gam}, and not quite as numerically robust.
......
......@@ -249,7 +249,7 @@ void pivoter(double *x,int *r,int *c,int *pivot, int *col, int *reverse)
If `col' is non zero then columns are un/pivoted, otherwise rows.
Typical applications are to pivot matrices int the same way that a
Typical applications are to pivot matrices in the same way that a
qr decomposition has been pivoted, or to reverse such a pivoting.
*/
{ double *dum,*px,*pd,*pd1,*p,*p1;
......@@ -1168,6 +1168,9 @@ void get_trA2(double *trA,double *trA1,double *trA2,double *P,double *K,double *
* This version uses only K and P, and is for the case where expressions involve weights which
are reciprocal variances, not the squares of weights which are reciprocal standard deviations.
* Note that tr(A) = tr(KK') and it is tempting to view diag(K'K) as giving the edfs
of the parameters, but this seems to be wrong. It gives the edfs for R \beta, where
R is (pseudo) inverse of P.
*/
......@@ -1591,7 +1594,7 @@ void undrop_cols(double *X,int r,int c, int *drop, int n_drop)
for (i = 0;i<n;i++,X--,Xs--) *X = *Xs;
for (i=0;i<r;i++,X--) *X = 0.0; /* insert 0s at col drop[n_drop-1] */
for (k=n_drop-1;k>0;k++) { /* work through drop */
for (k=n_drop-1;k>0;k--) { /* work through drop */
n = (drop[k] - drop[k-1]-1)*r; /* size of block between cols drop[k-1] and drop[k] */
for (i=0;i<n;i++,X--,Xs--) *X = *Xs;
for (i=0;i<r;i++,X--) *X = 0.0; /* insert 0s at col drop[k-1] */
......@@ -1779,7 +1782,7 @@ void gdi1(double *X,double *E,double *Es,double *rS,double *U1,
indicating the order of differentiation.
The arguments of this function point to the following:
* X is and n by q model matrix.
* X is and n by q model matrix. On output this will contain K.
* E is a q by Enrow square root of the total penalty matrix, so E'E=S
* Es is the square root of a "well scaled" version of the total penalty,
suitable for numerical determination of the theoretical rank of the problem.
......@@ -2287,21 +2290,13 @@ void gdi1(double *X,double *E,double *Es,double *rS,double *U1,
if (deriv2) for (j=0;j<*M;j++,p1++,p2++) *p2 += *p1;
}
}
/* unpivot P into rV and PKtz into beta */
/* PKtz into beta */
for (i=0;i< rank;i++) beta[pivot1[i]] = PKtz[i];
undrop_rows(beta,*q,1,drop,n_drop); /* zero rows inserted */
/* note that PP' and hence rV rV' are propto the cov matrix. */
for (p1=P,i=0;i < rank; i++) for (j=0;j<rank;j++,p1++) rV[pivot1[j] + i * rank] = *p1;
undrop_rows(rV,*q,rank,drop,n_drop); /* zero rows inserted */
p0 = rV + *q * rank;p1 = rV + *q * *q;
for (p2=p0;p2<p1;p2++) *p2 = 0.0; /* padding any trailing columns of rV with zeroes */
/* Now get the REML penalty */
if (*REML>0) { /* It's REML */
......@@ -2335,7 +2330,7 @@ void gdi1(double *X,double *E,double *Es,double *rS,double *U1,
if (neg_w) free(Vt);
free(pivot);free(work);free(PKtz);free(zz);
free(pivot1);
if (*deriv) {
free(b1);free(eta1);
free(eta2);
......@@ -2376,15 +2371,35 @@ void gdi1(double *X,double *E,double *Es,double *rS,double *U1,
for (p1=Q,p0=K,j=0;j<rank;j++,p1 += *Enrow) for (i=0;i<*n;i++,p1++,p0++) *p0 = *p1;
free(Q);free(WX);free(tau);
if (*deriv) pivoter(rS,&rank,&ScS,pivot,&FALSE,&FALSE); /* apply the latest pivoting to rows of rS */
free(pivot);
}
if (*REML) i=0; else i = *deriv;
get_trA2(trA,trA1,trA2,P,K,sp,rS,rSncol,Tfk,Tfkm,wf,n,&rank,&rank,M,&i);
/* unpivot P into rV.... */
/* note that PP' and hence rV rV' are propto the cov matrix. */
if (!*fisher) { /* first unpivot rows of P */
pivoter(P,&rank,&rank,pivot,&FALSE,&TRUE); /* unpivoting the rows of P */
free(pivot);
}
for (p1=P,i=0;i < rank; i++) for (j=0;j<rank;j++,p1++) rV[pivot1[j] + i * rank] = *p1;
undrop_rows(rV,*q,rank,drop,n_drop); /* zero rows inserted */
p0 = rV + *q * rank;p1 = rV + *q * *q;
for (p2=p0;p2<p1;p2++) *p2 = 0.0; /* padding any trailing columns of rV with zeroes */
/* Now unpack K into X -- useful for forming F = PK'W^.5X, diag of which is edf vector... */
for (p0=X,p1=K,p2=K + rank * *n;p1<p2;p0++,p1++) *p0 = *p1;
/* fill trailing columns with zero */
for (p0 = X + rank * *n,p1 = X + *q * *n;p0<p1;p0++) *p0 = 0.0;
if (n_drop) free(drop);
free(nulli);
free(pivot1);
free(P);free(K);
if (*deriv) {
free(Tk);free(Tkm);
......
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