Commit bd9ae805 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.7-24

parent 59a2c26e
Package: mgcv
Version: 1.7-23
Version: 1.7-24
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
......@@ -14,7 +14,7 @@ Suggests: nlme (>= 3.1-64), splines, Matrix, parallel
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2013-05-20 17:25:14 UTC; sw283
Packaged: 2013-06-05 07:28:12 UTC; sw283
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2013-05-20 22:23:35
Date/Publication: 2013-06-05 10:34:04
3833a46aaf23732d739f65e4143985ba *DESCRIPTION
c5fca2e3b064e4b18ad5c835193c7e19 *DESCRIPTION
5f405da54d89b1e236e272dad685cc2c *NAMESPACE
26423bcf418f00a9a89c06ddaa0646b5 *R/bam.r
e59ee1842b4546cb64ebe39e6a7f00ee *R/fast-REML.r
bed563ba6663fc3d2fa3e2e35c71130e *R/gam.fit3.r
2d679b4646af9013c3e6b11a61e171a6 *R/fast-REML.r
7873e574394156cd1ebb920c693bf3f0 *R/gam.fit3.r
902657a0ee2dedc3fdfa501bf3b37c5b *R/gam.sim.r
fd470608d9222676fc4c5e8c5e6530dd *R/gamm.r
51fed051fcdeb2d4159389f0a8e9eed2 *R/mgcv.r
6045d40194b509250344c55436845f8e *R/gamm.r
0421a1ed94ac15480a0a5a55e3886638 *R/mgcv.r
2370f3d746b680472cd51ccb9dc10284 *R/plots.r
7c562d025d05df4a5002cd0bd2c0030c *R/smooth.r
bc0821492c2db83bc43ecef857b455f9 *R/soap.r
fb66d6c18398411a99ffcb788b854f13 *R/sparse.r
146e60ee2eaed061319f27ef75a559b6 *changeLog
01a1ab28361fc9a574887a74a9309edd *changeLog
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
88d77139cc983317b6acd8c5f1252ab9 *gnugpl2.txt
......@@ -21,7 +21,7 @@ c51c9b8c9c73f81895176ded39b91394 *man/Predict.matrix.cr.smooth.Rd
e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
3e5e30b44947d3ddf00fcd55be219038 *man/Tweedie.Rd
27b7ae82bab920543222c4428d33cee1 *man/anova.gam.Rd
585d6b6a2e7816df83d5f9193195a233 *man/bam.Rd
7c6fcc602c8278f96ef8cacd28a9c8ea *man/bam.Rd
b385d6d5419d0d6aefe03af1a79d5c4e *man/bam.update.Rd
bb5ec26743d46d7ce1dbf128bceab35a *man/cSplineDes.Rd
9b4d616d1b6c4a46ca77d16cded3f806 *man/choose.k.Rd
......@@ -67,7 +67,7 @@ aba56a0341ba9526a302e39d33aa9042 *man/interpret.gam.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
bc9b89db7e7ff246749551c16f5f1f07 *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
f0363e3309cca00db840ad0fb8359c0f *man/negbin.Rd
bedf887cfc36712ff30e3129720dea96 *man/negbin.Rd
8b766a6ad848b0f1ca469e381ded0169 *man/new.name.Rd
7d8f62e182f7c428cc9d46ddd4d97d43 *man/notExp.Rd
445571cf0b115bf19a2d3474e70f59e7 *man/notExp2.Rd
......@@ -90,14 +90,14 @@ f6f1333be2587ffef5970905b13468ea *man/rig.Rd
7258dfc0149fff020054117fd2ee6bd8 *man/s.Rd
d07f7e4c812e8b798a969beb40dca8e2 *man/slanczos.Rd
b5a06918956fd10f23285b453c05bdb4 *man/smooth.construct.Rd
58198db8810ffe0f285d7393de893860 *man/smooth.construct.ad.smooth.spec.Rd
4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.smooth.spec.Rd
76013feaf70d00976bba0154b6f2c946 *man/smooth.construct.cr.smooth.spec.Rd
e571d0064d328b2c1698c4d06c2eba06 *man/smooth.construct.ds.smooth.spec.Rd
db75c958cbfb561914a3291ab58b9786 *man/smooth.construct.fs.smooth.spec.Rd
772e6c18d25cfaf3d9c194031ee042fe *man/smooth.construct.mrf.smooth.spec.Rd
abe15377f471a2d8957a59c19eeef0bb *man/smooth.construct.ps.smooth.spec.Rd
d202c6718fb1138fdd99e6102250aedf *man/smooth.construct.re.smooth.spec.Rd
b67bb662654b968a583301789ecaaf62 *man/smooth.construct.so.smooth.spec.Rd
53fc438bb9c623715ff323e60ca2c367 *man/smooth.construct.so.smooth.spec.Rd
ce311e544603f1089562747652abce20 *man/smooth.construct.sos.smooth.spec.Rd
8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
574641abfbf14339fbb01948235a2a19 *man/smooth.construct.tp.smooth.spec.Rd
......@@ -131,7 +131,7 @@ cd54024d76a9b53dc17ef26323fc053f *src/Makevars
ad5838098ac3fe9d20443fcad3e9a62d *src/mat.c
de0ae24ea5cb533640a3ab57e0383595 *src/matrix.c
0f8448f67d16668f9027084a2d9a1b52 *src/matrix.h
d6407d2677ddf834dbeda77cd0768cee *src/mgcv.c
cb199160696be2d8550b11502542b951 *src/mgcv.c
2b38840c84b6f8d6c2dce7603df8b1e9 *src/mgcv.h
b9328a6d6aa86380fe52a971deccd3d4 *src/misc.c
7e0ba698a21a01150fda519661ef9857 *src/qp.c
......
......@@ -523,7 +523,7 @@ Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NU
#list(reml=ldS$ldetS,reml1=ldS$ldet1,reml2=ldS$ldet2)
#list(reml=dift$rss,reml1=dift$rss1,reml2=dift$rss2)
#list(reml=dift$bSb,reml1=dift$bSb1,reml2=dift$bSb2)
list(reml=reml,reml1=reml1,reml2=reml2,beta=beta[rp],PP=PP,
list(reml=as.numeric(reml),reml1=reml1,reml2=reml2,beta=beta[rp],PP=PP,
rp=ldS$rp,rss=dift$rss+rss.extra,nobs=nobs)
}
......
......@@ -119,6 +119,8 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
nSp <- length(sp)
if (nSp==0) deriv.sp <- 0 else deriv.sp <- deriv
rank.tol <- .Machine$double.eps*100 ## tolerance to use for rank deficiency
xnames <- dimnames(x)[[2]]
ynames <- if (is.matrix(y))
rownames(y)
......@@ -304,7 +306,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
oo <- .C(C_pls_fit1,y=as.double(z),X=as.double(x[good,]),w=as.double(w),
E=as.double(Sr),Es=as.double(Eb),n=as.integer(sum(good)),
q=as.integer(ncol(x)),rE=as.integer(rows.E),eta=as.double(z),
penalty=as.double(1),rank.tol=as.double(.Machine$double.eps*100))
penalty=as.double(1),rank.tol=as.double(rank.tol))
if (!fisher&&oo$n<0) { ## likelihood indefinite - switch to Fisher for this step
z <- (eta - offset)[good] + (yg - mug)/mevg
......@@ -312,7 +314,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
oo <- .C(C_pls_fit1,y=as.double(z),X=as.double(x[good,]),w=as.double(w),
E=as.double(Sr),Es=as.double(Eb),n=as.integer(sum(good)),
q=as.integer(ncol(x)),rE=as.integer(rows.E),eta=as.double(z),
penalty=as.double(1),rank.tol=as.double(.Machine$double.eps*100))
penalty=as.double(1),rank.tol=as.double(rank.tol))
}
start <- oo$y[1:ncol(x)];
......@@ -525,7 +527,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
V2=as.double(V2),V3=as.double(V3),beta=as.double(coef),D1=as.double(D1),
D2=as.double(D2),P=as.double(dum),P1=as.double(P1),P2=as.double(P2),
trA=as.double(dum),trA1=as.double(trA1),trA2=as.double(trA2),
rV=as.double(rV),rank.tol=as.double(.Machine$double.eps*100),
rV=as.double(rV),rank.tol=as.double(rank.tol),
conv.tol=as.double(control$epsilon),rank.est=as.integer(1),n=as.integer(length(z)),
p=as.integer(ncol(x)),M=as.integer(nSp),Mp=as.integer(Mp),Enrow = as.integer(rows.E),
rSncol=as.integer(rSncol),deriv=as.integer(deriv.sp),
......
......@@ -966,21 +966,42 @@ extract.lme.cov2<-function(b,data,start.level=1)
{ if (!inherits(b,"lme")) stop("object does not appear to be of class lme")
grps <- nlme::getGroups(b) # labels of the innermost groupings - in data frame order
n <- length(grps) # number of data
n.levels <- length(b$groups) # number of levels of grouping
if (n.levels<start.level) ## then examine correlation groups
{ if (is.null(b$modelStruct$corStruct)) n.corlevels <- 0 else
n.levels <- length(b$groups) # number of levels of grouping (random effects only)
# if (n.levels<start.level) { ## then examine correlation groups
if (is.null(b$modelStruct$corStruct)) n.corlevels <- 0 else
n.corlevels <-
length(all.vars(nlme::getGroupsFormula(b$modelStruct$corStruct)))
} else n.corlevels <- 0 ## used only to signal irrelevance
# } else n.corlevels <- 0 ## used only to signal irrelevance
## so at this stage n.corlevels > 0 iff it determines the coarsest grouping
## level if > start.level.
if (n.levels<n.corlevels) ## then cor groups are finest
grps <- nlme::getGroups(b$modelStruct$corStruct) # replace grps (but not n.levels)
if (n.levels<n.corlevels) { ## then cor groups are finest
## correlation groups are awkward. The groups returned by
## getGroups(b$modelStruct$corStruct) are not in data frame
## order, but in a sorted order. This is odd because if you
## take the original uninitialized corStruct and Initialize
## using what is stored in b$data, and then apply getGroups,
## the groups are in data frame order (suggesting one route
## for getting the groups into data frame order). The approach
## taken here simply re-constructs the groups in data frame
## order, thereby avoiding the need to store the uninitialized
## corStruct. Note that covariates of corStructs do not
## cause re-ordering: only grouping variables do that.
## The above has been tested at nlme_3.1-109
## grps <- nlme::getGroups(b$modelStruct$corStruct) # replace grps (but not n.levels)
getGroupsFormula(b$modelStruct$corStruct)
vnames <- all.vars(nlme:::getGroupsFormula(b$modelStruct$corStruct))
## attr(b$modelStruct$corStruct,"formula") # would include covariates
lab <- paste(eval(parse(text=vnames[1]),envir=b$data))
if (length(vnames)>1) for (i in 2:length(vnames)) {
lab <- paste(lab,"/",eval(parse(text=vnames[i]),envir=b$data),sep="")
}
grps <- factor(lab)
}
if (n.levels >= start.level||n.corlevels >= start.level)
{ if (n.levels >= start.level)
Cgrps <- nlme::getGroups(b,level=start.level) # outer grouping labels (dforder)
else Cgrps <- nlme::getGroups(b$modelStruct$corStruct) # ditto
else Cgrps <- grps
#Cgrps <- nlme::getGroups(b$modelStruct$corStruct) # ditto
Cind <- sort(as.numeric(Cgrps),index.return=TRUE)$ix
# Cind[i] is where row i of sorted Cgrps is in original data frame order
rCind <- 1:n; rCind[Cind] <- 1:n
......
......@@ -353,7 +353,7 @@ gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5,with.pen=FALSE)
## so now each unique variable name has an associated array of
## the smooths of which it is an argument, arranged in ascending
## order of dimension.
if (maxDim==1) stop("model has repeated 1-d smooths of same variable.")
if (maxDim==1) warning("model has repeated 1-d smooths of same variable.")
## Now set things up to enable term specific model matrices to be
## augmented with square root penalties, on the fly...
......@@ -368,7 +368,8 @@ gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5,with.pen=FALSE)
}
nobs <- nrow(sm[[1]]$X) ## number of observations
for (d in 2:maxDim) { ## work up through dimensions
# for (d in 2:maxDim) { ## work up through dimensions
for (d in 1:maxDim) { ## work up through dimensions
for (i in 1:m) { ## work through smooths
if (sm[[i]]$dim == d&&sm[[i]]$side.constrain) { ## check for nesting
if (with.pen) X1 <- matrix(c(rep(1,nobs),rep(0,np)),nobs+np,as.integer(intercept)) else
......
** denotes quite substantial/important changes
*** denotes really big changes
1.7-24
* Examples pruned in negbin, smooth.construct.ad.smooth.spec and bam help
files to reduce CRAN checking load.
* gam.side now warns if only repeated 1-D smooths of the same variable are
encountered, but does not halt.
* Bug fix in C code for "cr" basis, that could cause a memory violation during
prediction, when an extrapolation was immediately followed by a prediction
that lay exactly on the upper boundary knot. Thanks to Keith Woolner for
reporting this.
* Fix for bug in fast REML code that could cause bam to fail with ti/te only
models. Thanks to Martijn Wieling.
* Fix of bug in extract.lme.cov2, which could cause gamm to fail when
a correlation structure was nested inside a grouping factor finer than
the finest random effect grouping factor.
* Fix for an interesting feature of lme that getGroups applied to the
corStruct that is part of the fitted lme object returns groups in
sorted order, not data frame order, and without an index from one order
to the other. (Oddly, the same corStruct Initialized outside lme has its
groups in data frame order.) This feature could cause gamm to fail,
complaining that the grouping factors for the correlation did not appear
to be nested inside the grouping structure of the random effects. A
bunch of ordering sensitivity tests have been added to the mgcv test suite.
Thanks to Dave Miller for reporting the bug.
1.7-23
*** Fix of severe bug introduced with R 2.15.2 LAPACK change. The shipped
version of dsyevr can fail to produce orthogonal eigenvectors when
uplo='U' (upper triangule of symmetric matrix used), as opposed to 'L'.
uplo='U' (upper triangle of symmetric matrix used), as opposed to 'L'.
This led to a substantial number of gam smoothing parameter estimation
convergence failures, as the key stabilizing re-parameterization was
substantially degraded. The issue did not affect gaussian additive models
......@@ -29,7 +59,7 @@
essentially unchanged by the change.
* summary.gam now deals gracefully with terms such as "fs" smooths
estimated using gamm, for which p-values can no be computed. (thanks to
estimated using gamm, for which p-values can not be computed. (thanks to
Gavin Simpson).
* gam.check/qq.gam now uses a normal QQ-plot when the model has been fitted
......
......@@ -200,9 +200,11 @@ The negbin family is only supported for the *known theta* case.
\examples{
library(mgcv)
## Some not very large examples...
## Some examples are marked 'Not run' purely to keep
## checking load on CRAN down. Sample sizes are small for
## the same reason.
dat <- gamSim(1,n=40000,dist="normal",scale=20)
dat <- gamSim(1,n=20000,dist="normal",scale=20)
bs <- "cr";k <- 20
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=dat)
......@@ -210,14 +212,15 @@ summary(b)
plot(b,pages=1,rug=FALSE) ## plot smooths, but not rug
plot(b,pages=1,rug=FALSE,seWithMean=TRUE) ## `with intercept' CIs
\dontrun{
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 GCV
summary(ba)
summary(ba)}
## A Poisson example...
k <- 15
dat <- gamSim(1,n=15000,dist="poisson",scale=.1)
dat <- gamSim(1,n=11000,dist="poisson",scale=.1)
system.time(b1 <- 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()))
b1
......@@ -235,8 +238,10 @@ system.time(b2 <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+
## ... first call has startup overheads, repeat shows speed up...
\dontrun{
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(),cluster=cl))
}
fv <- predict(b2,cluster=cl) ## parallel prediction
......@@ -244,11 +249,12 @@ if (!is.null(cl)) stopCluster(cl)
b2
## Sparse smoother example...
\dontrun{
dat <- gamSim(1,n=10000,dist="poisson",scale=.1)
system.time( b3 <- bam(y ~ te(x0,x1,bs="ps",k=10,np=FALSE)+
s(x2,bs="ps",k=30)+s(x3,bs="ps",k=30),data=dat,
method="REML",family=poisson(),sparse=TRUE))
b3
b3}
}
......
......@@ -133,10 +133,11 @@ b2a$family$getTheta()
## unknown theta via outer iteration and AIC search
## over a discrete set of values...
\dontrun{
b3<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=negbin(2:10/2),
data=dat)
plot(b3,pages=1)
print(b3)
print(b3)}
## another example...
set.seed(1)
......
......@@ -91,7 +91,9 @@ b <- gam(y~s(x,bs="ad",k=40,m=5,xt=list(bs="cr")))
plot(b,shade=TRUE,main=round(cor(fitted(b),mu),digits=4))
lines(x,mu-mean(mu),col=2)
## A 2D example....
## A 2D example (marked, 'Not run' purely to reduce
## checking load on CRAN).
\dontrun{
par(mfrow=c(2,2),mar=c(1,1,1,1))
x <- seq(-.5, 1.5, length= 60)
z <- x
......@@ -120,7 +122,7 @@ vis.gam(b1,theta=30,phi=30,ticktype="detailed")
b <- gam(y~s(x,z,bs="ad",k=15,m=3),gamma=1.4)
vis.gam(b,theta=30,phi=30,ticktype="detailed")
cor(fitted(b0),f);cor(fitted(b),f)
}
}
......
......@@ -203,9 +203,10 @@ fsb[[1]]$f <- fs.test(fsb[[1]]$v,fsb[[1]]$w,b=1,exclude=FALSE)
bk <- gam(y~s(v,w,bs="so",xt=list(bnd=fsb,nmax=nmax)),knots=knots)
plot(bk) ## default plot
#############################
## tensor product example....
#############################
##########################################
## tensor product example (marked
## 'Not run' to reduce CRAN checking load)
##########################################
\dontrun{
n <- 10000
v <- runif(n)*5-1;w<-runif(n)*2-1
......
......@@ -303,8 +303,11 @@ void crspl(double *x,int *n,double *xk, int *nk,double *X,double *S, double *F,i
if (xi < kmin||xi>kmax) {
extrapolate=1;
} else if (i>0 && fabs(xlast-xi) < 2*h) { /* use simple direct search */
while (xi<xk[j]&&j>0) j--;
while (xi>xk[j+1] && j+1 < *nk-1) j++;
while (xi <= xk[j] && j > 0) j--;
while (xi > xk[j+1] && j < *nk-2) j++;
/* next line should not be needed, except under dodgy use of
fpu registers during optimization... */
if (j<0) j=0;if (j > *nk-2) j = *nk - 2;
/* now xk[j] <= x[i] <= xk[j+1] */
} else { /* bisection search required */
j=0;jup=*nk-1;
......
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