Commit 3c9df280 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-4

parent ffd38c32
** denotes quite substantial/important changes
*** denotes really big changes
1.8-4
** JAGS/BUGS support added, enabling auto-generation of code and data
required to used mgcv type GAMs with JAGS. Useful for complex random
effects structures, for example.
* smoothCon failed if selection penalties requested, but term was unpenalized.
Now fixed (no selection penalties on unpenalized terms.)
* gam.check would fail for tensor product smooths with by variables - fixed.
* predict.gam would fail when predicting for more data than the blocksize
but selecting only some terms. Fixed thanks to Scott Kostyshak.
* smoothCon now has an argument `diagonal.penalty' allowing single penalty
smooths to be re-parameterized in order to diagonalize the penalty matrix.
PredictMat is modified to apply the same reparameterization, making it
user transparent. Facilitates the setup of smooths for export to other
packages.
* predict.bam now exported in response to a request from another
package maintainer.
* 1.8 allows some prediction tasks for some families (e.g. cox.ph) to
require response variables to be supplied. NAs in these then messed up
prediction when they were not needed (e.g. if response variables with
NAs were provided to predict.gam for a simple exponential family GAM).
Response NAs now passed to the family specific prediction code, restoring
the previous behaviour for most models. Thanks Casper Wilestofte Berg.
* backend parallel QR code used by gam modified to use a pivoted block
algorithm.
* nthreads argument added to 'bam' to allow for parallel computation
for computations in the main process (serial on any cluster nodes).
e.g. QR based combination of results from cluster nodes is now
parallel.
* fREML computation now partly in parallel (controlled by 'nthreads'
argument to 'bam')
* slanczos now accepts an nt argument allowing parallel computation of
main O(n^2) step.
* fix to newton logic problem, which could cause an attempt to use 'score2'
before definition.
* fix to fREML code which could cause matrix square root to lose dimensions
and cause an error.
* initial.sp could perform very poorly for very low basis dimensions - could
set initial sp to effective infinity.
1.8-3
* Fix of two illegal read/write bugs with extended family models with no
......@@ -49,8 +102,9 @@
* bam modified so that choleski based fitting works properly with rank
deficient model matrix (without regularization).
* fix of 1.8-0 bug - gam prior weights mishandled in computation of cov matrix.
Thanks Fabian Scheipl.
* fix of 1.8-0 bug - gam prior weights mishandled in computation of cov matrix,
resulting in incorrect variance estimates (even without prior weights
specified). Thanks Fabian Scheipl.
1.8-0
......@@ -127,6 +181,11 @@
* 'bfgs' now uses a finite difference approximation for the initial inverse
Hessian.
1.7-29
* Single character change to Makevars file so that openMP multi-threading
actually works.
1.7-28
* exclude.too.far updated to use kd-tree instead of inefficient search for
......
Package: mgcv
Version: 1.8-3
Version: 1.8-4
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: splines, parallel, survival, MASS
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2014-08-29 10:10:03 UTC; sw283
Packaged: 2014-11-27 09:49:02 UTC; sw283
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2014-08-29 22:07:35
Date/Publication: 2014-11-27 13:16:04
c0324227202a615e286343615936dd1b *ChangeLog
149d2f3b36649e882a15d0092cf01f55 *DESCRIPTION
0ea61afceef0bf25b8c9d74103884936 *ChangeLog
1a88c7e031e86b1a4d2f0a3733714046 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
870396939733c58c9d523cbc41f9da59 *NAMESPACE
b5f01a8f420b0d296352d500b75e4c35 *R/bam.r
9d5bb48cc0072b9ba638567986928a73 *R/coxph.r
68be7e47206e4b5ee210fb6de27f8169 *NAMESPACE
27f3c91b6ebcefafbf8f370ea38d08c9 *R/bam.r
52ea7081e235df890d3e51198bcc3092 *R/coxph.r
d6f40f3d2acd3e30f3a59cbc85da6bc3 *R/efam.r
2fb8891700c1a813dc387cdce8dd817b *R/fast-REML.r
e72ee5107056471bd2ac183003ab5e25 *R/gam.fit3.r
a824a6e49020786a1eba60b99939eff1 *R/fast-REML.r
065215b64c502afab215f6e97e9ec915 *R/gam.fit3.r
6586de2772276237525c646362926f8b *R/gam.fit4.r
e63276b78a4ff2566af2e651bfc6aa9c *R/gam.sim.r
bd651b8eec80d83d6e0de5f47715a0de *R/gamlss.r
449b21801a62d0d6b753e753d4861570 *R/gamm.r
ee07c79f181f6cc17edeafc8a33ae51b *R/mgcv.r
ed2cdcfa5c0f772bb75eee723b7294cd *R/misc.r
739156695988beabea60d2e259073821 *R/jagam.r
b5887ec9a16b1842918901bc5e4ecc94 *R/mgcv.r
a0d1da72334aefd546cad58d1b2062b1 *R/misc.r
c287514d201a02f060fee98aab1e93c8 *R/mvam.r
4b0df0a2310fba6b38ea8164f76d0ccd *R/plots.r
3f10d9fbae29ae28bf89499b874f9aa4 *R/smooth.r
8e3b23448498fd03e668165ca3b01a60 *R/plots.r
92302835adba3710624e434ad52173bf *R/smooth.r
90f72ef91dfc2a25c262fc002f781bf7 *R/soap.r
0c5d2acdc985d2af9de6bb736d4df2ab *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
......@@ -28,7 +29,7 @@ e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
f12468625253dd3603907de233762fd6 *man/Rrank.Rd
f6cadda5999c35800fd65c08d6812f7b *man/Tweedie.Rd
80f8763baa4987579e2aa56073a9e94e *man/anova.gam.Rd
9c392eb571b80b00fc6fada9d502b69f *man/bam.Rd
ef84108378322022dfe6ca5d5c9c917a *man/bam.Rd
71bde8b8caa24a36529ce7e0ac3165d8 *man/bam.update.Rd
a2beb811b1093c5e82ef32d7de1f7d32 *man/cSplineDes.Rd
9b4d616d1b6c4a46ca77d16cded3f806 *man/choose.k.Rd
......@@ -44,7 +45,7 @@ c7561a0f2e114247955e212aeabc6ae9 *man/formXtViX.Rd
d87ff6287d7e343ea24d2053f4724b35 *man/formula.gam.Rd
4da4d585b329769eb44f0c7a6e7dd554 *man/fs.test.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
b113877121ece5b7e227cd07138b7623 *man/gam.Rd
a1492616b069844ea05762687473c4fd *man/gam.Rd
fe61dd0efab9e920c17335faf3d5764c *man/gam.check.Rd
a65bc22f606e45d185bc375fbf5698a1 *man/gam.control.Rd
44db24b66ce63bc16d2c8bc3f5b42ac5 *man/gam.convergence.Rd
......@@ -59,7 +60,7 @@ b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
a66a814cc4c6f806e824751fda519ae0 *man/gam2objective.Rd
717401fd7efa3b39d90418a5d1d0c216 *man/gamObject.Rd
a2593990d6a87f7b783d0435062dec02 *man/gamSim.Rd
3acf983ba10da78af067d4d3a1794026 *man/gamm.Rd
1a69eb8822eac2d99f92fee8dfe0a818 *man/gamm.Rd
df4a6db696749986fd5e20586fc9b718 *man/gaulss.Rd
39ae109127110152e97dc9f79e08bb14 *man/get.var.Rd
a2ea1233d43fac89e0cacbc09a8d31e2 *man/in.out.Rd
......@@ -67,6 +68,7 @@ a2ea1233d43fac89e0cacbc09a8d31e2 *man/in.out.Rd
2f222eeeb3d7bc42f93869bf8c2af58a *man/influence.gam.Rd
bafea2eef12fdc819f8ac1fb41d8b914 *man/initial.sp.Rd
c00bcfe2d0b44b2ea955f3934421807c *man/interpret.gam.Rd
ae18c335a5bb7192ad608a5a500c2794 *man/jagam.Rd
07d2c259b9edf164f42935170b4fccd0 *man/ldTweedie.Rd
58e73ac26b93dc9d28bb27c8699e12cf *man/linear.functional.terms.Rd
93035193b0faa32700e1421ce8c1e9f6 *man/logLik.gam.Rd
......@@ -74,8 +76,8 @@ c00bcfe2d0b44b2ea955f3934421807c *man/interpret.gam.Rd
19a3f7da03e76d7d2fc0a4e544564c72 *man/magic.Rd
5169af4be5fccf9fa79b6de08e9ea035 *man/magic.post.proc.Rd
59672d8599211ce9cf44252f92110920 *man/mgcv-FAQ.Rd
f4c4902dc931807a14b54142d56cc0ea *man/mgcv-package.Rd
1448336c74a9512d298f44f709a308b4 *man/mgcv-parallel.Rd
01926b8c2bda6fc74b51f6ec7f5cc620 *man/mgcv-package.Rd
934b90e0845a4971b97f1d67d12d61ec *man/mgcv-parallel.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
bc9b89db7e7ff246749551c16f5f1f07 *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
......@@ -103,7 +105,7 @@ c523210ae95cb9aaa0aaa1c37da1a4c5 *man/residuals.gam.Rd
f6f1333be2587ffef5970905b13468ea *man/rig.Rd
7258dfc0149fff020054117fd2ee6bd8 *man/s.Rd
c438eb3e41cb1f51e72283380145d210 *man/scat.Rd
d07f7e4c812e8b798a969beb40dca8e2 *man/slanczos.Rd
6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
b5a06918956fd10f23285b453c05bdb4 *man/smooth.construct.Rd
4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.smooth.spec.Rd
76013feaf70d00976bba0154b6f2c946 *man/smooth.construct.cr.smooth.spec.Rd
......@@ -118,10 +120,10 @@ a8d5eb4772641567fee15714f01a4727 *man/smooth.construct.so.smooth.spec.Rd
8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
3e6d88ef6a8ab21bd6f120120602dcf6 *man/smooth.construct.tp.smooth.spec.Rd
75a042d74460dd0c6e0c0a95c1f5d3f7 *man/smooth.terms.Rd
3064042473b909ba6bbd160c476059d3 *man/smoothCon.Rd
8120512d52c741f8c27c2689f69b63d3 *man/smoothCon.Rd
b55a396da77559dac553613146633f97 *man/sp.vcov.Rd
83bd8e097711bf5bd0fff09822743d43 *man/spasm.construct.Rd
700699103b50f40d17d3824e35522c85 *man/step.gam.Rd
b9394812e5398ec95787c65c1325a027 *man/step.gam.Rd
f0791d830687d6155efb8a73db787401 *man/summary.gam.Rd
47d22ea4fe93bc0fa0f47bd262979941 *man/t2.Rd
63e2050f9a01b12096d640e87916b78a *man/te.Rd
......@@ -147,18 +149,18 @@ bd5ec325242ca6e41056c2f1f7f05bfe *po/pl.po
56d501c682dc6699ff0803f25533e849 *src/coxph.c
114701c0d407f946ac693da55da864da *src/gdi.c
49af97195accb65adc75620183d39a4c *src/general.h
94a7b06d798739744ef9342a6a7051ec *src/init.c
761581e96a73e1bec00a9553356eaaa6 *src/init.c
6dbd109048a3bd18b67db6230c064c21 *src/magic.c
a6527e26dde591b838d4c3c3f5a4e036 *src/mat.c
f45aa90d4ce295f0deebc54798dad330 *src/mat.c
02122df2193f86e2a08f68fe5b8ff972 *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
7a12ef0cbd8f5312cd67669a36d75ebf *src/mgcv.c
0245c7d4b550d27e6c89ddc00e224a8b *src/mgcv.h
f06b1075027c76e9a7ce483c3d70c6ad *src/mgcv.h
00f8a024faef17f90ed04f01e736df71 *src/misc.c
08c1706ffeec4277c484435a0644b0e3 *src/mvn.c
90fe2ad1997db56f2293ad16c7110f88 *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h
15ed450da804a9dbd0833459af1f2474 *src/soap.c
b3ff61207cd26308abf386d6a0666320 *src/sparse-smooth.c
7512172319ca54cbb8384f8ebc8d20bf *src/tprs.c
75d9bf91f5c2dc536918055b6a00d93f *src/tprs.c
5bd85bf0319a7b7c755cf49c91a7cd94 *src/tprs.h
......@@ -8,19 +8,19 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
gam2objective,
gamm, gam.check, gam.control,gam.fit3,
gam.fit, gam.outer,gam.vcomp, gamSim ,
gaulss,
gaulss,gam.side,get.var,
influence.gam,
in.out,inSide,interpret.gam,
gam.side,
get.var,ldTweedie,
initial.sp,logLik.gam,ls.size,
in.out,inSide,interpret.gam,initial.sp,
jagam,
ldTweedie,
logLik.gam,ls.size,
magic, magic.post.proc, model.matrix.gam,
mono.con, mroot, mvn, nb, negbin, new.name,
mono.con, mroot, mvn, nb, negbin, new.name,
notExp,notExp2,notLog,notLog2,pcls,null.space.dimension,
ocat,
pen.edf,pdIdnot,pdTens,
place.knots, plot.gam, polys.plot,print.anova.gam,
print.gam,print.summary.gam,predict.gam,
print.gam,print.summary.gam,predict.gam,predict.bam,
PredictMat,Predict.matrix,Predict.matrix2,
Predict.matrix.cr.smooth,
Predict.matrix.duchon.spline,
......@@ -40,6 +40,7 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
qq.gam,
residuals.gam,rig,rTweedie,
Rrank,s,scat,
sim2jam,
slanczos,
smoothCon,smooth.construct,smooth.construct2,
smooth.construct.cc.smooth.spec,
......@@ -81,14 +82,18 @@ S3method(formula, gam)
S3method(logLik, gam)
S3method(model.matrix,gam)
S3method(plot, gam)
S3method(plot, jam)
S3method(predict, gam)
S3method(predict, bam)
S3method(predict, jam)
S3method(print, anova.gam)
S3method(print, gam)
S3method(print, jam)
S3method(print, summary.gam)
S3method(residuals, gam)
S3method(summary, gam)
S3method(vcov,gam)
S3method(vcov,jam)
S3method(coef,pdTens)
S3method(pdConstruct,pdTens)
......
This diff is collapsed.
......@@ -97,6 +97,7 @@ cox.ph <- function (link = "identity") {
predict <- function(family,se=FALSE,eta=NULL,y=NULL,
X=NULL,beta=NULL,off=NULL,Vb=NULL) {
## prediction function.
if (sum(is.na(y))>0) stop("NA times supplied for cox.ph prediction")
ii <- order(y,decreasing=TRUE) ## C code expects non-increasing
n <- nrow(X)
oo <- .C("coxpred",as.double(X[ii,]),t=as.double(y[ii]),as.double(beta),as.double(Vb),
......
......@@ -304,7 +304,7 @@ ldetS <- function(Sl,rho,fixed,np,root=FALSE) {
}
} ## end of multi-S block branch
} ## end of block loop
if (root) E <- E[rowSums(abs(E))!=0,] ## drop zero rows.
if (root) E <- E[rowSums(abs(E))!=0,,drop=FALSE] ## drop zero rows.
list(ldetS=ldS,ldet1=d1.ldS,ldet2=d2.ldS,Sl=Sl,rp=rp,E=E)
} ## end ldetS
......@@ -396,7 +396,7 @@ Sl.mult <- function(Sl,A,k = 0,full=TRUE) {
A
} ## end Sl.mult
Sl.termMult <- function(Sl,A,full=FALSE) {
Sl.termMult <- function(Sl,A,full=FALSE,nt=1) {
## returns a list containing the product of each element S of Sl
## with A. If full==TRUE then the results include the zero rows
## otherwise these are stripped out, but in that case each element
......@@ -426,12 +426,16 @@ Sl.termMult <- function(Sl,A,full=FALSE) {
k <- k + 1
if (full) { ## return answer with all zeroes in place
B <- A*0
if (Amat) B[ind,] <- Sl[[b]]$Srp[[i]]%*%A[ind,,drop=FALSE] else
B[ind] <- Sl[[b]]$Srp[[i]]%*%A[ind]
if (Amat) {
B[ind,] <- if (nt==1) Sl[[b]]$Srp[[i]]%*%A[ind,,drop=FALSE] else
pmmult(Sl[[b]]$Srp[[i]],A[ind,,drop=FALSE],nt=nt)
} else B[ind] <- Sl[[b]]$Srp[[i]]%*%A[ind]
SA[[k]] <- B
} else { ## strip zero rows from answer
if (Amat) SA[[k]] <- Sl[[b]]$Srp[[i]]%*%A[ind,,drop=FALSE] else
SA[[k]] <- as.numeric(Sl[[b]]$Srp[[i]]%*%A[ind])
if (Amat) {
SA[[k]] <- if (nt==1) Sl[[b]]$Srp[[i]]%*%A[ind,,drop=FALSE] else
pmmult(Sl[[b]]$Srp[[i]],A[ind,,drop=FALSE],nt=nt)
} else SA[[k]] <- as.numeric(Sl[[b]]$Srp[[i]]%*%A[ind])
attr(SA[[k]],"ind") <- ind
}
} ## end of S loop for block b
......@@ -440,10 +444,10 @@ Sl.termMult <- function(Sl,A,full=FALSE) {
SA
} ## end Sl.termMult
d.detXXS <- function(Sl,PP) {
d.detXXS <- function(Sl,PP,nt=1) {
## function to obtain derivatives of log |X'X+S| given unpivoted PP' where
## P is inverse of R from the QR of the augmented model matrix.
SPP <- Sl.termMult(Sl,PP,full=FALSE) ## SPP[[k]] is S_k PP'
SPP <- Sl.termMult(Sl,PP,full=FALSE,nt=nt) ## SPP[[k]] is S_k PP'
nd <- length(SPP)
d1 <- rep(0,nd);d2 <- matrix(0,nd,nd)
for (i in 1:nd) {
......@@ -494,7 +498,7 @@ Sl.ift <- function(Sl,R,X,y,beta,piv,rp) {
list(rss =sum(rsd^2),bSb=sum(beta*Sb),rss1=rss1,bSb1=bSb1,rss2=rss2,bSb2=bSb2,d1b=db)
} ## end Sl.ift
Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NULL,Mp=0) {
Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NULL,Mp=0,nt=1) {
## fits penalized regression model with model matrix X and
## initialised block diagonal penalty Sl to data in y, given
## log smoothing parameters rho.
......@@ -507,8 +511,8 @@ Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NU
ldS <- ldetS(Sl,rho,fixed,np,root=TRUE)
## apply resulting stable re-parameterization to X...
X <- Sl.repara(ldS$rp,X)
## get pivoted QR decomp of augmented model matrix
qrx <- qr(rbind(X,ldS$E),LAPACK=TRUE)
## get pivoted QR decomp of augmented model matrix (in parallel if nt>1)
qrx <- if (nt>1) pqr2(rbind(X,ldS$E),nt=nt) else qr(rbind(X,ldS$E),LAPACK=TRUE)
rp <- qrx$pivot;rp[rp] <- 1:np ## reverse pivot vector
## find pivoted \hat beta...
R <- qr.R(qrx)
......@@ -518,10 +522,12 @@ Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NU
## get component derivatives based on IFT...
dift <- Sl.ift(ldS$Sl,R,X,y,beta,qrx$pivot,rp)
## and the derivatives of log|X'X+S|...
P <- backsolve(R,diag(np))[rp,] ## invert R and row unpivot
PP <- tcrossprod(P) ## PP'
P <- pbsi(R,nt=nt,copy=TRUE) ## invert R
## P <- backsolve(R,diag(np))[rp,] ## invert R and row unpivot
## crossprod and unpivot (don't unpivot if unpivoting P above)
PP <- if (nt==1) tcrossprod(P)[rp,rp] else pRRt(P,nt)[rp,rp] ## PP'
ldetXXS <- 2*sum(log(abs(diag(R)))) ## log|X'X+S|
dXXS <- d.detXXS(ldS$Sl,PP) ## derivs of log|X'X+S|
dXXS <- d.detXXS(ldS$Sl,PP,nt=nt) ## derivs of log|X'X+S|
## all ingredients are now in place to form REML score and
## its derivatives....
reml <- (rss.bSb/phi + (nobs-Mp)*log(2*pi*phi) +
......@@ -547,7 +553,7 @@ Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NU
} ## Sl.fit
fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
rss.extra=0,nobs=NULL,Mp=0,conv.tol=.Machine$double.eps^.5) {
rss.extra=0,nobs=NULL,Mp=0,conv.tol=.Machine$double.eps^.5,nt=1) {
## estimates log smoothing parameters rho, by optimizing fast REML
## using Newton's method. On input Sl is a block diagonal penalty
## structure produced by Sl.setup, while X is a model matrix
......@@ -560,7 +566,7 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
if (is.null(nobs)) nobs <- nrow(X)
np <- ncol(X)
if (nrow(X) > np) { ## might as well do an initial QR step
qrx <- qr(X,LAPACK=TRUE)
qrx <- if (nt>1) pqr2(X,nt=nt) else qr(X,LAPACK=TRUE)
rp <- qrx$pivot
rp[rp] <- 1:np
X <- qr.R(qrx)[,rp]
......@@ -580,7 +586,7 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
fixed <- rep(FALSE,nrow(L))
best <- Sl.fit(Sl,X,y,L%*%rho+rho.0,fixed,log.phi,phi.fixed,rss.extra,nobs,Mp)
best <- Sl.fit(Sl,X,y,L%*%rho+rho.0,fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt)
## get a typical scale for the reml score...
reml.scale <- abs(best$reml) + best$rss/best$nobs
......@@ -624,12 +630,12 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
step[uconv.ind] <- uc.step ## step includes converged
## try out the step...
rho1 <- L%*%(rho + step)+rho.0; if (!phi.fixed) log.phi <- rho1[nr+1]
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp)
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt)
k <- 0
while (trial$reml>best$reml && k<35) { ## step half until improvement
step <- step/2;k <- k + 1
rho1 <- L%*%(rho + step)+rho.0; if (!phi.fixed) log.phi <- rho1[nr+1]
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp)
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt)
}
if (k==35 && trial$reml>best$reml) { ## step has failed
step.failed <- TRUE
......@@ -667,7 +673,7 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
best ## return the best fit (note that it will need post-processing to be useable)
} ## end fast.REML.fit
ident.test <- function(X,E) {
ident.test <- function(X,E,nt=1) {
## routine to identify structurally un-identifiable coefficients
## for model with model matrix X and scaled sqrt penalty matrix E
## lambda is smoothing parameter vector corresponding to E,
......@@ -679,7 +685,7 @@ ident.test <- function(X,E) {
## then beta.full[undrop] <- beta, is the full, zero padded
## coeff vector, with dropped coefs re-nstated as zeroes.
Xnorm <- norm(X,type="F")
qrx <- qr(rbind(X/Xnorm,E),LAPACK=TRUE) ## pivoted QR
qrx <- if (nt>1) pqr2(rbind(X/Xnorm,E),nt=nt) else qr(rbind(X/Xnorm,E),LAPACK=TRUE) ## pivoted QR
rank <- Rrank(qr.R(qrx),tol=.Machine$double.eps^.75)
drop <- qrx$pivot[-(1:rank)] ## index of un-identifiable coefs
undrop <- 1:ncol(X)
......@@ -745,13 +751,13 @@ Sl.drop <- function(Sl,drop,np) {
Sl
} ## Sl.drop
Sl.Xprep <- function(Sl,X) {
Sl.Xprep <- function(Sl,X,nt=1) {
## Sl is block diag object from Sl.setup, X is a model matrix
## this routine applies preliminary Sl transformations to X
## tests for structural identifibility problems and drops
## un-identifiable parameters.
X <- Sl.initial.repara(Sl,X) ## apply re-para used in Sl to X
id <- ident.test(X,attr(Sl,"E")) ## deal with structural identifiability
id <- ident.test(X,attr(Sl,"E"),nt=nt) ## deal with structural identifiability
## id contains drop, undrop, lambda
if (length(id$drop)>0) { ## then there is something to do here
Sl <- Sl.drop(Sl,id$drop,ncol(X)) ## drop unidentifiable from Sl
......
......@@ -1323,7 +1323,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
} else if (score1>=score) { ## initial step failed to improve score, try step halving ...
step <- Nstep ## start with the (pseudo) Newton direction
##sc.extra <- 1e-4*sum(grad*step) ## -ve sufficient decrease
while (score1>score && ii < maxHalf) {
while (score1>=score && ii < maxHalf) {
if (ii==3&&i<10) { ## Newton really not working - switch to SD, but keeping step length
s.length <- min(sum(step^2)^.5,maxSstep)
step <- Sstep*s.length/sum(Sstep^2)^.5 ## use steepest descent direction
......@@ -1349,7 +1349,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
score1 <- b1$UBRE
} else score1 <- b1$GCV
##sc.extra <- 1e-4*sum(grad*Nstep) ## -ve sufficient decrease
if (score1 <= score) { ## accept
if (score1 < score) { ## accept
if (pdef||!sd.unused) { ## then accept and compute derivatives
b <- gam.fit3(x=X, y=y, sp=L%*%lsp1+lsp0,Eb=Eb,UrS=UrS,
offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=2,
......@@ -1383,7 +1383,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
}
score1 <- score - abs(score) - 1 ## make sure that score1 < score
} # end of if (score1<= score ) # accept
if (score1>score) ii <- ii + 1
if (score1>=score) ii <- ii + 1
} ## end while (score1>score && ii < maxHalf)
if (!pdef&&sd.unused&&ii<maxHalf) score1 <- score2
} ## end of step halving branch
......
This diff is collapsed.
......@@ -15,7 +15,7 @@ Rrank <- function(R,tol=.Machine$double.eps^.9) {
rank
}
slanczos <- function(A,k=10,kl=-1,tol=.Machine$double.eps^.5) {
slanczos <- function(A,k=10,kl=-1,tol=.Machine$double.eps^.5,nt=1) {
## Computes truncated eigen decomposition of symmetric A by
## Lanczos iteration. If kl < 0 then k largest magnitude
## eigenvalues returned, otherwise k highest and kl lowest.
......@@ -36,7 +36,7 @@ slanczos <- function(A,k=10,kl=-1,tol=.Machine$double.eps^.5) {
if (n != ncol(A)) stop("A not square")
if (m>n) stop("Can not have more eigenvalues than nrow(A)")
oo <- .C(C_Rlanczos,A=as.double(A),U=as.double(rep(0,n*m)),D=as.double(rep(0,m)),
n=as.integer(n),m=as.integer(k),ml=as.integer(kl),tol=as.double(tol))
n=as.integer(n),m=as.integer(k),ml=as.integer(kl),tol=as.double(tol),nt=as.integer(nt))
list(values = oo$D,vectors = matrix(oo$U,n,m),iter=oo$n)
}
......@@ -708,7 +708,8 @@ gam.setup.list <- function(formula,pterms,
gam.setup <- function(formula,pterms,
data=stop("No data supplied to gam.setup"),knots=NULL,sp=NULL,
min.sp=NULL,H=NULL,absorb.cons=TRUE,sparse.cons=0,select=FALSE,idLinksBases=TRUE,
scale.penalty=TRUE,paraPen=NULL,gamm.call=FALSE,drop.intercept=FALSE)
scale.penalty=TRUE,paraPen=NULL,gamm.call=FALSE,drop.intercept=FALSE,
diagonal.penalty=FALSE)
## set up the model matrix, penalty matrices and auxilliary information about the smoothing bases
## needed for a gam fit.
## elements of returned object:
......@@ -846,12 +847,12 @@ gam.setup <- function(formula,pterms,
id <- split$smooth.spec[[i]]$id
if (is.null(id)||!idLinksBases) { ## regular evaluation
sml <- smoothCon(split$smooth.spec[[i]],data,knots,absorb.cons,scale.penalty=scale.penalty,
null.space.penalty=select,sparse.cons=sparse.cons)
null.space.penalty=select,sparse.cons=sparse.cons,diagonal.penalty=diagonal.penalty)
} else { ## it's a smooth with an id, so basis setup data differs from model matrix data
names(id.list[[id]]$data) <- split$smooth.spec[[i]]$term ## give basis data suitable names
sml <- smoothCon(split$smooth.spec[[i]],id.list[[id]]$data,knots,
absorb.cons,n=nrow(data),dataX=data,scale.penalty=scale.penalty,
null.space.penalty=select,sparse.cons=sparse.cons)
null.space.penalty=select,sparse.cons=sparse.cons,diagonal.penalty=diagonal.penalty)
}
for (j in 1:length(sml)) {
newm <- newm + 1
......@@ -2372,8 +2373,8 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
object$Vp <- object$Vc
}
if (type!="link"&&type!="terms"&&type!="iterms"&&type!="response"&&type!="lpmatrix")
{ warning("Unknown type, reset to terms.")
if (type!="link"&&type!="terms"&&type!="iterms"&&type!="response"&&type!="lpmatrix") {
warning("Unknown type, reset to terms.")
type<-"terms"
}
if (!inherits(object,"gam")) stop("predict.gam can only be used to predict from gam objects")
......@@ -2415,10 +2416,21 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
## get names of required variables, less response, but including offset variable
## see ?terms.object and ?terms for more information on terms objects
yname <- all.vars(object$terms)[attr(object$terms,"response")]
if (!is.null(newdata[[yname]])) { ## response provided...
naresp <- FALSE
if (!is.null(object$family$predict)&&!is.null(newdata[[yname]])) {
## response provided, and potentially needed for prediction (e.g. Cox PH family)
if (!is.null(object$pred.formula)) object$pred.formula <- attr(object$pred.formula,"full")
response <- TRUE
Terms <- terms(object)
resp <- newdata[[yname]]
if (sum(is.na(resp))>0) {
naresp <- TRUE ## there are NAs in supplied response
## replace them with a numeric code, so that rows are not dropped below
rar <- range(resp,na.rm=TRUE)
thresh <- rar[1]*1.01-rar[2]*.01
resp[is.na(resp)] <- thresh
newdata[[yname]] <- thresh
}
} else { ## response not provided
response <- FALSE
Terms <- delete.response(terms(object))
......@@ -2427,10 +2439,12 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
if (length(allNames) > 0) {
ff <- if (is.null(object$pred.formula)) reformulate(allNames) else object$pred.formula
if (sum(!(allNames%in%names(newdata)))) {
warning("not all required variables have been supplied in newdata!\n")}
warning("not all required variables have been supplied in newdata!\n")
}
## note that `xlev' argument not used here, otherwise `as.factor' in
## formula can cause a problem ... levels reset later.
newdata <- eval(model.frame(ff,data=newdata,na.action=na.act),parent.frame())
newdata <- eval(model.frame(ff,data=newdata,na.action=na.act),parent.frame())
if (naresp) newdata[[yname]][newdata[[yname]]<=thresh] <- NA ## reinstate as NA
} ## otherwise it's intercept only and newdata can be left alone
na.act <- attr(newdata,"na.action")
response <- if (response) newdata[[yname]] else NULL
......@@ -2651,16 +2665,7 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
se <- se[,order==1,drop=FALSE]
colnames(se) <- term.labels }
}
if (!is.null(terms)) { # return only terms requested via `terms'
if (sum(!(terms %in%colnames(fit))))
warning("non-existent terms requested - ignoring")
else {
fit <- fit[,terms,drop=FALSE]
if (se.fit) {## se <- as.matrix(as.matrix(se)[,terms])
se <- se[,terms,drop=FALSE]
}
}
}
} else { ## "link" or "response" case
fam <- object$family
k <- attr(attr(object$model,"terms"),"offset")
......@@ -2736,6 +2741,18 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
rm(X)
} ## end of prediction block loop
if ((type=="terms"||type=="iterms")&&!is.null(terms)) { # return only terms requested via `terms'
if (sum(!(terms %in%colnames(fit))))
warning("non-existent terms requested - ignoring")
else {
fit <- fit[,terms,drop=FALSE]
if (se.fit) {
se <- se[,terms,drop=FALSE]
}
}
}
if (type=="response"&&!is.null(fit1)) {
fit <- fit1
if (se.fit) se <- se1
......@@ -3953,8 +3970,9 @@ initial.sp <- function(X,S,off,expensive=FALSE)
maS <- max(abs(S[[i]]))
rsS <- rowMeans(abs(S[[i]]))
csS <- colMeans(abs(S[[i]]))
thresh <- .Machine$double.eps*maS*10
ind <- rsS > thresh & csS > thresh # only these columns really penalize
dS <- diag(abs(S[[i]])) ## new 1.8-4
thresh <- .Machine$double.eps^.8 * maS ## .Machine$double.eps*maS*10
ind <- rsS > thresh & csS > thresh & dS > thresh # only these columns really penalize
ss <- diag(S[[i]])[ind] # non-zero elements of l.d. S[[i]]
start <- off[i];finish <- start+ncol(S[[i]])-1
xx <- ldxx[start:finish]
......
......@@ -49,10 +49,10 @@ pinv <- function(X,svd=FALSE) {
X
} ## end pinv
pqr2 <- function(x,nt=1) {
pqr2 <- function(x,nt=1,nb=30) {
## Function for parallel pivoted qr decomposition of a matrix using LAPACK
## householder routines...
## library(mgcv); n <- 10000;p<-1000;x <- matrix(runif(n*p),n,p)
## householder routines. Currently uses a block algorithm.
## library(mgcv); n <- 10000;p<-500;x <- matrix(runif(n*p),n,p)
## system.time(qrx <- qr(x,LAPACK=TRUE))
## system.time(qrx2 <- mgcv:::pqr2(x,2))
## system.time(qrx3 <- mgcv:::pqr(x,2))
......@@ -60,12 +60,54 @@ pqr2 <- function(x,nt=1) {
p <- ncol(x)
beta <- rep(0.0,p)
piv <- as.integer(rep(0,p))
xc <- x*1
rank <- .Call(C_mgcv_Rpiqr,xc,beta,piv,nt)
ret <- list(qr=xc,rank=rank,qraux=beta,pivot=piv+1)
## need to force a copy of x, otherwise x will be over-written
## by .Call *in environment from which function is called*
x <- x*1
rank <- .Call(C_mgcv_Rpiqr,x,beta,piv,nt,nb)
ret <- list(qr=x,rank=rank,qraux=beta,pivot=piv+1)
attr(ret,"useLAPACK") <- TRUE
class(ret) <- "qr"
ret
} ## pqr2
pbsi <- function(R,nt=1,copy=TRUE) {
## parallel back substitution inversion of upper triangular R
## library(mgcv); n <- 5000;p<-4000;x <- matrix(runif(n*p),n,p)
## qrx <- mgcv:::pqr2(x,2);R <- qr.R(qrx)
## system.time(Ri <- mgcv:::pbsi(R,2))
## system.time(Ri2 <- backsolve(R,diag(p)));range(Ri-Ri2)
if (copy) R <- R * 1 ## ensure that R modified only within pbsi
.Call(C_mgcv_Rpbsi,R,nt)
R
} ## pbsi
pchol <- function(A,nt=1,nb=30) {
## parallel Choleski factorization.
## library(mgcv);
## set.seed(2);n <- 200;r <- 190;A <- tcrossprod(matrix(runif(n*r),n,r))
## system.time(R <- chol(A,pivot=TRUE));system.time(L <- mgcv:::pchol(A));range(R[1:r,]-L[1:r,])
## system.time(L <- mgcv:::pchol(A,nt=2,nb=30))
## piv <- attr(L,"pivot");attr(L,"rank");range(crossprod(L)-A[piv,piv])
## should nb be obtained from 'ILAENV' as page 23 of Lucas 2004??
piv <- as.integer(rep(0,ncol(A)))
A <- A*1 ## otherwise over-write in calling env!
rank <- .Call(C_mgcv_Rpchol,A,piv,nt,nb)
attr(A,"pivot") <- piv+1;attr(A,"rank") <- rank
A
}
pRRt <- function(R,nt=1) {
## parallel RR' for upper triangular R
## following creates index of lower triangular elements...
## n <- 4000;a <- rep(1:n,n);b <- rep(1:n,each=n)-1;which(a-b>0) -> ii;a[ii]+b[ii]*n->ii
## library(mgcv);R <- matrix(0,n,n);R[ii] <- runif(n*(n+1)/2)
## Note: A[a-b<=0] <- 0 zeroes upper triangle
## system.time(A <- mgcv:::pRRt(R,2))
## system.time(A2 <- tcrossprod(R));range(A-A2)
n <- nrow(R)
A <- matrix(0,n,n)
.Call(C_mgcv_RPPt,A,R,nt)
A
}
block.reorder <- function(x,n.blocks=1,reverse=FALSE) {
......@@ -82,7 +124,8 @@ block.reorder <- function(x,n.blocks=1,reverse=FALSE) {
oo <- .C(C_row_block_reorder,x=as.double(x),as.integer(r),as.integer(cols),
as.integer(nb),as.integer(reverse));
matrix(oo$x,r,cols)
}
} ## block.reorder
pqr <- function(x,nt=1) {
## parallel QR decomposition, using openMP in C, and up to nt threads (only if worthwhile)
......
......@@ -192,6 +192,7 @@ k.check <- function(b,subsample=5000,n.rep=400) {
nr <- length(rsd)
for (k in 1:m) { ## work through smooths
ok <- TRUE
b$smooth[[k]]$by <- "NA" ## can't deal with by variables
dat <- ExtractData(b$smooth[[k]],modf,NULL)$data
if (!is.null(attr(dat,"index"))||!is.null(attr(dat[[1]],"matrix"))||is.matrix(dat[[1]])) ok <- FALSE
if (ok) dat <- as.data.frame(dat)
......
This diff is collapsed.
......@@ -18,7 +18,8 @@ bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
na.action=na.omit, offset=NULL,method="fREML",control=list(),
scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL,paraPen=NULL,
chunk.size=10000,rho=0,AR.start=NULL,sparse=FALSE,cluster=NULL,
gc.level=1,use.chol=FALSE,samfrac=1,drop.unused.levels=TRUE,...)
nthreads=NA,gc.level=1,use.chol=FALSE,samfrac=1,
drop.unused.levels=TRUE,...)
}
%- maybe also `usage' for other objects documented here.
......@@ -41,7 +42,10 @@ covariates required by the formula. By default the variables are taken
from \code{environment(formula)}: typically the environment from
which \code{gam} is called.}
\item{weights}{ prior weights on the data.}
\item{weights}{ prior weights on the contribution of the data to the log likelihood. Note that a weight of 2, for example,