Commit 997f667a authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

New upstream version 1.8-28

parent d9ac53eb
......@@ -17,6 +17,18 @@ present, and reparameterization needs checking (also for bam).
* OpenBLAS 0.3.2/3 appears not to be thread safe at present - can
cause problems when opemmp multithreading.
1.8-28
* fix of obscure sp naming bug.
* changed some contour default colours from green to blue (they overlay
heatmaps, so green was not clever).
* Tweedie likelihood evaluation code made slightly more robust - for a model
with machine zero scale parameter estimate it could segfault, as series
maximum location could then overflow integer storage. Fixed + upper limit
imposed on series length (warning if it's not enough).
1.8-27
** Added routine 'ginla' for fully Bayesian inference based on an integrated
......@@ -25,7 +37,7 @@ present, and reparameterization needs checking (also for bam).
* Tweedie location scale family added: 'twlss'.
* gam.fit5 modified to distinguish more carefully between +ve semi definite
and +ve definite. Previously could fail, claiming indefinteness when it
and +ve definite. Previously could fail, claiming indefiniteness when it
should not have. Affects general families.
* bam was ignoring supplied scale parameter in extended family cases - fixed.
......@@ -71,7 +83,7 @@ present, and reparameterization needs checking (also for bam).
sarse matrices). Useful for other packages using constructors and
smooth2random, for 'fs' smooths.
* The mrf smooth constructor contained a obsolete hack in which the term
* The mrf smooth constructor contained an obsolete hack in which the term
dimension was set to 2 to avoid plotting when used as a te marginal. This
messed up side constraints for terms where a mrf smooth was a main effect
and te marginal. Fixed.
......
Package: mgcv
Version: 1.8-27
Version: 1.8-28
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with Automatic Smoothness
......@@ -20,6 +20,6 @@ LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2019-02-06 10:57:24 UTC; sw283
Packaged: 2019-03-21 08:56:32 UTC; sw283
Repository: CRAN
Date/Publication: 2019-02-06 15:00:03 UTC
Date/Publication: 2019-03-21 11:40:07 UTC
fd791b31993ade7c507f85afff3242e2 *ChangeLog
136bd8b6da7fe438d10fa69479e60c31 *DESCRIPTION
5fa8085364f960a26ac7690edfd6e423 *ChangeLog
1ca0566710cd19035e444ca30d94eda8 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
a7edea9aaf9e16a95b7844f9161515d6 *NAMESPACE
d35108707f4ccc3b5237626aa6dbe546 *R/bam.r
fac2c9f6811b7e760da5d55bac9b6270 *R/bam.r
2c132dc64eedf0c17d34f7061346f2d0 *R/coxph.r
5e3e5d8699554ed9a3593195055d99ba *R/efam.r
efcb37e8ceeba5da36cb1e3ac36a2bbe *R/fast-REML.r
4eaaa4960d578cd66a4201e6d0dc1d3a *R/gam.fit3.r
9ba6e588b9bbd50a3165fdfb8c16d23a *R/gam.fit4.r
36b00131cb4d9d98b45d8d4acb4e5953 *R/efam.r
b6953debd773ebbaab28afe97047e34c *R/fast-REML.r
81a226e30a5e69e0620cfbb6029a283c *R/gam.fit3.r
922b0786419a82a2127a94b062befc4d *R/gam.fit4.r
1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
021b3c1566fa81b6260c17b199b011c8 *R/gamlss.r
8e457764d72567f936b70d5811473f9e *R/gamlss.r
ce69039c15d4db9cfda5d1309af8b242 *R/gamm.r
f4c20e36b4f1d42c5bf902270ba3515a *R/inla.r
10facb791e4cfd123d183f05660119c6 *R/jagam.r
1a82805404b2027b46a124f020da30e2 *R/mgcv.r
05722ab28faea920a01ca48b07ff30ba *R/mgcv.r
97e439c8d1f8a0562d1383d4e4877fe2 *R/misc.r
b73809727c6eb8f7c7b0a20e9625fddc *R/mvam.r
eed19ceca4c4114023f12cda3f622659 *R/plots.r
1d3c2d5dae39bdb462ed75411829957d *R/plots.r
ae5ad4af1b30efb65b5deb05f6b8e778 *R/smooth.r
7398607a9ba7e85276bcf09df0d9344b *R/soap.r
bde1774ce7903cabc57b3196f8872ea8 *R/sparse.r
......@@ -44,7 +44,7 @@ fd0cfd64be579f9fbecdbb7f2b8ec1eb *man/Sl.initial.repara.Rd
60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd
69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd
8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd
4114840d3ed1b238b0430eed07fbd44c *man/bam.Rd
9e06549c2c3ca2e4945dfb31693512a4 *man/bam.Rd
ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd
cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd
745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd
......@@ -81,14 +81,14 @@ d828d474a4b069f9b9e3ebe5f05b70ec *man/gam.selection.Rd
b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
eb8648cc6b3c9374b899abd2f2b12f7b *man/gam2objective.Rd
ed77ce6e1b941625485706d7e307b816 *man/gamObject.Rd
a2593990d6a87f7b783d0435062dec02 *man/gamSim.Rd
89148f2dc12caff5073ac70c8873b33f *man/gamSim.Rd
e5d2541f32dab56972f58b0773eba50c *man/gamlss.etamu.Rd
c7f140d128d1d1d76909499900faf49e *man/gamlss.gH.Rd
623851672072e852e0c0ff03a96de0d6 *man/gamm.Rd
222535dd19201cfd929efd3305b13f43 *man/gaulss.Rd
398a5c12285401c1d37a8edb58780bc3 *man/get.var.Rd
a62dd487f34f476f7f830ed9d1bc58dc *man/gevlss.Rd
3158643d44410c71f15075fab0e516e0 *man/ginla.Rd
09dba82ee7459800056dbd0511788ebf *man/ginla.Rd
39b47f30a7ea45382d62ca1753d876a8 *man/identifiability.Rd
4f96476abbf9692f52030d3859580a31 *man/in.out.Rd
6c33026ebb458483d34a04548c05d664 *man/inSide.Rd
......@@ -98,7 +98,7 @@ a62dd487f34f476f7f830ed9d1bc58dc *man/gevlss.Rd
4bc4d96708fc6454ad81207299608d28 *man/jagam.Rd
d37e837db3089e3c0eb105669c04f0c8 *man/k.check.Rd
87d942b17bee05bb662270b894b183e6 *man/ldTweedie.Rd
c7073b72336b0385080f2f4a40694416 *man/ldetS.Rd
adc52f0e82ccbbe115bf15045b0a15a4 *man/ldetS.Rd
1b314907562e9236747ec1e344ebe491 *man/linear.functional.terms.Rd
ab35592c0a29d1574579811ea1c6ec39 *man/logLik.gam.Rd
b1c95a20afd6eb0825c00b46b8c3cbfa *man/ls.size.Rd
......@@ -142,9 +142,9 @@ c45d0a4edfa63ff362bd34195b3059ca *man/scat.Rd
898e7cc2def2ee234475e68d0b904b29 *man/sdiag.Rd
d54f4042e212fca7704cf8428bdaea38 *man/single.index.Rd
6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
8020154bd5c709d11f0e7cf043df886d *man/smooth.construct.Rd
40a5e35173c808dde102e2bfc492ac9d *man/smooth.construct.Rd
3f3e0cd76b77e207ee5a6ff89e5f7a9f *man/smooth.construct.ad.smooth.spec.Rd
fc7fba34b89fdf29f66744d1fdadda79 *man/smooth.construct.bs.smooth.spec.Rd
d2f3fb49a85066ef08d395a244d08a4d *man/smooth.construct.bs.smooth.spec.Rd
2f0463d1aca0b8798da6e681bd4c6e54 *man/smooth.construct.cr.smooth.spec.Rd
f5e6d0f5122f61c336827b3615482157 *man/smooth.construct.ds.smooth.spec.Rd
bd58515156b5e07c006e69f29aa830d1 *man/smooth.construct.fs.smooth.spec.Rd
......@@ -171,7 +171,7 @@ a6feff25ec8241bf5afb3d9fe219d26d *man/totalPenaltySpace.Rd
f22f1cee0ff2b70628846d1d0f8e9a66 *man/trichol.Rd
87e6b4437d00fab4fc814f4cefa3795c *man/trind.generator.Rd
6e975eef6b1214e0be93fc641f982f67 *man/twlss.Rd
96c48dd705710f639d76a0d0cc3fb128 *man/uniquecombs.Rd
bc350bfd3f4f8316d3b29b247292a16d *man/uniquecombs.Rd
a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
cbb69f16706da27b44fc62c92dcae161 *man/ziP.Rd
......@@ -195,12 +195,12 @@ e16c691700bbb44215597c6b0c7e6d2e *src/coxph.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
87f4add06d78bdaebf588eadbdbffbbf *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c
613c7bccb3d307824cf278b829190087 *src/mat.c
cd7abd2d8c808315207aeb566065efbb *src/mat.c
e4cef7f1753153fbab242d1c4d4f7e7f *src/matrix.c
de37b0972199b796654405efc007f25b *src/matrix.h
8df4b96961491d76989b50856237ee2d *src/mgcv.c
c7b7d5e90758bb042f206c2e9e93b129 *src/mgcv.h
8e3efe7a0d6b9619671cf2f7471f2112 *src/misc.c
0e9733504235e94635c6688ce6bb6bb1 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
563938b7bb6504ab10df5376c4360220 *src/qp.c
073a4b5b0bc6e869c5b35478c47facf1 *src/qp.h
......
## routines for very large dataset generalized additive modelling.
## (c) Simon N. Wood 2009-2017
## (c) Simon N. Wood 2009-2019
ls.size <- function(x) {
......@@ -249,7 +249,6 @@ discrete.mf <- function(gp,mf,names.pmf,m=NULL,full=TRUE) {
nmarg <- if (inherits(gp$smooth.spec[[i]],"tensor.smooth.spec")) length(gp$smooth.spec[[i]]$margin) else 1
maxj <- if (gp$smooth.spec[[i]]$by=="NA") nmarg else nmarg + 1
mi <- if (is.null(m)||length(m)==1) m else m[i]
# if (!is.null(gp$smooth.spec[[i]]$id)) stop("discretization can not handle smooth ids")
j <- 0
for (jj in 1:maxj) { ## loop through marginals
if (jj==1&&maxj!=nmarg) termi <- gp$smooth.spec[[i]]$by else {
......@@ -599,24 +598,15 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
w <- dd$EDeta2 * .5
z <- (eta-offset) - dd$Deta.EDeta2
}
#if (rho!=0) {
# ind <- which(w<0)
# if (length(ind)>0) { ## substitute Fisher weights
# w[ind] <- dd$EDeta2[ind] * .5
# z[ind] <- (eta[ind]-offset[ind]) - dd$Deta.EDeta2[ind]
# }
#}
good <- is.finite(z)&is.finite(w)
w[!good] <- 0 ## drop if !good
z[!good] <- 0 ## irrelevant
#dev <- sum(family$dev.resids(G$y,mu,G$w,theta))
} else { ## exponential family
mu.eta.val <- mu.eta(eta)
good <- mu.eta.val != 0
mu.eta.val[!good] <- .1 ## irrelvant as weight is zero
z <- (eta - offset) + (G$y - mu)/mu.eta.val
w <- (G$w * mu.eta.val^2)/variance(mu)
#dev <- sum(family$dev.resids(G$y,mu,G$w))
}
......@@ -655,7 +645,7 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
## block diagonal penalty object, Sl, set up before loop
if (iter==1) { ## need to get initial smoothing parameters
lambda.0 <- initial.sp(qrx$R,G$S,G$off,XX=TRUE) ## note that this uses the untrasformed X'X in qrx$R
lambda.0 <- initial.sp(qrx$R,G$S,G$off,XX=TRUE) ## note that this uses the untransformed X'X in qrx$R
## convert intial s.p.s to account for L
lsp0 <- log(lambda.0) ## initial s.p.
if (!is.null(G$L)) lsp0 <-
......
......@@ -935,7 +935,7 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) {
#if (vec) lsth1 <- LS[,c(4,2)]
LS <- colSums(LS)
#if (!vec) lsth1 <- c(LS[4],LS[2])
lsth1 <- c(LS[4],LS[2])
lsth1 <- c(LS[4],LS[2]) ## deriv w.r.t. p then log scale
lsth2 <- matrix(c(LS[5],LS[6],LS[6],LS[3]),2,2)
list(ls=LS[1],lsth1=lsth1,lsth2=lsth2)
}
......
......@@ -401,7 +401,7 @@ ldetSblock <- function(rS,rho,deriv=2,root=FALSE,nt=1) {
list(det = 2*sum(log(diag(R))+log(d[piv])),det1=dS1,det2=dS2,E=E)
} ## ldetSblock
ldetS <- function(Sl,rho,fixed,np,root=FALSE,repara=TRUE,nt=1) {
ldetS <- function(Sl,rho,fixed,np,root=FALSE,repara=TRUE,nt=1,deriv=2) {
## Get log generalized determinant of S stored blockwise in an Sl list.
## If repara=TRUE multi-term blocks will be re-parameterized using gam.reparam, and
## a re-parameterization object supplied in the returned object.
......@@ -445,13 +445,13 @@ ldetS <- function(Sl,rho,fixed,np,root=FALSE,repara=TRUE,nt=1) {
ind <- k.sp:(k.sp+m-1) ## index for smoothing parameters
## call gam.reparam to deal with this block
## in a stable way...
grp <- if (repara) gam.reparam(Sl[[b]]$rS,lsp=rho[ind],deriv=2) else
ldetSblock(Sl[[b]]$rS,rho[ind],deriv=2,root=root,nt=nt)
grp <- if (repara) gam.reparam(Sl[[b]]$rS,lsp=rho[ind],deriv=deriv) else
ldetSblock(Sl[[b]]$rS,rho[ind],deriv=deriv,root=root,nt=nt)
Sl[[b]]$lambda <- exp(rho[ind])
ldS <- ldS + grp$det
## next deal with the derivatives...
grp$det1 <- grp$det1[!fixed[ind]] ## discard derivatives for fixed components
grp$det2 <- grp$det2[!fixed[ind],!fixed[ind]]
grp$det2 <- if (deriv>1) grp$det2[!fixed[ind],!fixed[ind]] else 0 ##NULL
nd <- length(grp$det1)
if (nd>0) { ## then not all sp's are fixed
dind <- k.deriv:(k.deriv+nd-1)
......@@ -735,7 +735,7 @@ Sl.termMult <- function(Sl,A,full=FALSE,nt=1) {
SA
} ## end Sl.termMult
d.detXXS <- function(Sl,PP,nt=1) {
d.detXXS <- function(Sl,PP,nt=1,deriv=2) {
## 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,nt=nt) ## SPP[[k]] is S_k PP'
......@@ -744,11 +744,13 @@ d.detXXS <- function(Sl,PP,nt=1) {
for (i in 1:nd) {
indi <- attr(SPP[[i]],"ind")
d1[i] <- sum(diag(SPP[[i]][,indi,drop=FALSE]))
for (j in i:nd) {
indj <- attr(SPP[[j]],"ind")
d2[i,j] <- d2[j,i] <- -sum(t(SPP[[i]][,indj,drop=FALSE])*SPP[[j]][,indi,drop=FALSE])
}
d2[i,i] <- d2[i,i] + d1[i]
if (deriv==2) {
for (j in i:nd) {
indj <- attr(SPP[[j]],"ind")
d2[i,j] <- d2[j,i] <- -sum(t(SPP[[i]][,indj,drop=FALSE])*SPP[[j]][,indi,drop=FALSE])
}
d2[i,i] <- d2[i,i] + d1[i]
}
}
list(d1=d1,d2=d2)
} ## end d.detXXS
......@@ -803,10 +805,12 @@ Sl.iftChol <- function(Sl,XX,R,d,beta,piv,nt=1) {
## alternative all in one code - matches loop results, but
## timing close to identical - modified for parallel exec
D <- matrix(unlist(Skb),nrow(XX),nd)
D <- matrix(unlist(Skb),length(beta),nd)
bSb1 <- colSums(beta*D)
D1 <- .Call(C_mgcv_Rpforwardsolve,R,D[piv,]/d[piv],nt) ## note R transposed internally unlike forwardsolve
db[piv,] <- -.Call(C_mgcv_Rpbacksolve,R,D1,nt)/d[piv]
if (is.null(XX)) return(list(bSb1=bSb1,db=db)) ## return early
## XX.db <- XX%*%db
XX.db <- .Call(C_mgcv_pmmult2,XX,db,0,0,nt)
......@@ -825,7 +829,7 @@ Sl.fitChol <- function(Sl,XX,f,rho,yy=0,L=NULL,rho0=0,log.phi=0,phi.fixed=TRUE,
nobs=0,Mp=0,nt=c(1,1),tol=0,gamma=1) {
## given X'WX in XX and f=X'Wy solves the penalized least squares problem
## with penalty defined by Sl and rho, and evaluates a REML Newton step, the REML
## gradiant and the the estimated coefs bhat. If phi.fixed=FALSE then we need
## gradient and the the estimated coefs bhat. If phi.fixed=FALSE then we need
## yy = y'Wy in order to get derivsatives w.r.t. phi.
## NOTE: with an optimized BLAS nt==1 is likely to be much faster than
## nt > 1
......@@ -1198,7 +1202,7 @@ Sl.Xprep <- function(Sl,X,nt=1) {
}
rank <- 0
for (b in 1:length(Sl)) rank <- rank + Sl[[b]]$rank ## the total penalty rank
## Also add Mp, the total mull space dimension to return list.
## Also add Mp, the total null space dimension to return list.
list(X=X,Sl=Sl,undrop=id$undrop,rank=rank,Mp=ncol(X)-rank)
} ## end Sl.Xprep
......
......@@ -28,7 +28,7 @@ gam.reparam <- function(rS,lsp,deriv)
rSncol <- unlist(lapply(rS,ncol))
M <- length(lsp)
if (length(rS)>M) fixed.penalty <- TRUE else fixed.penalty <- FALSE
d.tol <- .Machine$double.eps^.3 ## group `similar sized' penalties, to save work
r.tol <- .Machine$double.eps^.75 ## This is a bit delicate -- too large and penalty range space can be supressed.
......@@ -1672,8 +1672,6 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
} ## newton
bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
control,gamma,scale,conv.tol=1e-6,maxNstep=3,maxSstep=2,
maxHalf=30,printWarn=FALSE,scoreType="GCV",start=NULL,
......@@ -1810,6 +1808,8 @@ bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
initial$dVkk <- diag(t(L0) %*% b$dVkk %*% L0)
initial$score <- score;initial$grad <- grad;
initial$scale.est <- b$scale.est
start0 <- coef(b)
mustart0 <- fitted(b)
rm(b)
B <- diag(length(initial$grad)) ## initial Hessian
......@@ -1819,7 +1819,7 @@ bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
b <- gam.fit3(x=X, y=y, sp=L%*%ilsp+lsp0,Eb=Eb,UrS=UrS,
offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=1,
control=control,gamma=gamma,scale=scale,printWarn=FALSE,
start=start,mustart=mustart,
start=start0,mustart=mustart0,
scoreType=scoreType,null.coef=null.coef,
pearson.extra=pearson.extra,dev.extra=dev.extra,n.true=n.true,Sl=Sl,...)
if (reml) {
......@@ -2497,7 +2497,7 @@ fix.family.var <- function(fam)
} ## fix.family.var
fix.family.ls<-function(fam)
fix.family.ls <- function(fam)
# adds ls the log saturated likelihood and its derivatives
# w.r.t. the scale parameter to the family object.
{ if (!inherits(fam,"family")) stop("fam not a family object")
......@@ -2951,6 +2951,12 @@ ldTweedie <- function(y,mu=y,p=1.5,phi=1,rho=NA,theta=NA,a=1.001,b=1.999,all.der
w2pp=as.double(y*0),y=as.double(y),eps=as.double(.Machine$double.eps^2),n=as.integer(length(y)),
th=as.double(theta[ind]),rho=as.double(rho[ind]),a=as.double(a),b=as.double(b))
}
if (oo$eps < -.5) {
if (oo$eps < -1.5) { ## failure of series in C code
oo$w2 <- oo$w1 <- oo$w2p <- oo$w1p <- oo$w2pp <- rep(NA,length(y))
}
else warning("Tweedie density may be unreliable - series not fully converged")
}
phii <- phi[ind]
if (!work.param) { ## transform working param derivatives to p/phi derivs...
if (length(dthp1)!=n) dthp1 <- array(dthp1,dim=n)
......@@ -2960,7 +2966,7 @@ ldTweedie <- function(y,mu=y,p=1.5,phi=1,rho=NA,theta=NA,a=1.001,b=1.999,all.der
oo$w1 <- oo$w1/phii
oo$w2p <- oo$w2p*dthp1i^2 + dthp2[ind] * oo$w1p
oo$w1p <- oo$w1p*dthp1i
oo$w2pp <- oo$w2pp*dthp1i/phii ## this appears to be wrong
oo$w2pp <- oo$w2pp*dthp1i/phii
}
......@@ -3014,7 +3020,7 @@ ldTweedie <- function(y,mu=y,p=1.5,phi=1,rho=NA,theta=NA,a=1.001,b=1.999,all.der
ld[ind,6] <- ld[ind,6] + oo$w2pp ## d2 logf / dphi dp
}
if (FALSE) { ## DEBUG disconnetion of density terms
if (FALSE) { ## DEBUG disconnection of density terms
ld[ind,1] <- oo$w ## log density
ld[ind,2] <- oo$w1 ## d log f / dphi
ld[ind,3] <- oo$w2 ## d2 logf / dphi2
......
......@@ -326,13 +326,6 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
## if (!is.null(start)) null.coef <- start - can be on edge of feasible - not good
null.eta <- as.numeric(x%*%null.coef + as.numeric(offset))
#old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef
#if (!is.null(start)) { ## check it's at least better than null.coef
# pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start
# if (pdev>old.pdev) start <- mustart <- etastart <- NULL
#}
## call the families initialization code...
if (is.null(mustart)) {
......@@ -372,7 +365,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
mu <- linkinv(eta);etaold <- eta
conv <- boundary <- FALSE
dd <- dDeta(y,mu,weights,theta,family,0) ## derivatives of deviance w.r.t. eta
w <- dd$Deta2 * .5;
w <- dd$Deta2 * .5
wz <- w*(eta-offset) - .5*dd$Deta
z <- (eta-offset) - dd$Deta.Deta2
good <- is.finite(z)&is.finite(w)
......@@ -907,7 +900,6 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
if (is.null(weights)) weights <- rep.int(1, nobs)
if (is.null(offset)) offset <- rep.int(0, nobs)
## get log likelihood, grad and Hessian (w.r.t. coefs - not s.p.s) ...
llf <- family$ll
ll <- llf(y,x,coef,weights,family,offset=offset,deriv=1)
......@@ -919,6 +911,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
check.deriv <- FALSE; eps <- 1e-5
drop <- NULL;bdrop <- rep(FALSE,q) ## by default nothing dropped
perturbed <- 0 ## counter for number of times perturbation tried on possible saddle
for (iter in 1:(2*control$maxit)) { ## main iteration
## get Newton step...
if (check.deriv) { ## code for checking derivatives when debugging
......
......@@ -195,144 +195,6 @@ gamlss.etamu <- function(l1,l2,l3=NULL,l4=NULL,ig1,g2,g3=NULL,g4=NULL,i2,i3=NULL
} # gamlss.etamu
gamlss.gH0 <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=NULL,D=NULL) {
## X[,jj[[i]]] is the ith model matrix.
## lj contains jth derivatives of the likelihood for each datum,
## columns are w.r.t. different combinations of parameters.
## ij is the symmetric array indexer for the jth order derivs...
## e.g. l4[,i4[i,j,l,m]] contains derivatives with
## respect to parameters indexed by i,j,l,m
## d1b and d2b are first and second derivatives of beta w.r.t. sps.
## fh is a factorization of the penalized hessian, while D contains the corresponding
## Diagonal pre-conditioning weights.
## deriv: 0 - just grad and Hess
## 1 - diagonal of first deriv of Hess wrt sps
## 2 - first deriv of Hess wrt sps
## 3 - everything.
K <- length(jj)
p <- ncol(X);n <- nrow(X)
trHid2H <- d1H <- d2H <- NULL ## defaults
## the gradient...
lb <- rep(0,p)
for (i in 1:K) { ## first derivative loop
lb[jj[[i]]] <- colSums(l1[,i]*X[,jj[[i]],drop=FALSE])
}
## the Hessian...
lbb <- matrix(0,p,p)
for (i in 1:K) for (j in i:K) {
lbb[jj[[i]],jj[[j]]] <- t(X[,jj[[i]],drop=FALSE])%*%(l2[,i2[i,j]]*X[,jj[[j]],drop=FALSE])
lbb[jj[[j]],jj[[i]]] <- t(lbb[jj[[i]],jj[[j]]])
}
if (deriv>0) {
## the first derivative of the Hessian, using d1b
## the first derivates of the coefficients wrt the sps
m <- ncol(d1b) ## number of smoothing parameters
## stack the derivatives of the various linear predictors on top
## of each other...
d1eta <- matrix(0,n*K,m)
ind <- 1:n
for (i in 1:K) {
d1eta[ind,] <- X[,jj[[i]],drop=FALSE]%*%d1b[jj[[i]],]
ind <- ind + n
}
}
if (deriv==1) {
d1H <- matrix(0,p,m) ## only store diagonals of d1H
for (l in 1:m) {
for (i in 1:K) {
v <- rep(0,n);ind <- 1:n
for (q in 1:K) {
v <- v + l3[,i3[i,i,q]] * d1eta[ind,l]
ind <- ind + n
}
d1H[jj[[i]],l] <- colSums(X[,jj[[i]],drop=FALSE]*(v*X[,jj[[i]],drop=FALSE]))
}
}
} ## if deriv==1
if (deriv>1) {
d1H <- list()
for (l in 1:m) {
d1H[[l]] <- matrix(0,p,p)
for (i in 1:K) for (j in i:K) {
v <- rep(0,n);ind <- 1:n
for (q in 1:K) {
v <- v + l3[,i3[i,j,q]] * d1eta[ind,l]
ind <- ind + n
}
## d1H[[l]][jj[[j]],jj[[i]]] <-
d1H[[l]][jj[[i]],jj[[j]]] <- t(X[,jj[[i]],drop=FALSE])%*%(v*X[,jj[[j]],drop=FALSE])
d1H[[l]][jj[[j]],jj[[i]]] <- t(d1H[[l]][jj[[i]],jj[[j]]])
}
}
} ## if deriv>1
if (deriv>2) {
## need tr(Hp^{-1} d^2H/drho_k drho_j)
## First form the expanded model matrix...
VX <- Xe <- matrix(0,K*n,ncol(X))
ind <- 1:n
for (i in 1:K) {
Xe[ind,jj[[i]]] <- X[,jj[[i]]]
ind <- ind + n
}
## Now form Hp^{-1} Xe'...
if (is.list(fh)) { ## then the supplied factor is an eigen-decomposition
d <- fh$values;d[d>0] <- 1/d[d>0];d[d<=0] <- 0
Xe <- t(D*((fh$vectors%*%(d*t(fh$vectors)))%*%(D*t(Xe))))
} else { ## the supplied factor is a choleski factor
ipiv <- piv <- attr(fh,"pivot");ipiv[piv] <- 1:p
Xe <- t(D*(backsolve(fh,forwardsolve(t(fh),(D*t(Xe))[piv,]))[ipiv,]))
}
## now compute the required trace terms
d2eta <- matrix(0,n*K,ncol(d2b))
ind <- 1:n
for (i in 1:K) {
d2eta[ind,] <- X[,jj[[i]],drop=FALSE]%*%d2b[jj[[i]],]
ind <- ind + n
}
trHid2H <- rep(0,ncol(d2b))
kk <- 0 ## counter for second derivatives
for (k in 1:m) for (l in k:m) { ## looping over smoothing parameters...
kk <- kk + 1
for (i in 1:K) for (j in 1:K) {
v <- rep(0,n);ind <- 1:n
for (q in 1:K) { ## accumulate the diagonal matrix for X_i'diag(v)X_j
v <- v + d2eta[ind,kk]*l3[,i3[i,j,q]]
ins <- 1:n
for (s in 1:K) {
v <- v + d1eta[ind,k]*d1eta[ins,l]*l4[,i4[i,j,q,s]]
ins <- ins + n
}
ind <- ind + n
}
if (i==j) {
rind <- 1:n + (i-1)*n
VX[rind,jj[[i]]] <- v * X[,jj[[i]]]
} else {
rind1 <- 1:n + (i-1)*n
rind2 <- 1:n + (j-1)*n
VX[rind2,jj[[i]]] <- v * X[,jj[[i]]]
VX[rind1,jj[[j]]] <- v * X[,jj[[j]]]
}
}
trHid2H[kk] <- sum(Xe*VX)
}
} ## if deriv>2
list(lb=lb, ## grad w.r.t. coefs
lbb=lbb, ## hess w.r.t. coefs, H
d1H=d1H, ## grad of H wrt sps
## d2H=d2H,
trHid2H=trHid2H) ## tr(H^{-1}d^2H/dspsp)
} ## end of gamlss.gH0
gamlss.gH <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=NULL,D=NULL) {
## X[,jj[[i]]] is the ith model matrix.
## lj contains jth derivatives of the likelihood for each datum,
......
......@@ -1080,7 +1080,7 @@ gam.setup <- function(formula,pterms,
# idea here is that terms are set up in accordance with information given in split$smooth.spec
# appropriate basis constructor is called depending on the class of the smooth
# constructor returns penalty matrices model matrix and basis specific information
## sm[[i]] <- smoothCon(split$smooth.spec[[i]],data,knots,absorb.cons,scale.penalty=scale.penalty,sparse.cons=sparse.cons) ## old code
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,
......@@ -1126,11 +1126,12 @@ gam.setup <- function(formula,pterms,
if (is.null(sm[[i]]$L)) Li <- diag(length.S) else Li <- sm[[i]]$L
if (length.S > 0) { ## there are smoothing parameters to name
if (length.S == 1) spn <- sm[[i]]$label else {
if (length.S == 1) lspn <- sm[[i]]$label else {
Sname <- names(sm[[i]]$S)
if (is.null(Sname)) spn <- paste(sm[[i]]$label,1:length.S,sep="") else
spn <- paste(sm[[i]]$label,Sname,sep="")
lspn <- if (is.null(Sname)) paste(sm[[i]]$label,1:length.S,sep="") else
paste(sm[[i]]$label,Sname,sep="") ## names for all sp's
}
spn <- lspn[1:ncol(Li)] ## names for actual working sps
}
## extend the global L matrix...
......@@ -1143,7 +1144,7 @@ gam.setup <- function(formula,pterms,
cbind(matrix(0,nrow(Li),ncol(L)),Li))
if (length.S > 0) { ## there are smoothing parameters to name
sp.names <- c(sp.names,spn) ## extend the sp name vector
lsp.names <- c(lsp.names,spn) ## extend full.sp name vector
lsp.names <- c(lsp.names,lspn) ## extend full.sp name vector
}
} else { ## it's a repeat id => shares existing sp's
L0 <- matrix(0,nrow(Li),ncol(L))
......@@ -1153,7 +1154,7 @@ gam.setup <- function(formula,pterms,
L0[,idx[[id]]$c:(idx[[id]]$c+ncol(Li)-1)] <- Li
L <- rbind(L,L0)
if (length.S > 0) { ## there are smoothing parameters to name
lsp.names <- c(lsp.names,spn) ## extend full.sp name vector
lsp.names <- c(lsp.names,lspn) ## extend full.sp name vector
}
}
}
......
......@@ -32,6 +32,10 @@ fix.family.qf <- function(fam) {
}
} else if (family=="binomial") {
fam$qf <- function(p,mu,wt,scale) {
if (all.equal(wt,ceiling(wt))!=TRUE) {
wt <- ceiling(wt)
warning("non-integer binomial denominator: quantiles incorrect")
}
qbinom(p,wt,mu)/(wt + as.numeric(wt==0))
}
} else if (family=="Gamma") {
......@@ -440,7 +444,7 @@ plot.sos.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=2,se2.mult=1,
pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,hcolors=heat.colors(100),
contour.col=3,...) {
contour.col=4,...) {
## plot method function for sos.smooth terms
if (scheme>1) return(plot.mgcv.smooth(x,P=P,data=data,label=label,se1.mult=se1.mult,se2.mult=se2.mult,
partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,
......@@ -702,7 +706,7 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=2,se2.mult=1,
pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,hcolors=heat.colors(50),
contour.col=3,...) {
contour.col=4,...) {
## default plot method for smooth objects `x' inheriting from "mgcv.smooth"
## `x' is a smooth object, usually part of a `gam' fit. It has an attribute
## 'coefficients' containing the coefs for the smooth, but usually these
......@@ -1551,7 +1555,7 @@ vis.gam <- function(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,color="he
surf.col<-surf.col/(max.z-min.z)
surf.col<-round(surf.col*nCol)
con.col <-1
if (color=="heat") { pal<-heat.colors(nCol);con.col<-3;}
if (color=="heat") { pal<-heat.colors(nCol);con.col<-4;}
else if (color=="topo") { pal<-topo.colors(nCol);con.col<-2;}
else if (color=="cm") { pal<-cm.colors(nCol);con.col<-1;}
else if (color=="terrain") { pal<-terrain.colors(nCol);con.col<-2;}
......
......@@ -13,7 +13,7 @@ are designed for datasets containing upwards of several tens of thousands of dat
of \code{bam} is much lower memory footprint than \code{\link{gam}}, but it can also be much faster,
for large datasets. \code{bam} can also compute on a cluster set up by the \link[parallel]{parallel} package.
An alternative fitting approach (Wood et al. 2017) is provided by the \code{discrete==TRUE} method. In this case a method based on discretization of covariate values and C code level parallelization (controlled by the \code{nthreads} argument instead of the \code{cluster} argument) is used. This extends both the data set and model size that are practical.
An alternative fitting approach (Wood et al. 2017, Li and Wood, 2019) is provided by the \code{discrete==TRUE} method. In this case a method based on discretization of covariate values and C code level parallelization (controlled by the \code{nthreads} argument instead of the \code{cluster} argument) is used. This extends both the data set and model size that are practical.
}
\usage{
......@@ -187,7 +187,7 @@ can be used in the 1D case, and tensor product smooths (\code{te}) are typically
If \code{cluster} is provided as a cluster set up using \code{\link[parallel]{makeCluster}} (or \code{\link[parallel]{makeForkCluster}}) from the \code{parallel} package, then the rate limiting QR decomposition of the model matrix is performed in parallel using this cluster. Note that the speed ups are often not that great. On a multi-core machine it is usually best to set the cluster size to the number of physical cores, which is often less than what is reported by \code{\link[parallel]{detectCores}}. Using more than the number of physical cores can result in no speed up at all (or even a slow down). Note that a highly parallel BLAS may negate all advantage from using a cluster of cores. Computing in parallel of course requires more memory than computing in series. See examples.
When \code{discrete=TRUE} the covariate data are first discretized. Discretization takes place on a smooth by smooth basis, or in the case of tensor product smooths (or any smooth that can be represented as such, such as random effects), separately for each marginal smooth. The required spline bases are then evaluated at the discrete values, and stored, along with index vectors indicating which original observation they relate to. Fitting is by a version of performance oriented iteration/PQL using REML smoothing parameter selection on each iterative working model (as for the default method). The iteration is based on the derivatives of the REML score, without computing the score itself, allowing the expensive computations to be reduced to one parallel block Cholesky decomposition per iteration (plus two basic operations of equal cost, but easily parallelized). Unlike standard POI/PQL, only one step of the smoothing parameter update for the working model is taken at each step (rather than iterating to the optimal set of smoothing parameters for each working model). At each step a weighted model matrix crossproduct of the model matrix is required - this is efficiently computed from the pre-computed basis functions evaluated at the discretized covariate values. Efficient computation with tensor product terms means that some terms within a tensor product may be re-ordered for maximum efficiency.
When \code{discrete=TRUE} the covariate data are first discretized. Discretization takes place on a smooth by smooth basis, or in the case of tensor product smooths (or any smooth that can be represented as such, such as random effects), separately for each marginal smooth. The required spline bases are then evaluated at the discrete values, and stored, along with index vectors indicating which original observation they relate to. Fitting is by a version of performance oriented iteration/PQL using REML smoothing parameter selection on each iterative working model (as for the default method). The iteration is based on the derivatives of the REML score, without computing the score itself, allowing the expensive computations to be reduced to one parallel block Cholesky decomposition per iteration (plus two basic operations of equal cost, but easily parallelized). Unlike standard POI/PQL, only one step of the smoothing parameter update for the working model is taken at each step (rather than iterating to the optimal set of smoothing parameters for each working model). At each step a weighted model matrix crossproduct of the model matrix is required - this is efficiently computed from the pre-computed basis functions evaluated at the discretized covariate values. Efficient computation with tensor product terms means that some terms within a tensor product may be re-ordered for maximum efficiency. See Wood et al (2017) and Li and Wood (2019) for full details.
When \code{discrete=TRUE} parallel computation is controlled using the \code{nthreads} argument. For this method no cluster computation is used, and the \code{parallel} package is not required. Note that actual speed up from parallelization depends on the BLAS installed and your hardware. With the (R default) reference BLAS using several threads can make a substantial difference, but with a single threaded tuned BLAS, such as openblas, the effect is less marked (since cache use is typically optimized for one thread, and is then sub optimal for several). However the tuned BLAS is usually much faster than using the reference BLAS, however many threads you use. If you have a multi-threaded BLAS installed then you should leave \code{nthreads} at 1, since calling a multi-threaded BLAS from multiple threads usually slows things down: the only exception to this is that you might choose to form discrete matrix cross products (the main cost in the fitting routine) in a multi-threaded way, but use single threaded code for other computations: this can be achieved by e.g. \code{nthreads=c(2,1)}, which would use 2 threads for discrete inner products, and 1 for most code calling BLAS. Not that the basic reason that multi-threaded performance is often disappointing is that most computers are heavily memory bandwidth limited, not flop rate limited. It is hard to get data to one core fast enough, let alone trying to get data simultaneously to several cores.
......@@ -204,6 +204,10 @@ Wood, S.N., Goude, Y. & Shaw S. (2015) Generalized additive models for large dat
Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association. 112(519):1199-1210
\url{http://dx.doi.org/10.1080/01621459.2016.1195744}
Li, Z & S.N. Wood (2019) Faster model matrix crossproducts for large generalized linear models with discretized covariates. Statistics and Computing.
\url{https://doi.org/10.1007/s11222-019-09864-2}
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}
......
......@@ -16,7 +16,7 @@ the length of the example code sections.
\item{n}{ number of data to simulate.}
\item{dist}{character string which may be used to spcify the distribution of
\item{dist}{character string which may be used to specify the distribution of
the response.}
\item{scale}{Used to set noise level.}
......
......@@ -2,7 +2,7 @@
\alias{ginla}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{GAM Integrated Nested Laplace Approximation Newton Enhanced}
\description{Apply Integrated Nested Laplace Approximation (INLA, Rue et al. 2009) to models estimable by \code{\link{gam}} or \code{\link{bam}}, using the INLA variant described in Wood (2018). Produces marginal posterior densities for each coefficient, selected coefficients or linear transformations of the coefficient vector.
\description{Apply Integrated Nested Laplace Approximation (INLA, Rue et al. 2009) to models estimable by \code{\link{gam}} or \code{\link{bam}}, using the INLA variant described in Wood (2019). Produces marginal posterior densities for each coefficient, selected coefficients or linear transformations of the coefficient vector.
}
\usage{
ginla(G,A=NULL,nk=16,nb=100,J=1,interactive=FALSE,int=0,approx=0)
......@@ -36,7 +36,7 @@ Note that for many models the INLA estimates are very close to the usual Gaussia
\references{
Rue, H, Martino, S. & Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series B. 71: 319-392.
Wood (2018, submitted) Simplified Integrated Laplace Approximation.
Wood (2019) Simplified Integrated Laplace Approximation. In press Biometrika.
}
\section{WARNINGS}{
......
......@@ -4,7 +4,7 @@
\alias{ldetS}
\title{Getting log generalized determinant of penalty matrices}
\usage{
ldetS(Sl, rho, fixed, np, root = FALSE, repara = TRUE, nt = 1)
ldetS(Sl, rho, fixed, np, root = FALSE, repara = TRUE, nt = 1,deriv=2)
}
\arguments{
\item{Sl}{the output of \code{Sl.setup}.}
......@@ -21,6 +21,8 @@ ldetS(Sl, rho, fixed, np, root = FALSE, repara = TRUE, nt = 1)
a re-parameterization object supplied in the returned object.}