Commit d9ac53eb authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

New upstream version 1.8-27

parent 93a4e91b
......@@ -17,11 +17,44 @@ 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-27
** Added routine 'ginla' for fully Bayesian inference based on an integrated
nested Laplace approximation (INLA) approach. See ?ginla.
* 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
should not have. Affects general families.
* bam was ignoring supplied scale parameter in extended family cases - fixed.
* work around in list formula handling for reformulate sometimes setting
response to a name in place of a call.
* preinitialize in general families is now a function, not an expression. See
cox.ph for an example.
* added routine cholup for rank one modification of Cholesky factor.
* two stage gam/bam fitting now allows 'sp' to be modified.
* predict.gam could fail with type="response" for families requiring the
response to be provided in this case (e.g. cox.ph). Fixed.
* sp.vcov defaults to extracting edge corrected log sp cov matrix, if
gam(...,gam.control(edge.control=TRUE)) used for fitting.
* gam(...,gam.control(edge.correct=TRUE)) could go into infinite loop if
sp was effectively zero. Corrected.
1.8-26
* LINPACK dependency removed.
* Added service routine chol.down to down date a Cholesky factor on row/col
* Added service routine choldrop to down date a Cholesky factor on row/col
deletion.
* liu2 had a check messed up when vectorized. Fix to stop vector being
......
Package: mgcv
Version: 1.8-26
Version: 1.8-27
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with Automatic Smoothness
......@@ -7,17 +7,19 @@ Title: Mixed GAM Computation Vehicle with Automatic Smoothness
Description: Generalized additive (mixed) models, some of their extensions and
other generalized ridge regression with multiple smoothing
parameter estimation by (Restricted) Marginal Likelihood,
Generalized Cross Validation and similar. Includes a gam()
function, a wide variety of smoothers, JAGS support and
distributions beyond the exponential family.
Generalized Cross Validation and similar, or using iterated
nested Laplace approximation for fully Bayesian inference. See
Wood (2017) <doi:10.1201/9781315370279> for an overview.
Includes a gam() function, a wide variety of smoothers, 'JAGS'
support and distributions beyond the exponential family.
Priority: recommended
Depends: R (>= 2.14.0), nlme (>= 3.1-64)
Imports: methods, stats, graphics, Matrix
Suggests: splines, parallel, survival, MASS
Imports: methods, stats, graphics, Matrix, splines, utils
Suggests: parallel, survival, MASS
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2018-11-20 20:40:48 UTC; sw283
Packaged: 2019-02-06 10:57:24 UTC; sw283
Repository: CRAN
Date/Publication: 2018-11-21 11:50:02 UTC
Date/Publication: 2019-02-06 15:00:03 UTC
27ebfd695d348a41dde36452454d4f80 *ChangeLog
8c7c5ba5650690f9409b109c928718b1 *DESCRIPTION
fd791b31993ade7c507f85afff3242e2 *ChangeLog
136bd8b6da7fe438d10fa69479e60c31 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
6fbddaad9c9fad126ce709cbd09e56d3 *NAMESPACE
459045bffc77934a4f287258f31f8aa6 *R/bam.r
de69a45fe00e28a0309193df283c5829 *R/coxph.r
6c79693fe31558339cd694883d0beeb1 *R/efam.r
8fb9318d7a272467924c1ef4cd4eb97e *R/fast-REML.r
a00617b7da9f9de5fff27f13b5086c76 *R/gam.fit3.r
5ae9e2ab9b0b30811a14e6be88f4588c *R/gam.fit4.r
a7edea9aaf9e16a95b7844f9161515d6 *NAMESPACE
d35108707f4ccc3b5237626aa6dbe546 *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
1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
f2368e22ccf048cf3d1e210e970013dd *R/gamlss.r
021b3c1566fa81b6260c17b199b011c8 *R/gamlss.r
ce69039c15d4db9cfda5d1309af8b242 *R/gamm.r
f4c20e36b4f1d42c5bf902270ba3515a *R/inla.r
10facb791e4cfd123d183f05660119c6 *R/jagam.r
081f37414ae038f5608a1a15e34bb34a *R/mgcv.r
a186527989eee0303f966f80a5fd7b84 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r
1a82805404b2027b46a124f020da30e2 *R/mgcv.r
97e439c8d1f8a0562d1383d4e4877fe2 *R/misc.r
b73809727c6eb8f7c7b0a20e9625fddc *R/mvam.r
eed19ceca4c4114023f12cda3f622659 *R/plots.r
530e432792d072e2587cbd0896b620cf *R/smooth.r
ae5ad4af1b30efb65b5deb05f6b8e778 *R/smooth.r
7398607a9ba7e85276bcf09df0d9344b *R/soap.r
bde1774ce7903cabc57b3196f8872ea8 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
......@@ -33,6 +34,7 @@ e43fa21d71e3a38fdf798e720aa8a7d1 *inst/po/de/LC_MESSAGES/R-mgcv.mo
fdc7570ff0761a731464fad88b23c062 *inst/po/pl/LC_MESSAGES/R-mgcv.mo
8d3d346c7f6d142b940e9ab751e1db71 *inst/po/pl/LC_MESSAGES/mgcv.mo
1f716ff54bb49f75188f92b830b49428 *man/Beta.Rd
e7c91ae36c5d1ec720d210fd23140991 *man/FFdes.Rd
9d12d96a7f2dfe7faa1c112bf6631a8e *man/Predict.matrix.Rd
9c2573f0134dc2cbac7e11998e220558 *man/Predict.matrix.cr.smooth.Rd
d0fa291cbbcef359c61e426f6ba38fbb *man/Predict.matrix.soap.film.Rd
......@@ -42,17 +44,17 @@ fd0cfd64be579f9fbecdbb7f2b8ec1eb *man/Sl.initial.repara.Rd
60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd
69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd
8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd
ca103588880287accfbde4265cc56f54 *man/bam.Rd
4114840d3ed1b238b0430eed07fbd44c *man/bam.Rd
ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd
cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd
745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd
530b8b4bacffa9353561f19ecfccfe19 *man/cSplineDes.Rd
05d2b5d4b20621a356a4dff2f88c3fb5 *man/chol.down.Rd
8133260bd3c26231322cd7cfbfa3b421 *man/chol.down.Rd
1e4a88fca144eda0f985b362ff5b3972 *man/choose.k.Rd
c03748964ef606621418e428ae49b103 *man/columb.Rd
9906a83ce29a3b17044bc2a3c9940cee *man/concurvity.Rd
04d4285d6b0cea7f8b3f4cb7ec9a466c *man/coxph.Rd
9a1cabe4d6abc94326a04e90c218bc93 *man/coxpht.Rd
13d652273fb27fa3ecb9bb704ec2da02 *man/coxph.Rd
239e4c9f917ff2d94d02972fa6c31e4d *man/coxpht.Rd
b78faf4ab9477183e7a3fbbd8801afeb *man/dDeta.Rd
0a6d4b3858cbb69a5d375ecd09282256 *man/exclude.too.far.Rd
3add7e72948f304246a412b3a91b9fad *man/extract.lme.cov.Rd
......@@ -63,7 +65,7 @@ e75719779e18c723ee1fd17e44e7901b *man/formXtViX.Rd
88888e966394c9f9792d9537341d053c *man/formula.gam.Rd
4da4d585b329769eb44f0c7a6e7dd554 *man/fs.test.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
22ddcbad33f08332cd15a6cb88b0a510 *man/gam.Rd
e7bd0528ea893d692736c360ad9e6436 *man/gam.Rd
42d669d0f18eba79234b2e849d724863 *man/gam.check.Rd
6d1a075aab054f1e14e6f0f1f432626a *man/gam.control.Rd
afd2fdc49ac57b4b94ee233fa5da1d64 *man/gam.convergence.Rd
......@@ -86,6 +88,7 @@ c7f140d128d1d1d76909499900faf49e *man/gamlss.gH.Rd
222535dd19201cfd929efd3305b13f43 *man/gaulss.Rd
398a5c12285401c1d37a8edb58780bc3 *man/get.var.Rd
a62dd487f34f476f7f830ed9d1bc58dc *man/gevlss.Rd
3158643d44410c71f15075fab0e516e0 *man/ginla.Rd
39b47f30a7ea45382d62ca1753d876a8 *man/identifiability.Rd
4f96476abbf9692f52030d3859580a31 *man/in.out.Rd
6c33026ebb458483d34a04548c05d664 *man/inSide.Rd
......@@ -102,7 +105,7 @@ b1c95a20afd6eb0825c00b46b8c3cbfa *man/ls.size.Rd
9a2c8f14c7a56eca69f4a59bef27a9bf *man/magic.Rd
5169af4be5fccf9fa79b6de08e9ea035 *man/magic.post.proc.Rd
5c4061016779a1504b6721cdf5182074 *man/mgcv-FAQ.Rd
894a7c1420d43c4dd32e0500809e28a4 *man/mgcv-package.Rd
1efa4f05a9b01dce18e0f75bca778259 *man/mgcv-package.Rd
e85bd20fe881288da19e187c4c5d7a16 *man/mgcv-parallel.Rd
1f085fc302c7c0b10a44d59c439c3ec5 *man/mini.roots.Rd
8a89166312f624d58d79c1f86ae43e30 *man/missing.data.Rd
......@@ -123,7 +126,7 @@ ee9352ba4c531a8def16deddcab9a9fd *man/pdIdnot.Rd
8bc429d92aa9f58c4c43f2852e1f8123 *man/pdTens.Rd
1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd
931c3aefb8b5e42aa230cfedd281bed1 *man/place.knots.Rd
8280e8aff6d5a004449d9e19661e3920 *man/plot.gam.Rd
1724c9c2f4ded94d4b855d474673e72a *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
1a9d83c9fc67e5f0fc85d66d3112f4ef *man/predict.bam.Rd
650cd7383fe7d54920e2e3273745780e *man/predict.gam.Rd
......@@ -157,7 +160,7 @@ c522c270c217e5b83cf8f3e95220a03f *man/smooth.construct.tp.smooth.spec.Rd
ae5e27524e37d57505754639455f18a5 *man/smooth.terms.Rd
f642b1caecf3d2bcdbbc0a884e1d3fa5 *man/smooth2random.Rd
844f9653d74441293d05a24dd3e2876a *man/smoothCon.Rd
b55a396da77559dac553613146633f97 *man/sp.vcov.Rd
08a186c579a81d75d9855e53743db03e *man/sp.vcov.Rd
83bd8e097711bf5bd0fff09822743d43 *man/spasm.construct.Rd
b9394812e5398ec95787c65c1325a027 *man/step.gam.Rd
62a42d898c2f1ccd7a64aef33d07b3a1 *man/summary.gam.Rd
......@@ -167,6 +170,7 @@ d1358e5c7f1f9d9a56072a77787803d2 *man/te.Rd
a6feff25ec8241bf5afb3d9fe219d26d *man/totalPenaltySpace.Rd
f22f1cee0ff2b70628846d1d0f8e9a66 *man/trichol.Rd
87e6b4437d00fab4fc814f4cefa3795c *man/trind.generator.Rd
6e975eef6b1214e0be93fc641f982f67 *man/twlss.Rd
96c48dd705710f639d76a0d0cc3fb128 *man/uniquecombs.Rd
a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
......@@ -185,17 +189,17 @@ ccbcb019a355ee59fcae8fa82b12bf2d *po/R-ko.po
b85b3ef1a39c4e1346b93173be8914aa *po/mgcv.pot
fdbb7d250e7baa7d73f26837246480f3 *po/pl.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
031c78f98f40e2615bf4794cb8fbae91 *src/coxph.c
e16c691700bbb44215597c6b0c7e6d2e *src/coxph.c
4801447cdbc0638d0aaf8ff723ba4b6b *src/discrete.c
0aa17445bd208f941d79e5dd70181297 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
818b7134dc51a935b07f1f2792757635 *src/init.c
87f4add06d78bdaebf588eadbdbffbbf *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c
6c0695593ed52375475c9bd15493d38a *src/mat.c
613c7bccb3d307824cf278b829190087 *src/mat.c
e4cef7f1753153fbab242d1c4d4f7e7f *src/matrix.c
de37b0972199b796654405efc007f25b *src/matrix.h
8df4b96961491d76989b50856237ee2d *src/mgcv.c
3cc1f17f3c2f04e0ba21639453959038 *src/mgcv.h
c7b7d5e90758bb042f206c2e9e93b129 *src/mgcv.h
8e3efe7a0d6b9619671cf2f7471f2112 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
563938b7bb6504ab10df5376c4360220 *src/qp.c
......@@ -205,7 +209,7 @@ d5673b88f6f3d85c62a1337f49abba24 *src/soap.c
fe0444bece322bc229e46b3d1c150779 *src/tprs.c
5bd85bf0319a7b7c755cf49c91a7cd94 *src/tprs.h
38e593a85a6fd0bb4fbed836f3361406 *tests/bam.R
3677e01cc0f90ae5ef9cc44583b045c8 *tests/coxpht.R
98029c1a07ed9b4b00be53cdd67d154a *tests/coxpht.R
fefd6fe58a089c4692652bc5c0bcc65c *tests/gam.R
fa2508c443bdc759a734df0d00ed735e *tests/mgcv-parallel.R
9c1a01e6ea9ce8855f5489bc67762ecb *tests/missing.data.R
......
useDynLib(mgcv, .registration = TRUE, .fixes = "C_")
export(anova.gam, bam, bam.update,bandchol, betar,
choldrop,cox.ph,concurvity,
choldrop,cholup,cox.ph,concurvity,
cSplineDes,dDeta,
exclude.too.far,extract.lme.cov, extract.lme.cov2,
exclude.too.far,extract.lme.cov, extract.lme.cov2,FFdes,
formXtViX, full.score, formula.gam,fixDependence,fix.family.link,
fix.family.var, fix.family.ls, fix.family.qf,fix.family.rd,
fs.test,fs.boundary,gam, gam2derivative,
......@@ -12,7 +12,7 @@ export(anova.gam, bam, bam.update,bandchol, betar,
gam.fit,gam.fit5.post.proc,
gamlss.etamu,gamlss.gH,
gam.outer,gam.reparam, gam.vcomp, gamSim ,
gaulss,gam.side,get.var,gevlss,
gaulss,gam.side,get.var,gevlss,ginla,
influence.gam,
in.out,inSide,interpret.gam,initial.sp,
jagam,k.check,ldetS,
......@@ -72,13 +72,16 @@ export(anova.gam, bam, bam.update,bandchol, betar,
spasm.construct,spasm.sp,spasm.smooth,
t2,te,ti,tensor.prod.model.matrix,tensor.prod.penalties,
totalPenaltySpace,trichol,trind.generator,
Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
Tweedie,tw,twlss,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
importFrom(grDevices,cm.colors,dev.interactive,devAskNewPage,gray,grey,
heat.colors,terrain.colors,topo.colors,axisTicks)
importFrom(graphics,abline,axis,axTicks,box,contour,hist,image,lines,
mtext, par, persp,plot,points,
polygon,rect,strheight,strwidth,text,title)
importFrom(utils, setTxtProgressBar, txtProgressBar)
importFrom(splines,interpSpline)
importFrom(stats,.checkMFClasses,.getXlevels,anova,approx,as.formula,
binomial,coef,contrasts,"contrasts<-",cooks.distance,cor,cov,
......
......@@ -436,7 +436,8 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
## and to control the step length to ensure that at the end of the step we
## are not going uphill w.r.t. the REML criterion...
y <- mf[[gp$response]]
#y <- mf[[gp$response]]
y <- G$y
weights <- G$w
conv <- FALSE
nobs <- nrow(mf)
......@@ -447,7 +448,8 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
pini <- if (is.null(G$family$preinitialize)) NULL else G$family$preinitialize(y,G$family)
if (!is.null(pini$Theta)) G$family$putTheta(pini$Theta)
if (!is.null(pini$y)) y <- pini$y
scale <- if (is.null(G$family$scale)) 1 else G$family$scale
if (is.null(G$family$scale)) scale <- 1 else scale <- if (G$family$scale<0) scale else G$family$scale
scale1 <- scale
if (scale < 0) scale <- var(y) *.1 ## initial guess
} else efam <- FALSE
......@@ -578,9 +580,9 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
if (efam) { ## extended family
if (iter>1) { ## estimate theta
scale1 <- if (!is.null(family$scale)) family$scale else scale
if (family$n.theta>0||scale<0) theta <- estimate.theta(theta,family,y,mu,scale=scale1,wt=G$w,tol=1e-7)
if (!is.null(family$scale) && family$scale<0) {
#scale1 <- if (!is.null(family$scale)) family$scale else scale
if (family$n.theta>0||scale1<0) theta <- estimate.theta(theta,family,y,mu,scale=scale1,wt=G$w,tol=1e-7)
if (!is.null(family$scale) && scale1<0) {
scale <- exp(theta[family$n.theta+1])
theta <- theta[1:family$n.theta]
}
......@@ -791,8 +793,9 @@ regular.Sb <- function(S,off,sp,beta) {
bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etastart = NULL,
mustart = NULL, offset = rep(0, nobs), control = gam.control(), intercept = TRUE,
cl = NULL,gc.level=0,use.chol=FALSE,nobs.extra=0,samfrac=1,npt=1)
{ y <- mf[[gp$response]]
cl = NULL,gc.level=0,use.chol=FALSE,nobs.extra=0,samfrac=1,npt=1) {
#y <- mf[[gp$response]]
y <- G$y
weights <- G$w
conv <- FALSE
nobs <- nrow(mf)
......@@ -805,7 +808,8 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
pini <- if (is.null(G$family$preinitialize)) NULL else G$family$preinitialize(y,G$family)
if (!is.null(pini$Theta)) G$family$putTheta(pini$Theta)
if (!is.null(pini$y)) y <- pini$y
scale <- if (is.null(G$family$scale)) 1 else G$family$scale
if (is.null(G$family$scale)) scale <- 1 else scale <- if (G$family$scale<0) scale else G$family$scale
scale1 <-scale
if (scale < 0) scale <- var(y) *.1 ## initial guess
} else efam <- FALSE
......@@ -1086,9 +1090,9 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
}
if (efam && iter>1) { ## estimate theta
scale1 <- if (!is.null(family$scale)) family$scale else scale
if (family$n.theta>0||scale<0) theta <- estimate.theta(theta,family,G$y,linkinv(eta),scale=scale1,wt=G$w,tol=1e-7)
if (!is.null(family$scale) && family$scale<0) {
#scale1 <- if (!is.null(family$scale)) family$scale else scale
if (family$n.theta>0||scale1<0) theta <- estimate.theta(theta,family,G$y,linkinv(eta),scale=scale1,wt=G$w,tol=1e-7)
if (!is.null(family$scale) && scale1<0) {
scale <- exp(theta[family$n.theta+1])
theta <- theta[1:family$n.theta]
}
......@@ -1377,9 +1381,9 @@ predict.bam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclu
bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
cl=NULL,gc.level=0,use.chol=FALSE,npt=1)
cl=NULL,gc.level=0,use.chol=FALSE,npt=1) {
## function that does big additive model fit in strictly additive case
{ ## first perform the QR decomposition, blockwise....
## first perform the QR decomposition, blockwise....
n <- nrow(mf)
if (rho!=0) { ## AR1 error model
ld <- 1/sqrt(1-rho^2) ## leading diagonal of root inverse correlation
......@@ -1526,7 +1530,7 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
yX.last <- res[[n.threads]]$yX.last
}
G$n <- n
G$y <- mf[[gp$response]]
#G$y <- mf[[gp$response]]
} else { ## n <= chunk.size
if (rho==0) qrx <- qr.update(sqrt(G$w)*G$X,sqrt(G$w)*(G$y-G$offset),use.chol=use.chol,nt=npt) else {
......@@ -1946,9 +1950,6 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
if (inherits(family,"extended.family")) {
family <- fix.family.link(family); efam <- TRUE
} else efam <- FALSE
#else {
#if (inherits(family,"extended.family")) stop("used bam(...,discrete=TRUE) with extended families")
#}
if (method%in%c("fREML")&&!is.null(min.sp)) {
min.sp <- NULL
......@@ -1995,11 +1996,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
if (gc.level>0) gc()
mf <- eval(mf, parent.frame()) # the model frame now contains all the data
# if ("matrix"%in%unlist(lapply(mf,class))) {
# mfattr <- attributes(mf)
# mf <- lapply(mf,drop) # avoid single column matrices
# mfattr -> attributes(mf)
# }
if (nrow(mf)<2) stop("Not enough (non-NA) data to do anything meaningful")
terms <- attr(mf,"terms")
if (gc.level>0) gc()
......@@ -2121,7 +2118,6 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
rind <- k:(k+dt[kb] - 1 - by.present)
dk$nr[rind] <- dk$nr[k+G$smooth[[i]]$rind-1]
G$ks[rind,] <- G$ks[k+G$smooth[[i]]$rind-1,] # either this line or next not both
#G$kd[,rind] <- G$kd[,k+G$smooth[[i]]$rind-1]
}
for (j in 1:nmar) {
G$Xd[[k]] <- G$smooth[[i]]$margin[[j]]$X[1:dk$nr[k],,drop=FALSE]
......@@ -2175,7 +2171,8 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
if (is.null(mf$"(weights)")) G$w<-rep(1,n)
else G$w<-mf$"(weights)"
G$y <- mf[[gp$response]]
G$offset <- model.offset(mf)
if (is.null(G$offset)) G$offset <- rep(0,n)
......@@ -2198,17 +2195,28 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
environment(G$formula) <- .BaseNamespaceEnv
} else { ## G supplied
scale <- G$scale
if (scale<=0) scale <- G$scale
efam <- G$efam
mf <- G$mf; G$mf <- NULL
gp <- G$gp; G$gp <- NULL
na.action <- G$na.action; G$na.action <- NULL
if (!is.null(sp)&&any(sp>=0)) { ## request to modify smoothing parameters
if (is.null(G$L)) G$L <- diag(length(G$sp))
if (length(sp)!=ncol(G$L)) stop('length of sp must be number of free smoothing parameters in original model')
ind <- sp>=0 ## which smoothing parameters are now fixed
spind <- log(sp[ind]);
spind[!is.finite(spind)] <- -30 ## set any zero parameters to effective zero
G$lsp0 <- G$lsp0 + drop(G$L[,ind,drop=FALSE] %*% spind) ## add fix to lsp0
G$L <- G$L[,!ind,drop=FALSE] ## drop the cols of G
G$sp <- rep(-1,ncol(G$L))
}
} ## end of G setup
if (!fit) {
G$efam <- efam
G$scale <- scale
G$mf <- mf;G$na.action <- na.action;G$gp <- gp
class(G) <- "bam.prefit"
return(G)
}
......
......@@ -25,11 +25,29 @@ cox.ph <- function (link = "identity") {
## initialization is tough here... need data frame in reverse time order,
## and intercept removed from X...
preinitialize <- expression({
# preinitialize <- expression({
## code to evaluate in estimate.gam...
## sort y (time) into decending order, and
## re-order weights and rows of X accordingly
## matrix y has strat as second column
# G$family$data <- list()
# y.order <- if (is.matrix(G$y)) order(G$y[,2],G$y[,1],decreasing=TRUE) else
# order(G$y,decreasing=TRUE)
# G$family$data$y.order <- y.order
# G$y <- if (is.matrix(G$y)) G$y[y.order,] else G$y[y.order]
# attrX <- attributes(G$X)
# G$X <- G$X[y.order,,drop=FALSE]
# attributes(G$X) <- attrX
# G$w <- G$w[y.order]
# G$offset <- G$offset[y.order]
# })
preinitialize <- function(G) {
## G is a gam pre-fit object. Pre-initialize can manipulate some of its
## elements, returning a named list of the modified ones.
## sort y (time) into decending order, and
## re-order weights and rows of X accordingly
## matrix y has strat as second column
G$family$data <- list()
y.order <- if (is.matrix(G$y)) order(G$y[,2],G$y[,1],decreasing=TRUE) else
order(G$y,decreasing=TRUE)
......@@ -40,7 +58,8 @@ cox.ph <- function (link = "identity") {
attributes(G$X) <- attrX
G$w <- G$w[y.order]
G$offset <- G$offset[y.order]
})
list(family=G$family,y=G$y,X=G$X,w=G$w,offset=G$offset)
} ## preinitialize
postproc <- expression({
## code to evaluate in estimate.gam, to do with data ordering and
......@@ -59,6 +78,7 @@ cox.ph <- function (link = "identity") {
object$null.deviance <- ## sum of squares of null deviance residuals
2*sum(abs((object$prior.weights + log(s0.base) + object$prior.weights*(log(-log(s0.base))))))
## and undo the re-ordering...
y.order <- G$family$data$y.order
object$linear.predictors[y.order] <- object$linear.predictors
object$fitted.values[y.order] <- object$fitted.values
if (is.matrix(object$y)) object$y[y.order,] <- object$y else object$y[y.order] <- object$y
......
......@@ -45,11 +45,19 @@ estimate.theta <- function(theta,family,y,mu,scale=1,wt=1,tol=1e-7,attachH=FALSE
} else H <- NULL
list(nll=nll,g=g,H=H)
} ## nlogl
if (scale>=0&&family$n.theta==0) stop("erroneous call to estimate.theta - no free parameters")
n.theta <- length(theta) ## dimension of theta vector (family$n.theta==0 => all fixed)
del.ind <- 1:n.theta
if (scale<0) theta <- c(theta,log(var(y)*.1))
nll <- nlogl(theta,family,y,mu,scale,wt,2)
g <- if (family$n.theta==0) nll$g[-del.ind] else nll$g
H <- if (family$n.theta==0) nll$H[-del.ind,-del.ind,drop=FALSE] else nll$H
step.failed <- FALSE
for (i in 1:100) { ## main Newton loop
eh <- eigen(nll$H,symmetric=TRUE)
#H <- if (family$n.theta==0) nll$H[-del.ind,-del.ind,drop=FALSE] else nll$H
#g <- if (family$n.theta==0) nll$g[-del.ind] else nll$g
eh <- eigen(H,symmetric=TRUE)
pdef <- sum(eh$values <= 0)==0
if (!pdef) { ## Make the Hessian pos def
eh$values <- abs(eh$values)
......@@ -57,10 +65,11 @@ estimate.theta <- function(theta,family,y,mu,scale=1,wt=1,tol=1e-7,attachH=FALSE
eh$values[eh$values<thresh] <- thresh
}
## compute step = -solve(H,g) (via eigen decomp)...
step <- - drop(eh$vectors %*% ((t(eh$vectors) %*% nll$g)/eh$values))
step <- - drop(eh$vectors %*% ((t(eh$vectors) %*% g)/eh$values))
if (family$n.theta==0) step <- c(rep(0,n.theta),step)
## limit the step length...
ms <- max(abs(step))
if (ms>4) step <- step*4/ms
if (ms>4) step <- step*4/ms
nll1 <- nlogl(theta+step,family,y,mu,scale,wt,2)
iter <- 0
......@@ -76,11 +85,13 @@ estimate.theta <- function(theta,family,y,mu,scale=1,wt=1,tol=1e-7,attachH=FALSE
theta <- theta + step ## accept updated theta
## accept log lik and derivs ...
nll <- if (iter>0) nlogl(theta,family,y,mu,scale,wt,2) else nll1
g <- if (family$n.theta==0) nll$g[-del.ind] else nll$g
H <- if (family$n.theta==0) nll$H[-del.ind,-del.ind,drop=FALSE] else nll$H
## convergence checking...
if (sum(abs(nll$g) > tol*abs(nll$nll))==0) break
if (sum(abs(g) > tol*abs(nll$nll))==0) break
} ## main Newton loop
if (step.failed) warning("step failure in theta estimation")
if (attachH) attr(theta,"H") <- nll$H
if (attachH) attr(theta,"H") <- H#nll$H
theta
} ## estimate.theta
......
......@@ -123,6 +123,7 @@ Sl.setup <- function(G) {
Sl[[b]]$start <- G$smooth[[i]]$first.para + sbStart[j]-1
Sl[[b]]$stop <- G$smooth[[i]]$first.para + sbStop[j]-1
Sl[[b]]$rank <- G$smooth[[i]]$rank[j]
Sl[[b]]$lambda <- 1 ## dummy here
Sl[[b]]$repara <- TRUE ## signals ok to linearly reparameterize
if (!is.null(G$smooth[[i]]$g.index)) { ## then some parameters are non-linear - can't re-param
if (any(G$smooth[[i]]$g.index[ind])) Sl[[b]]$repara <- FALSE
......@@ -774,7 +775,6 @@ Sl.ift <- function(Sl,R,X,y,beta,piv,rp) {
}
XX.db <- t(X)%*%(X%*%db)
S.db <- Sl.mult(Sl,db,k=0)
## Sk.db <- Sl.termMult(Sl,db,full=TRUE) ## Sk.db[[j]][,k] is S_j d beta / d rho_k
rss2 <- bSb2 <- matrix(0,nd,nd)
for (k in 1:nd) { ## second derivative loop
......@@ -805,39 +805,16 @@ Sl.iftChol <- function(Sl,XX,R,d,beta,piv,nt=1) {
## timing close to identical - modified for parallel exec
D <- matrix(unlist(Skb),nrow(XX),nd)
bSb1 <- colSums(beta*D)
#D <- D[piv,]/d[piv]
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]
#db[piv,] <- -backsolve(R,forwardsolve(t(R),D))/d[piv]
## original serial - a bit slow with very large numbers of smoothing
## parameters....
#Rt <- t(R) ## essential to do this first, or t(R) dominates cost in loop
#for (i in 1:nd) { ## compute the first derivatives
# db[piv,i] <- -backsolve(R,forwardsolve(Rt,Skb[[i]][piv]/d[piv]))/d[piv] ## d beta/ d rho_i
# bSb1[i] <- sum(beta*Skb[[i]]) ## d b'Sb / d_rho_i
#}
## XX.db <- XX%*%db
XX.db <- .Call(C_mgcv_pmmult2,XX,db,0,0,nt)
#XX.db[piv,] <- d[piv]*(t(R)%*%(R%*%(d[piv]*db[piv,]))) ## X'Xdb
S.db <- Sl.mult(Sl,db,k=0)
##Sk.db <- Sl.termMult(Sl,db,full=TRUE) ## Sk.db[[j]][,k] is S_j d beta / d rho_k
## following loop very slow if large numbers of smoothing parameters...
#rss2 <- bSb2 <- matrix(0,nd,nd)
#for (k in 1:nd) { ## second derivative loop
# for (j in k:nd) {
# rss2[j,k] <- rss2[k,j] <- 2 * sum(db[,j]*XX.db[,k])
# bSb2[j,k] <- bSb2[k,j] <- (k==j)*sum(beta*Skb[[k]]) + 2*(sum(db[,k]*(Skb[[j]]+S.db[,j])) +
# sum(db[,j]*Skb[[k]]))
# }
#}
## rss2 <- 2 * t(db) %*% XX.db
rss2 <- 2 * .Call(C_mgcv_pmmult2,db,XX.db,1,0,nt)
bSb2 <- diag(x=colSums(beta*D),nrow=nd)
## bSb2 <- bSb2 + 2*(t(db)%*%(D+S.db) + t(D)%*%db)
bSb2 <- bSb2 + 2 * (.Call(C_mgcv_pmmult2,db,D+S.db,1,0,nt) + .Call(C_mgcv_pmmult2,D,db,1,0,nt))
list(bSb=sum(beta*Sb),bSb1=bSb1,bSb2=bSb2,
d1b=db ,rss1=rss1,rss2=rss2)
......
......@@ -1262,6 +1262,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
## initial fit
initial.lsp <- lsp ## used if edge correcting to set direction of correction
b<-gam.fit3(x=X, y=y, sp=L%*%lsp+lsp0,Eb=Eb,UrS=UrS,
offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=2,
control=control,gamma=gamma,scale=scale,printWarn=FALSE,start=start,
......@@ -1627,18 +1628,20 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
REML <- b$REML
alpha <- if (is.logical(edge.correct)) .02 else abs(edge.correct) ## target RE/ML change per sp
b1 <- b; lsp1 <- lsp
if (length(flat)) for (i in flat) {
REML <- b1$REML + alpha
while (b1$REML < REML) {
lsp1[i] <- lsp1[i] - 1
b1 <- gam.fit3(x=X, y=y, sp=L%*%lsp1+lsp0,Eb=Eb,UrS=UrS,
if (length(flat)) {
step <- as.numeric(initial.lsp - lsp)*2-1
for (i in flat) {
REML <- b1$REML + alpha
while (b1$REML < REML) {
lsp1[i] <- lsp1[i] + step[i]
b1 <- 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=0,
control=control,gamma=gamma,scale=scale,printWarn=FALSE,start=start,
mustart=mustart,scoreType=scoreType,null.coef=null.coef,
pearson.extra=pearson.extra,dev.extra=dev.extra,n.true=n.true,Sl=Sl,...)
}
}
}
} ## if length(flat)
b1 <- 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,
control=control,gamma=gamma,scale=scale,printWarn=FALSE,start=start,
......@@ -2837,7 +2840,7 @@ ldTweedie <- function(y,mu=y,p=1.5,phi=1,rho=NA,theta=NA,a=1.001,b=1.999,all.der
## evaluates log Tweedie density for 1<=p<=2, using series summation of
## Dunn & Smyth (2005) Statistics and Computing 15:267-280.
n <- length(y)
if (!is.na(rho)&&!is.na(theta)) { ## use rho and theta and get derivs w.r.t. these
if (all(!is.na(rho))&&all(!is.na(theta))) { ## use rho and theta and get derivs w.r.t. these
#if (length(rho)>1||length(theta)>1) stop("only scalar `rho' and `theta' allowed.")
if (a>=b||a<=1||b>=2) stop("1<a<b<2 (strict) required")
work.param <- TRUE
......
......@@ -689,7 +689,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
names(residuals) <- ynames
wtdmu <- sum(weights * y)/sum(weights) ## has to then be corrected when this is incorrect
## wtdmu <- sum(weights * mu)/sum(weights) ## changed from y
nulldev <- sum(dev.resids(y, rep(wtdmu,length(y)), weights))
nulldev <- sum(dev.resids(y, rep(wtdmu,length(y)), weights)) ## this will be corrected in family postproc
n.ok <- nobs - sum(weights == 0)
nulldf <- n.ok
ww <- wt <- rep.int(0, nobs)
......@@ -921,7 +921,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
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) {
if (check.deriv) { ## code for checking derivatives when debugging
fdg <- ll$lb*0; fdh <- ll$lbb*0
for (k in 1:length(coef)) {
coef1 <- coef;coef1[k] <- coef[k] + eps
......@@ -929,14 +929,22 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
fdg[k] <- (ll.fd$l-ll$l)/eps
fdh[,k] <- (ll.fd$lb-ll$lb)/eps
}
}
} ## derivative checking end
grad <- ll$lb - St%*%coef
Hp <- -ll$lbb+St
D <- diag(Hp)
if (sum(!is.finite(D))>0) stop("non finite values in Hessian")
indefinite <- FALSE
if (sum(D <= 0)) { ## Hessian indefinite, for sure
D <- rep(1,ncol(Hp))
if (min(D)<0) { ## 2/2/19 replaces any D<0 indicating indef
Dthresh <- max(D)*sqrt(.Machine$double.eps)
if (-min(D) < Dthresh) { ## could be indef or +ve semi def
indefinite <- FALSE
D[D<Dthresh] <- Dthresh
} else indefinite <- TRUE ## min D too negative to be semi def
} else indefinite <- FALSE ## On basis of D could be +ve def
if (indefinite) { ## Hessian indefinite, for sure
## D <- rep(1,ncol(Hp)) # moved to later otherwise Ip/Ib pointless 2/2/19
if (eigen.fix) {
eh <- eigen(Hp,symmetric=TRUE);
ev <- abs(eh$values)
......@@ -946,6 +954,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
Ip <- diag(rank)*abs(max(D)*.Machine$double.eps^.5)
Hp <- Hp + Ip + Ib
}
D <- rep(1,ncol(Hp)) ## 2/2/19
indefinite <- TRUE
} else { ## Hessian could be +ve def in which case Choleski is cheap!
D <- D^-.5 ## diagonal pre-conditioner
......@@ -1080,8 +1089,9 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
}
} ## end of main fitting iteration
## at this stage the Hessian (of pen lik. w.r.t. coefs) should be +ve definite,
## at this stage the Hessian (of pen lik. w.r.t. coefs) should be +ve semi definite,
## so that the pivoted Choleski factor should exist...
if (iter == 2*control$maxit&&converged==FALSE)
warning(gettextf("iteration limit reached: max abs grad = %g",max(abs(grad))))
......@@ -1163,11 +1173,6 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
k <- k + 1
d2ldetH[i,j] <- -sum(d1Hp[[i]]*t(d1Hp[[j]])) - llr$trHid2H[k]
if (i==j) { ## need to add term relating to smoothing penalty
#A <- t(Sl.mult(rp$Sl,diag(q),i,full=FALSE))
#bind <- rowSums(abs(A))!=0 ## FIX: abs 3/3/16
#ind <- which(bind)
#bind <- bind[!bdrop]
#A <- A[!bdrop,!bdrop[ind]]
A <- Sl.mult(rp$Sl,diag(q),i,full=TRUE)[!bdrop,!bdrop]
bind <- rowSums(abs(A))!=0 ## row/cols of non-zero block
A <- A[,bind] ## drop the zero columns
......@@ -1180,24 +1185,20 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
## Compute derivs of b'Sb...
if (deriv>0) {
# Sb <- St%*%coef
Skb <- Sl.termMult(rp$Sl,fcoef,full=TRUE)
d1bSb <- rep(0,m)