...
 
Commits (3)
......@@ -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)
for (i in 1:m) {
Skb[[i]] <- Skb[[i]][!bdrop]
d1bSb[i] <- # 2*sum(d1b[,i]*Sb) + # cancels
sum(coef*Skb[[i]])
d1bSb[i] <- sum(coef*Skb[[i]])
}
}
if (deriv>1) {
d2bSb <- matrix(0,m,m)
# k <- 0
for (i in 1:m) {
Sd1b <- St%*%d1b[,i]
for (j in i:m) {
k <- k + 1
d2bSb[j,i] <- d2bSb[i,j] <- 2*sum( # d2b[,k]*Sb + # cancels
d2bSb[j,i] <- d2bSb[i,j] <- 2*sum(
d1b[,i]*Skb[[j]] + d1b[,j]*Skb[[i]] + d1b[,j]*Sd1b)
}
d2bSb[i,i] <- d2bSb[i,i] + sum(coef*Skb[[i]])
......@@ -1217,12 +1218,13 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
if (!is.null(REML1)) cat("REML1 =",REML1,"\n")
}
REML2 <- if (deriv<2) NULL else -( (d2l - d2bSb/2)/gamma + rp$ldet2/2 - d2ldetH/2 )
## bSb <- t(coef)%*%St%*%coef
## Get possibly multiple linear predictors
lpi <- attr(x,"lpi")
if (is.null(lpi)) {
if (is.null(lpi)) { ## only one...
linear.predictors <- if (is.null(offset)) as.numeric(x%*%coef) else as.numeric(x%*%coef+offset)
fitted.values <- family$linkinv(linear.predictors)
} else {
} else { ## multiple...
fitted.values <- linear.predictors <- matrix(0,nrow(x),length(lpi))
if (!is.null(offset)) offset[[length(lpi)+1]] <- 0
for (j in 1:length(lpi)) {
......@@ -1280,6 +1282,7 @@ efsud <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,family,
fit <- gam.fit5(x=x,y=y,lsp=lsp,Sl=Sl,weights=weights,offset=offset,deriv=0,family=family,
control=control,Mp=Mp,start=start,gamma=1)
score.hist <- rep(0,200)
tiny <- .Machine$double.eps^.5 ## used to bound above zero
for (iter in 1:200) {
start <- fit$coefficients
## obtain Vb...
......@@ -1305,8 +1308,9 @@ efsud <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,family,
for (i in 1:length(Sb)) {
bSb[i] <- sum(start*Sb[[i]])
}
a <- pmax(0,fit$S1*exp(-lsp) - trVS)
r <- a/pmax(0,bSb)
a <- pmax(tiny,fit$S1*exp(-lsp) - trVS)
r <- a/pmax(tiny,bSb)
r[a==0&bSb==0] <- 1
r[!is.finite(r)] <- 1e6
lsp1 <- pmin(lsp + log(r)*mult,control$efs.lspmax)
......
This diff is collapsed.
This diff is collapsed.
......@@ -369,8 +369,13 @@ interpret.gam <- function(gf,extra.special=NULL) {
}
av <- unique(av) ## strip out duplicate variable names
pav <- unique(pav)
ret$fake.formula <- if (length(av)>0) reformulate(av,response=ret[[1]]$response) else
ret[[1]]$fake.formula ## create fake formula containing all variables
if (length(av)>0) {
## work around - reformulate with response = "log(x)" will treat log(x) as a name,
## not the call it should be...
fff <- formula(paste(ret[[1]]$response,"~ ."))
ret$fake.formula <- reformulate(av,response=ret[[1]]$response)
ret$fake.formula[[2]] <- fff[[2]] ## fix messed up response
} else ret$fake.formula <- ret[[1]]$fake.formula ## create fake formula containing all variables
ret$pred.formula <- if (length(pav)>0) reformulate(pav) else ~1 ## predictor only formula
ret$response <- ret[[1]]$response
ret$nlp <- nlp ## number of linear predictors
......@@ -1596,7 +1601,7 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,start=NU
method <- "REML" ## any method you like as long as it's REML
G$Sl <- Sl.setup(G) ## prepare penalty sequence
## Xorig <- G$X ## store original X incase it is needed by family - poor option pre=proc can manipulate G$X
G$X <- Sl.initial.repara(G$Sl,G$X,both.sides=FALSE) ## re-parameterize accordingly
if (!is.null(start)) start <- Sl.initial.repara(G$Sl,start,inverse=FALSE,both.sides=FALSE)
......@@ -1690,7 +1695,10 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,start=NU
## extended family may need to manipulate G...
if (!is.null(G$family$preinitialize)) {
if (inherits(G$family,"general.family")) eval(G$family$preinitialize) else {
if (inherits(G$family,"general.family")) {
Gmod <- G$family$preinitialize(G) ## modifies some elements of G
for (gnam in names(Gmod)) G[[gnam]] <- Gmod[[gnam]] ## copy these into G
} else {
## extended family - just initializes theta and possibly y
pini <- G$family$preinitialize(G$y,G$family)
if (!is.null(pini$Theta)) G$family$putTheta(pini$Theta)
......@@ -2014,9 +2022,23 @@ gam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
G$formula <- formula
G$pred.formula <- gp$pred.formula
environment(G$formula)<-environment(formula)
} else { ## G not 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))
}
}
if (!fit) return(G)
if (!fit) {
class(G) <- "gam.prefit"
return(G)
}
if (ncol(G$X)>nrow(G$X)) stop("Model has more coefficients than data")
......@@ -2932,7 +2954,7 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclu
#lpi <- list();pst <- c(pstart,ncol(X)+1)
#for (i in 1:(length(pst)-1)) lpi[[i]] <- pst[i]:(pst[i+1]-1)
attr(X,"lpi") <- lpi
ffv <- fam$predict(fam,se.fit,y=response,X=X,beta=object$coefficients,
ffv <- fam$predict(fam,se.fit,y=response[start:stop],X=X,beta=object$coefficients,
off=offs,Vb=object$Vp)
if (is.matrix(fit)&&!is.matrix(ffv[[1]])) {
fit <- fit[,1]; if (se.fit) se <- se[,1]
......@@ -2959,7 +2981,7 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclu
if (se.fit) se[start:stop]<-se[start:stop]*abs(dmu.deta(fit[start:stop]))
fit[start:stop] <- linkinv(fit[start:stop])
} else { ## family has its own prediction code for response case
ffv <- fam$predict(fam,se.fit,y=response,X=X,beta=object$coefficients,off=offs,Vb=object$Vp)
ffv <- fam$predict(fam,se.fit,y=response[start:stop],X=X,beta=object$coefficients,off=offs,Vb=object$Vp)
if (is.null(fit1)&&is.matrix(ffv[[1]])) {
fit1 <- matrix(0,np,ncol(ffv[[1]]))
if (se.fit) se1 <- fit1
......@@ -3840,11 +3862,17 @@ cooks.distance.gam <- function(model,...)
}
sp.vcov <- function(x) {
sp.vcov <- function(x,edge.correct=TRUE,reg=1e-3) {
## get cov matrix of smoothing parameters, if available
if (!inherits(x,"gam")) stop("argument is not a gam object")
if (x$method%in%c("ML","P-ML","REML","P-REML","fREML")&&!is.null(x$outer.info$hess)) {
return(solve(x$outer.info$hess))
hess <- x$outer.info$hess
p <- ncol(hess)
if (edge.correct&&!is.null(attr(hess,"hess1"))) {
V <- solve(attr(hess,"hess1")+diag(p)*reg)
attr(V,"lsp") <- attr(hess,"lsp1")
return(V)
} else return(solve(hess+reg))
} else return(NULL)
}
......
......@@ -242,6 +242,19 @@ choldrop <- function(R,k) {
Rup
} ## choldrop
cholup <- function(R,u,up=TRUE) {
## routine to update Cholesky factor R to the factor of R'R + uu' (up == TRUE)
## or R'R - uu' (up=FALSE).
n <- as.integer(ncol(R))
up <- as.integer(up)
eps <- as.double(.Machine$double.eps)
R1 <- R * 1.0
.Call(C_mgcv_chol_up,R1,u,n,up,eps)
if (up==0) if ((n>1 && R1[2,1] < -1)||(n==1&&u[1]>R[1])) stop("update not positive definite")
R1
} ## cholup
vcorr <- function(dR,Vr,trans=TRUE) {
## Suppose b = sum_k op(dR[[k]])%*%z*r_k, z ~ N(0,Ip), r ~ N(0,Vr). vcorr returns cov(b).
## dR is a list of p by p matrices. 'op' is 't' if trans=TRUE and I() otherwise.
......
......@@ -88,10 +88,44 @@ mvn <- function(d=2) {
## initialization has to add in the extra parameters of
## the cov matrix...
preinitialize <- expression({
# preinitialize <- expression({
## code to evaluate in estimate.gam...
## extends model matrix with dummy columns and
## finds initial coefficients
# ydim <- ncol(G$y) ## dimension of response
# nbeta <- ncol(G$X)
# ntheta <- ydim*(ydim+1)/2 ## number of cov matrix factor params
# lpi <- attr(G$X,"lpi")
# #offs <- attr(G$X,"offset")
# XX <- crossprod(G$X)
# G$X <- cbind(G$X,matrix(0,nrow(G$X),ntheta)) ## add dummy columns to G$X
# #G$cmX <- c(G$cmX,rep(0,ntheta)) ## and corresponding column means
# G$term.names <- c(G$term.names,paste("R",1:ntheta,sep="."))
# attr(G$X,"lpi") <- lpi
# #offs -> attr(G$X,"offset")
# attr(G$X,"XX") <- XX
# ## pad out sqrt of balanced penalty matrix to account for extra params
# attr(G$Sl,"E") <- cbind(attr(G$Sl,"E"),matrix(0,nbeta,ntheta))
# G$family$data <- list(ydim = ydim,nbeta=nbeta)
# G$family$ibeta = rep(0,ncol(G$X))
# ## now get initial parameters and store in family...
# for (k in 1:ydim) {
# sin <- G$off %in% lpi[[k]]
# #Sk <- G$S[sin]
# um <- magic(G$y[,k],G$X[,lpi[[k]]],rep(-1,sum(sin)),G$S[sin],
# match(G$off[sin],lpi[[k]]), ## turn G$off global indices into indices for this predictor
# nt=control$nthreads)
# G$family$ibeta[lpi[[k]]] <- um$b
# G$family$ibeta[nbeta+1] <- -.5*log(um$scale) ## initial log root precision
# nbeta <- nbeta + ydim - k + 1
# }
# })
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.
## extends model matrix with dummy columns and
## finds initial coefficients
ydim <- ncol(G$y) ## dimension of response
nbeta <- ncol(G$X)
ntheta <- ydim*(ydim+1)/2 ## number of cov matrix factor params
......@@ -105,7 +139,7 @@ mvn <- function(d=2) {
#offs -> attr(G$X,"offset")
attr(G$X,"XX") <- XX
## pad out sqrt of balanced penalty matrix to account for extra params
attr(G$Sl,"E") <- cbind(attr(G$Sl,"E"),matrix(0,nbeta,ntheta))
if (!is.null(G$Sl)) attr(G$Sl,"E") <- cbind(attr(G$Sl,"E"),matrix(0,nbeta,ntheta))
G$family$data <- list(ydim = ydim,nbeta=nbeta)
G$family$ibeta = rep(0,ncol(G$X))
## now get initial parameters and store in family...
......@@ -113,15 +147,17 @@ mvn <- function(d=2) {
sin <- G$off %in% lpi[[k]]
#Sk <- G$S[sin]
um <- magic(G$y[,k],G$X[,lpi[[k]]],rep(-1,sum(sin)),G$S[sin],
match(G$off[sin],lpi[[k]]), ## turn G$off global indices into indices for this predictor
nt=control$nthreads)
match(G$off[sin],lpi[[k]])) # , ## turn G$off global indices into indices for this predictor
#nt=control$nthreads)
G$family$ibeta[lpi[[k]]] <- um$b
G$family$ibeta[nbeta+1] <- -.5*log(um$scale) ## initial log root precision
nbeta <- nbeta + ydim - k + 1
}
})
postproc <- expression({
list(X=G$X,term.names=G$term.names,family=G$family)
} ## preinitialize
postproc <- expression({
## code to evaluate in estimate.gam, to do with estimated factor of
## precision matrix, etc...
ydim <- G$family$data$ydim
......
## R routines for the package mgcv (c) Simon Wood 2000-2016
## R routines for the package mgcv (c) Simon Wood 2000-2019
## This file is primarily concerned with defining classes of smoother,
## via constructor methods and prediction matrix methods. There are
......
mgcv (1.8-27-1) unstable; urgency=medium
* New upstream release
* debian/control: Set Build-Depends: to current R version
* debian/control: Set Standards-Version: to current version
-- Dirk Eddelbuettel <edd@debian.org> Sat, 09 Feb 2019 08:47:04 -0600
mgcv (1.8-26-1) unstable; urgency=medium
* New upstream release
......
......@@ -2,8 +2,8 @@ Source: mgcv
Section: gnu-r
Priority: optional
Maintainer: Dirk Eddelbuettel <edd@debian.org>
Build-Depends: debhelper (>= 7.0.0), r-base-dev (>= 3.5.1), r-cran-nlme, r-cran-matrix, dh-r
Standards-Version: 4.2.1
Build-Depends: debhelper (>= 7.0.0), r-base-dev (>= 3.5.2), r-cran-nlme, r-cran-matrix, dh-r
Standards-Version: 4.3.0
Vcs-Browser: https://salsa.debian.org/edd/r-cran-mgcv
Vcs-Git: https://salsa.debian.org/edd/r-cran-mgcv.git
Homepage: https://cran.r-project.org/package=mgcv
......
\name{FFdes}
\alias{FFdes}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Level 5 fractional factorial designs}
\description{Computes level 5 fractional factorial designs for up to 120 factors
using the agorithm of Sanchez and Sanchez (2005), and optionally central composite designs.
}
\usage{
FFdes(size=5,ccd=FALSE)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{size}{number of factors up to 120.}
\item{ccd}{if \code{TRUE}, adds points along each axis at the same distance from the origin as the points in the
fractional factorial design, to create the outer points of a central composite design. Add central points to complete.}
}
\details{Basically a translation of the code provided in the appendix of Sanchez and Sanchez (2005).
}
\references{
Sanchez, S. M. & Sanchez, P. J. (2005) Very large fractional factorial and central composite designs.
ACM Transactions on Modeling and Computer Simulation. 15: 362-377
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}}
\examples{
require(mgcv)
plot(rbind(0,FFdes(2,TRUE)),xlab="x",ylab="y",
col=c(2,1,1,1,1,4,4,4,4),pch=19,main="CCD")
FFdes(5)
FFdes(5,TRUE)
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..
......@@ -151,7 +151,7 @@ be reduced by first fitting a model to a random sample of the data, and using th
involving factor variables you might want to turn this off. Only do so if you know what you are doing.}
\item{G}{if not \code{NULL} then this should be the object returned by a previous call to \code{bam} with
\code{fit=FALSE}. Causes all other arguments to be ignored except \code{chunk.size}, \code{gamma},\code{nthreads}, \code{cluster}, \code{rho}, \code{gc.level}, \code{samfrac}, \code{use.chol} and \code{method}.}
\code{fit=FALSE}. Causes all other arguments to be ignored except \code{sp}, \code{chunk.size}, \code{gamma},\code{nthreads}, \code{cluster}, \code{rho}, \code{gc.level}, \code{samfrac}, \code{use.chol}, \code{method} and \code{scale} (if >0).}
\item{fit}{if \code{FALSE} then the model is set up for fitting but not estimated, and an object is returned, suitable for passing as the \code{G} argument to \code{bam}.}
......
\name{choldrop}
\alias{choldrop}
\alias{cholup}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Down date a Cholesky factor on dropping a row/col}
\description{Given a Cholesky factor, \code{R}, of a matrix, \code{A}, finds the Cholesky factor of \code{A[-k,-k]},
where \code{k} is an integer.
\title{Deletion and rank one Cholesky factor update}
\description{Given a Cholesky factor, \code{R}, of a matrix, \code{A}, \code{choldrop} finds the Cholesky factor of \code{A[-k,-k]},
where \code{k} is an integer. \code{cholup} finds the factor of \eqn{A + uu^T}{A+uu'} (update) or \eqn{A - uu^T}{A-uu'} (downdate).
}
\usage{
choldrop(R,k)
cholup(R,u,up)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{R}{Cholesky factor of a matrix, \code{A}.}
\item{k}{row and column of \code{A} to drop.}
\item{u}{vector defining rank one update.}
\item{up}{if \code{TRUE} compute update, otherwise downdate.}
}
\details{If \code{R} is upper triangular then \code{t(R[,-k])\%*\%R[,-k] == A[-k,-k]}, but \code{R[,-k]} has elements on the first sub-diagonal, from its kth column onwards. To get from this to a triangular Cholesky factor of \code{A[-k,-k]} we can apply a sequence of Given rotations from the left to eliminate the sub-diagonal elements. The routine does this. If \code{R} is a lower triangular factor then Givens rotations from the right are needed to remove the extra elements. If \code{n} is the dimension of \code{R} then the update has O(n^2) computational cost.
\details{First consider \code{choldrop}. If \code{R} is upper triangular then \code{t(R[,-k])\%*\%R[,-k] == A[-k,-k]}, but \code{R[,-k]} has elements on the first sub-diagonal, from its kth column onwards. To get from this to a triangular Cholesky factor of \code{A[-k,-k]} we can apply a sequence of Givens rotations from the left to eliminate the sub-diagonal elements. The routine does this. If \code{R} is a lower triangular factor then Givens rotations from the right are needed to remove the extra elements. If \code{n} is the dimension of \code{R} then the update has \eqn{O(n^2)}{O(n^2)} computational cost.
Note that the update is vector oriented, and is hence not susceptible to speed up by use of an optimized BLAS. The update is set up to be relatively Cache friendly, in that in the upper triangular case successive Givens rotations are stored for sequential application columnwise, rather than being applied rowwise as soon as they are computed. Even so, the upper triangular update is slightly slower than the lower triangular update.
\code{cholup} (which assumes \code{R} is upper triangular) updates based on the observation that \eqn{ R^TR + uu^T = [u,R^T][u,R^T]^T = [u,R^T]Q^TQ[u,R^T]^T}{R'R + uu' = [u,R'][u,R']' = [u,R']Q'Q[u,R']'}, and therefore we can construct \eqn{Q}{Q} so that \eqn{Q[u,R^T]^T=[0,R_1^T]^T}{Q[u,R']'=[0,R1']'}, where \eqn{R_1}{R1} is the modified factor. \eqn{Q}{Q} is constructed from a sequence of Givens rotations in order to zero the elements of \eqn{u}{u}. Downdating is similar except that hyperbolic rotations have to be used in place of Givens rotations --- see Golub and van Loan (2013, section 6.5.4) for details. Downdating only works if \eqn{A - uu^T}{A-uu'} is positive definite. Again the computational cost is \eqn{O(n^2)}{O(n^2)}.
Note that the updates are vector oriented, and are hence not susceptible to speed up by use of an optimized BLAS. The updates are set up to be relatively Cache friendly, in that in the upper triangular case successive Givens rotations are stored for sequential application column-wise, rather than being applied row-wise as soon as they are computed. Even so, the upper triangular update is slightly slower than the lower triangular update.
}
\references{
Golub GH and CF Van Loan (2013) Matrix Computations (4th edition) Johns Hopkins
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}}
......@@ -39,6 +47,19 @@ Note that the update is vector oriented, and is hence not susceptible to speed u
Ld <- choldrop(L,k)
range(Ld-t(chol(A[-k,-k])))
Ld;t(chol(A[-k,-k]))
## Rank one update example
u <- runif(n)
R <- cholup(R0,u,TRUE)
Ru <- chol(A+u \%*\% t(u)) ## direct for comparison
R;Ru
range(R-Ru)
## Downdate - just going back from R to R0
Rd <- cholup(R,u,FALSE)
R0;Rd
range(R-Ru)
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..
......
......@@ -32,8 +32,8 @@ censoring.
Stratification is possible, allowing for different baseline hazards in different strata. In that case the response has two columns: the first is event/censoring time and the second is a numeric stratum index. See below for an example.
Prediction from the fitted model object (using the \code{predict} method) with \code{type="response"} will predict on the
survivor function scale. See example code below for extracting the cumulative baseline hazard/survival directly.
Martingale or deviance
survivor function scale. This requires evaluation times to be provided as well as covariates (see example). Also see example code
below for extracting the cumulative baseline hazard/survival directly. Martingale or deviance
residuals can be extracted. The \code{fitted.values} stored in the model object are survival function estimates for each
subject at their event/censoring time.
......
......@@ -51,13 +51,13 @@ tdpois <- function(dat,event="z",et="futime",t="day",status="status1",
## Now do the interpolation of covariates to event times...
um <- data.frame(lapply(X=di,FUN=app,t=di[[t]],to=tr))
## Mark the actual event...
if (um[[et]][1]==max(tr)&&um[[status]]==1) um[[event]][nrow(um)] <- 1
if (um[[et]][1]==max(tr)&&um[[status]][1]==1) um[[event]][nrow(um)] <- 1
um[[et]] <- tr ## reset time to relevant event times
dap[start:(start-1+nrow(um)),] <- um ## copy to dap
start <- start + nrow(um)
if (inter) setTxtProgressBar(prg, i)
}
close(prg)
if (inter) close(prg)
dap[1:(start-1),]
} ## tdpois
......
......@@ -137,7 +137,7 @@ would be required to fit is returned is returned. See argument \code{G}.}
\item{G}{Usually \code{NULL}, but may contain the object returned by a previous call to \code{gam} with
\code{fit=FALSE}, in which case all other arguments are ignored except for
\code{gamma}, \code{in.out}, \code{scale}, \code{control}, \code{method} \code{optimizer} and \code{fit}.}
\code{sp}, \code{gamma}, \code{in.out}, \code{scale}, \code{control}, \code{method} \code{optimizer} and \code{fit}.}
\item{in.out}{optional list for initializing outer iteration. If supplied
then this must contain two elements: \code{sp} should be an array of
......
\name{ginla}
\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.
}
\usage{
ginla(G,A=NULL,nk=16,nb=100,J=1,interactive=FALSE,int=0,approx=0)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{G}{A pre-fit gam object, as produced by \code{gam(...,fit=FALSE)} or \code{bam(...,discrete=TRUE,fit=FALSE)}.}
\item{A}{Either a matrix of transforms of the coefficients that are of interest, or an array of indices of the parameters of interest. If \code{NULL} then distributions are produced for all coefficients.}
\item{nk}{Number of values of each coefficient at which to evaluate its log marginal posterior density. These points are then spline interpolated.}
\item{nb}{Number of points at which to evaluate posterior density of coefficients for returning as a gridded function.}
\item{J}{How many determinant updating steps to take in the log determinant approximation step. Not recommended to increase this. }
\item{interactive}{If this is \code{>0} or \code{TRUE} then every approximate posterior is plotted in red, overlaid on the simple Gaussian approximate posterior. If \code{2} then waits for user to press return between each plot. Useful for judging whether anything is gained by using INLA approach. }
\item{int}{0 to skip integration and just use the posterior modal smoothing parameter. >0 for integration using the CCD approach proposed in Rue et al. (2009).}
\item{approx}{0 for full approximation; 1 to update Hessian, but use approximate modes; 2 as 1 and assume constant Hessian. See details.}
}
\value{
A list with elements \code{beta} and \code{density}, both of which are matrices. Each row relates to one coefficient (or linear coefficient combination) of interest. Both matrices have \code{nb} columns. If \code{int!=0} then a further element \code{reml} gives the integration weights used in the CCD integration, with the central point weight given first.
}
\details{Let \eqn{\beta}{b}, \eqn{\theta}{h} and \eqn{y}{y} denote the model coefficients, hyperparameters/smoothing parameters and response data, respectively. In principle, INLA employs Laplace approximations for \eqn{\pi(\beta_i|\theta,y)}{p(b_i|h,y)} and \eqn{\pi(\theta|y)}{p(h|y)} and then obtains the marginal posterior distribution \eqn{\pi(\beta_i|y)}{p(b_i|y)} by intergrating the approximations to \eqn{\pi(\beta_i|\theta,y)\pi(\theta|y)}{p(b_i|h,y)p(h|y)} w.r.t \eqn{\theta}{h} (marginals for the hyperparameters can also be produced). In practice the Laplace approximation for \eqn{\pi(\beta_i|\theta,y)}{p(b_i|h,y)} is too expensive to compute for each \eqn{\beta_i}{b_i} and must itself be approximated. To this end, there are two quantities that have to be computed: the posterior mode \eqn{\beta^*|\beta_i}{b*|b_i} and the determinant of the Hessian of the joint log density \eqn{\log \pi(\beta,\theta,y)}{log p(b,h,y)} w.r.t. \eqn{\beta}{b} at the mode. Rue et al. (2009) originally approximated the posterior conditional mode by the conditional mode implied by a simple Gaussian approximation to the posterior \eqn{\pi(\beta|y)}{p(b|y)}. They then approximated the log determinant of the Hessian as a function of \eqn{\beta_i}{b_i} using a first order Taylor expansion, which is cheap to compute for the sparse model representaiton that they use, but not when using the dense low rank basis expansions used by \code{\link{gam}}. They also offer a more expensive alternative approximation based on computing the log determiannt with respect only to those elements of \eqn{\beta}{b} with sufficiently high correlation with \eqn{\beta_i}{b_i} according to the simple Gaussian posterior approximation: efficiency again seems to rest on sparsity. Wood (2018) suggests computing the required posterior modes exactly, and basing the log determinant approximation on a BFGS update of the Hessian at the unconditional model. The latter is efficient with or without sparsity, whereas the former is a `for free' improvement. Both steps are efficient because it is cheap to obtain the Cholesky factor of \eqn{H[-i,-i]}{H[-i,-i]} from that of \eqn{H}{H} - see \code{\link{choldrop}}. This is the approach taken by this routine.
The \code{approx} argument allows two further approximations to speed up computations. For \code{approx==1} the exact posterior conditional modes are not used, but instead the conditional modes implied by the simple Gaussian posterior approximation. For \code{approx==2} the same approximation is used for the modes and the Hessian is assumed constant. The latter is quite fast as no log joint density gradient evaluations are required.
Note that for many models the INLA estimates are very close to the usual Gaussian approximation to the posterior, the \code{interactive} argument is useful for investigating this issue.
\code{\link{bam}} models are only supported with the \code{disrete=TRUE} option. The \code{discrete=FALSE} approach would be too inefficient. AR1 models are not supported (related arguments are simply ignored).
}
\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.
}
\section{WARNINGS}{
This routine is still somewhat experimental, so details are liable to change. Also currently not all steps are optimally efficient.
The routine is written for relatively expert users.
\code{ginla} is not designed to deal with rank deficient models.
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}}
\examples{
require(mgcv); require(MASS)
## example using a scale location model for the motorcycle data. A simple plotting
## routine is produced first...
plot.inla <- function(x,inla,k=1,levels=c(.025,.1,.5,.9,.975),
lcol = c(2,4,4,4,2),lwd = c(1,1,2,1,1),lty=c(1,1,1,1,1),
xlab="x",ylab="y",cex.lab=1.5) {
## a simple effect plotter, when distributions of function values of
## 1D smooths have been computed
require(splines)
p <- length(x)
betaq <- matrix(0,length(levels),p) ## storage for beta quantiles
for (i in 1:p) { ## work through x and betas
j <- i + k - 1
p <- cumsum(inla$density[j,])*(inla$beta[j,2]-inla$beta[j,1])
## getting quantiles of function values...
betaq[,i] <- approx(p,y=inla$beta[j,],levels)$y
}
xg <- seq(min(x),max(x),length=200)
ylim <- range(betaq)
ylim <- 1.1*(ylim-mean(ylim))+mean(ylim)
for (j in 1:length(levels)) { ## plot the quantiles
din <- interpSpline(x,betaq[j,])
if (j==1) {
plot(xg,predict(din,xg)$y,ylim=ylim,type="l",col=lcol[j],
xlab=xlab,ylab=ylab,lwd=lwd[j],cex.lab=1.5,lty=lty[j])
} else lines(xg,predict(din,xg)$y,col=lcol[j],lwd=lwd[j],lty=lty[j])
}
} ## plot.inla
## set up the model with a `gam' call...
G <- gam(list(accel~s(times,k=20,bs="ad"),~s(times)),
data=mcycle,family=gaulss(),fit=FALSE)
b <- gam(G=G,method="REML") ## regular GAM fit for comparison
## Now use ginla to get posteriors of estimated effect values
## at evenly spaced times. Create A matrix for this...
rat <- range(mcycle$times)
pd0 <- data.frame(times=seq(rat[1],rat[2],length=20))
X0 <- predict(b,newdata=pd0,type="lpmatrix")
X0[,21:30] <- 0
pd1 <- data.frame(times=seq(rat[1],rat[2],length=10))
X1 <- predict(b,newdata=pd1,type="lpmatrix")
X1[,1:20] <- 0
A <- rbind(X0,X1) ## A maps coefs to required function values
## call ginla. Set int to 1 for integrated version.
## Set interactive = 1 or 2 to plot marginal posterior distributions
## (red) and simple Gaussian approximation (black).
inla <- ginla(G,A,int=0)
par(mfrow=c(1,2),mar=c(5,5,1,1))
fv <- predict(b,se=TRUE) ## usual Gaussian approximation, for comparison
## plot inla mean smooth effect...
plot.inla(pd0$times,inla,k=1,xlab="time",ylab=expression(f[1](time)))
## overlay simple Gaussian equivalent (in grey) ...
points(mcycle$times,mcycle$accel,col="grey")
lines(mcycle$times,fv$fit[,1],col="grey",lwd=2)
lines(mcycle$times,fv$fit[,1]+2*fv$se.fit[,1],lty=2,col="grey",lwd=2)
lines(mcycle$times,fv$fit[,1]-2*fv$se.fit[,1],lty=2,col="grey",lwd=2)
## same for log sd smooth...
plot.inla(pd1$times,inla,k=21,xlab="time",ylab=expression(f[2](time)))
lines(mcycle$times,fv$fit[,2],col="grey",lwd=2)
lines(mcycle$times,fv$fit[,2]+2*fv$se.fit[,2],col="grey",lty=2,lwd=2)
lines(mcycle$times,fv$fit[,2]-2*fv$se.fit[,2],col="grey",lty=2,lwd=2)
## ... notice some real differences for the log sd smooth, especially
## at the lower and upper ends of the time interval.
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ..
......@@ -15,7 +15,7 @@ and the provision of a variety of smooths of more than one variable. User define
smooths can be added. A Bayesian approach to confidence/credible interval calculation is
provided. Linear functionals of smooths, penalization of parametric model terms and linkage
of smoothing parameters are all supported. Lower level routines for generalized ridge
regression and penalized linearly constrained least squares are also available.
regression and penalized linearly constrained least squares are also available. In addition to the main modelling functions, \code{\link{jagam}} provided facilities to ease the set up of models for use with JAGS, while \code{\link{ginla}} provides marginal inference via a version of Integrated Nested Laplace Approximation.
}
\details{ \code{mgcv} provides generalized additive modelling functions \code{\link{gam}},
......
......@@ -147,7 +147,8 @@ separation of the contours of the different plots, but this is not
always easy to achieve. Note that these plots can not be modified to the same extent as the other plot.
For 2-d smooths \code{scheme==1} produces a perspective plot, while \code{scheme==2} produces a heatmap,
with overlaid contours and \code{scheme==3} a greyscale heatmap.
with overlaid contours and \code{scheme==3} a greyscale heatmap (\code{contour.col} controls the
contour colour).
Smooths of 3 and 4 variables are displayed as tiled heatmaps with overlaid contours. In the 3 variable case the third variable is discretized and a contour plot of the first 2 variables is produced for each discrete value. The panels in the lower and upper rows are labelled with the corresponding third variable value. The lowest value is bottom left, and highest at top right. For 4 variables, two of the variables are coarsely discretized and a square array of image plots is produced for each combination of the discrete values. The first two arguments of the smooth are the ones used for the image/contour plots, unless a tensor product term has 2D marginals, in which case the first 2D marginal is image/contour plotted. \code{n3} controls the number of panels.
See also \code{\link{vis.gam}}.
......
......@@ -7,27 +7,30 @@ estimates from a (RE)ML estimated \code{gam} object, provided the fit was with a
that evaluated the required Hessian.
}
\usage{
sp.vcov(x)
sp.vcov(x,edge.correct=TRUE,reg=1e-3)
}
%- maybe also `usage' for other objects documented here.
\arguments{
\item{x}{ a fitted model object of class \code{gam} as produced by \code{gam()}.}
\item{edge.correct}{ if the model was fitted with \code{edge.correct=TRUE} (see \code{\link{gam.control}}), then thereturned covariance matrix will be for the edge corrected log smoothing parameters.}
\item{reg}{regularizer for Hessian - default is equivalent to prior variance of 1000 on log smoothing parameters.}
}
\details{ Just extracts the inverse of the hessian matrix of the negative (restricted) log likelihood w.r.t
the log smoothing parameters, if this has been obtained as part of fitting.
}
\value{ A matrix corresponding to the estimated covariance matrix of the log smoothing parameter estimators,
if this can be extracted, otherwise \code{NULL}. If the scale parameter has been (RE)ML estimated (i.e. if the method was
\code{"ML"} or \code{"REML"} and the scale parameter was unknown) then the
last row and column relate to the log scale parameter.
if this can be extracted, otherwise \code{NULL}. If the scale parameter has been (RE)ML estimated (i.e. if the method was \code{"ML"} or \code{"REML"} and the scale parameter was unknown) then the
last row and column relate to the log scale parameter. If \code{edge.correct=TRUE} and this was used in fitting then the edge corrected smoothing parameters are in attribute \code{lsp} of the returned matrix.
}
\author{Simon N. Wood \email{simon.wood@r-project.org}
}
\references{
Wood, S.N. (2006) On confidence intervals for generalized additive models based on penalized regression splines. Australian and New Zealand Journal of Statistics. 48(4): 445-464.
\references{Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and
model selection for general smooth models (with discussion).
Journal of the American Statistical Association 111, 1548-1575
\url{http://dx.doi.org/10.1080/01621459.2016.1180986}