Commit b47683b8 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-12

parent 3c859201
......@@ -3,10 +3,44 @@
Currently deprecated and liable to be removed:
- bam(...,sparse=TRUE) [1.8-5]
- negbin() with search for `theta' - use 'nb' instead [1.8-0]
- single penalty tensor product smooths.
- p.type!=0 in summary.gam.
1.8-12
** "bs" B-spline smoothing basis. B-splines with derivative based penalties
of various orders.
* 'gamm' now uses a fixed scale parameter in PQL estimation for Poisson and
binomial data via the `sigma' option in lmeControl.
* bam null deviance computation was wrong with prior weights (including
binomial other than binary case), and returned deviance was wrong for
non-binary binomial. Fixed (did not affect estimation).
* improvements to "bfgs" optimizer to better deal with `infinite' smoothing
parameters.
* changed scheme=3 in default 2-D plotting to grey scale version of
scheme=2.
* 'trichol' and 'bandchol' added for banded Cholesky decompositions, plus
'sdiag' functions added for extracting and setting matrix sub- and
super- diagonals.
* p-spline constructor and Predict.matrix.pspline.smooth now allow set
up of SCOP-spline monotonic smoothers, and derivatives of smooths. Not
used in modelling functions yet.
* s(...,bs="re") now allows known precision matrix structures to be defined
using the `xt' argument of 's' see ?smooth.construct.re.smooth.spec for
details and example.
* negbin() with a grid search for `theta' is no longer supported - use
'nb' instead.
* bug fix to bam aic computation with AR rho correction.
1.8-11
......
Package: mgcv
Version: 1.8-11
Version: 1.8-12
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness
......@@ -16,6 +16,6 @@ LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2016-01-22 10:21:27 UTC; sw283
Packaged: 2016-03-02 18:21:25 UTC; sw283
Repository: CRAN
Date/Publication: 2016-01-24 00:07:54
Date/Publication: 2016-03-03 07:57:37
42c7248ef130aa0acc82c3fd89dc5d46 *ChangeLog
38d226d8a82c189fc41ea959d7550a17 *DESCRIPTION
164ed459fdb2fc5b8c362d7509285707 *ChangeLog
05cb1cd52bd24d5ec3b6e370fac7fa86 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
3e32a2a94cab5c89a107d3a0f48b02dd *NAMESPACE
3f79e768d946cfe7b3cdcf520e692d36 *R/bam.r
52ea7081e235df890d3e51198bcc3092 *R/coxph.r
c0cdd9aeb293ae9be0e46458f0255bbf *NAMESPACE
dfb4c20ba65950cf6ef13e73512a59a1 *R/bam.r
9f7dfce7213f23c0ccddd7fc31279b2d *R/coxph.r
81dedc3a4d8529ae61099024a96fdb93 *R/efam.r
e977376bf0123ebd107132bbd185f612 *R/fast-REML.r
57f8b12b97304227870a4e43a33b5d75 *R/gam.fit3.r
472d6e64f3f403aad4f6550f921a4bba *R/gam.fit4.r
693a3d54326ccd997cf63a44633b11c5 *R/gam.fit3.r
4041a125a6cc574a31d0884f9f947957 *R/gam.fit4.r
1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
5aa7c3e6ce301ae6fbd678a9b008e7c5 *R/gamlss.r
ff163969e9ad7a38857907bf38a39ec0 *R/gamm.r
3b0d5cac4a59ef1a8cb325b249699476 *R/jagam.r
8e5975e278d2405747e21a67742bf436 *R/mgcv.r
c69a68d6d3d987a61979dd73c03efc3d *R/misc.r
66c24aed2f8fc9f6bce321794f8aff87 *R/mvam.r
76192b5b36d8828af6e50013823732b4 *R/plots.r
10bc146058121c709fa7dc69419e0a65 *R/smooth.r
52292bf0d5c6e0132eea4ab68b1d2b40 *R/gamm.r
88c39fa594ea5b4105f4654f9f156d7a *R/jagam.r
535b90098a6deb5341962afe7ad82e85 *R/mgcv.r
2feca0dc9d354de7bc707c67a5891364 *R/misc.r
e64b3f08d86be17b1b641694794aebd4 *R/mvam.r
04d1a101dad47a624357982f473f49a3 *R/plots.r
9fadbd09a55b6ec6284c222840be1620 *R/smooth.r
1dde3ff2f6c3227d1a4765e41d0cf92b *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
......@@ -39,8 +39,9 @@ d0fa291cbbcef359c61e426f6ba38fbb *man/Predict.matrix.soap.film.Rd
f12468625253dd3603907de233762fd6 *man/Rrank.Rd
f50d43782718aa1ea66e0077952afc93 *man/Tweedie.Rd
80f8763baa4987579e2aa56073a9e94e *man/anova.gam.Rd
6870afb13301f2ea16c57783f557fd2b *man/bam.Rd
0d6c9d148c9fa98ae107367a796c7e6d *man/bam.Rd
3b8ce471159dd843bf26ab7eb17a9f33 *man/bam.update.Rd
cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd
f4112e262b8280c024c80ff8fa02735f *man/bug.reports.mgcv.Rd
a2beb811b1093c5e82ef32d7de1f7d32 *man/cSplineDes.Rd
a72647229fd92763636941e61597995d *man/choose.k.Rd
......@@ -56,7 +57,7 @@ e75719779e18c723ee1fd17e44e7901b *man/formXtViX.Rd
5af8941596cacadd1c54ebe2e7a83c23 *man/formula.gam.Rd
4da4d585b329769eb44f0c7a6e7dd554 *man/fs.test.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
db48a1f19f831916df2882a5bd83bf8c *man/gam.Rd
bbdaff30f65445ab481bbfdca691aca8 *man/gam.Rd
adaf0bd8e82d9472823cf3f3fa05e111 *man/gam.check.Rd
49de68e2abeb557b994032e4d7b5407a *man/gam.control.Rd
44db24b66ce63bc16d2c8bc3f5b42ac5 *man/gam.convergence.Rd
......@@ -71,7 +72,7 @@ b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
eb8648cc6b3c9374b899abd2f2b12f7b *man/gam2objective.Rd
717401fd7efa3b39d90418a5d1d0c216 *man/gamObject.Rd
a2593990d6a87f7b783d0435062dec02 *man/gamSim.Rd
8f3cacbb5e4448989aad8e5a45218bad *man/gamm.Rd
327a7fff692c5dc42862555915ae7f8d *man/gamm.Rd
f30b9dc971521416a167a6b13302f06b *man/gaulss.Rd
398a5c12285401c1d37a8edb58780bc3 *man/get.var.Rd
4f96476abbf9692f52030d3859580a31 *man/in.out.Rd
......@@ -89,12 +90,13 @@ b1c95a20afd6eb0825c00b46b8c3cbfa *man/ls.size.Rd
e5cb91b2fd8646476b8f1114228a33cf *man/mgcv-FAQ.Rd
ba14a20f6fa77f066bac7cdfe88b8fff *man/mgcv-package.Rd
6db1ae82808e56bd44c04066e2ec09aa *man/mgcv-parallel.Rd
df702cea24d0f92044a973b66a57e21f *man/missing.data.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
2f2fdc722c5e9e58664da9378451cd4a *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
324839278497102d76a0905339cf7950 *man/multinom.Rd
fafa20038af329d4d9c2e96023f7ee5c *man/mvn.Rd
887b486433499791310a39af8c3f0dbd *man/negbin.Rd
d777fbaa68f4e7b7d628e39181afb45d *man/negbin.Rd
8a6a1926188511235f1e7406120c791e *man/new.name.Rd
00e39f302ab5efbe3b14265fffc16c18 *man/notExp.Rd
7a3280b766cab8424a12d6a8b1d5748e *man/notExp2.Rd
......@@ -112,23 +114,25 @@ c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
cf14ce6cf8e4147f0f5c6e5b93b2af73 *man/print.gam.Rd
6d0ce4e574fabceffdbedd46c91364cb *man/qq.gam.Rd
f77ca1471881d2f93c74864d076c0a0e *man/rTweedie.Rd
5ff3bd5034e8f2afa3174c0a2d989d97 *man/random.effects.Rd
f834603ecaffd45099a0af864d25f691 *man/random.effects.Rd
c523210ae95cb9aaa0aaa1c37da1a4c5 *man/residuals.gam.Rd
3c747a8066bcc28ae706ccf74f903d3e *man/rig.Rd
9f6f46f5c5da080bc82f9aa4685d364a *man/rmvn.Rd
a8ad211da595840e42f71abb405b20c1 *man/s.Rd
1eb8049d1020f6ef2d6deea0aa3da435 *man/scat.Rd
bec1ae0c4e77d9b0a6534671f364d76d *man/sdiag.Rd
9641bb5e2573b5fcbdc303fcf160f3e6 *man/single.index.Rd
6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
8020154bd5c709d11f0e7cf043df886d *man/smooth.construct.Rd
4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.smooth.spec.Rd
f0311da9a0c1e3c3c087e8c377a6a3a6 *man/smooth.construct.bs.smooth.spec.Rd
76013feaf70d00976bba0154b6f2c946 *man/smooth.construct.cr.smooth.spec.Rd
f5e6d0f5122f61c336827b3615482157 *man/smooth.construct.ds.smooth.spec.Rd
db75c958cbfb561914a3291ab58b9786 *man/smooth.construct.fs.smooth.spec.Rd
e38f6ea7f89cb068976b73a177878906 *man/smooth.construct.gp.smooth.spec.Rd
4aaa84b520992fbc32b0c37f7f63c1dd *man/smooth.construct.mrf.smooth.spec.Rd
abe15377f471a2d8957a59c19eeef0bb *man/smooth.construct.ps.smooth.spec.Rd
6aaff8575f68775ed21930893ea9e03d *man/smooth.construct.re.smooth.spec.Rd
2523b6cefa306210c00f3477853b7f07 *man/smooth.construct.ps.smooth.spec.Rd
b45d8e71bda4ceca8203dffea577e441 *man/smooth.construct.re.smooth.spec.Rd
224a08b5edcd7c8af6c995f03b17320a *man/smooth.construct.so.smooth.spec.Rd
0bfe981f2c3e6ea5b8d5372076ccde53 *man/smooth.construct.sos.smooth.spec.Rd
3cb4e59f915c8d64b90754eaeeb5a86f *man/smooth.construct.t2.smooth.spec.Rd
......@@ -143,6 +147,7 @@ f0791d830687d6155efb8a73db787401 *man/summary.gam.Rd
fd71e0a218ad840e6ce4b873f9fa083e *man/t2.Rd
41dbd27ba9ae4ac0960893942880d77b *man/te.Rd
6eebb6ef90374ee09453d6da6449ed79 *man/tensor.prod.model.matrix.Rd
f22f1cee0ff2b70628846d1d0f8e9a66 *man/trichol.Rd
42542181aa314eed22f05907a8546735 *man/uniquecombs.Rd
a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
......@@ -163,15 +168,15 @@ dfd4eec9edc7d1ab6354d47b6e2bd42f *po/pl.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
b0459e16b04844271cf5e6b53aca0e47 *src/coxph.c
91c7e18bb76056ed3d89541aae8ff561 *src/discrete.c
a6b9681fae3eeddce24b3151a9442f2a *src/gdi.c
75d9ded0436a549de5b84b88c0da97e6 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
33d2b72915de9474ae236c0f3de7ca1f *src/init.c
5a8f485ee3aeba83a201bb4f136f338c *src/init.c
7c1553b521f89ab43477245c44b9d006 *src/magic.c
1c84e40793603bd3b0e3eeb6161505ca *src/mat.c
92781b80e84fbb462c6e07afaf333d9b *src/mat.c
25921d60399f6aad99a4cd761fd72628 *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
cba80885f5885bd48227696a239c53bb *src/mgcv.c
a7225f1fae5be15131c23fa20b0730ce *src/mgcv.h
46b14087c9a2d8ca5a487cae292fde13 *src/mgcv.h
b679e9063bc032f364c3ec24d57ddb08 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c
......
useDynLib(mgcv, .registration = TRUE, .fixes = "C_")
export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
cSplineDes,
exclude.too.far,extract.lme.cov, extract.lme.cov2,
formXtViX, full.score, formula.gam,fixDependence,fix.family.link,
fix.family.var, fix.family.ls, fix.family.qf,fix.family.rd,
......@@ -22,6 +23,7 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
place.knots, plot.gam, polys.plot,print.anova.gam,
print.gam,print.summary.gam,predict.gam,predict.bam,
PredictMat,Predict.matrix,Predict.matrix2,
Predict.matrix.Bspline.smooth,
Predict.matrix.cr.smooth,
Predict.matrix.duchon.spline,
Predict.matrix.cs.smooth,
......@@ -40,10 +42,11 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
Predict.matrix.gp.smooth,
qq.gam,
residuals.gam,rig,rTweedie,rmvn,
Rrank,s,scat,
Rrank,s,scat,sdiag,"sdiag<-",
sim2jam,
slanczos,
smoothCon,smooth.construct,smooth.construct2,
smooth.construct.bs.smooth.spec,
smooth.construct.cc.smooth.spec,
smooth.construct.cp.smooth.spec,
smooth.construct.cr.smooth.spec,
......@@ -65,7 +68,7 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
summary.gam,sp.vcov,
spasm.construct,spasm.sp,spasm.smooth,
t2,te,ti,tensor.prod.model.matrix,tensor.prod.penalties,
Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
trichol,Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
importFrom(grDevices,cm.colors,dev.interactive,devAskNewPage,gray,grey,heat.colors,terrain.colors,topo.colors)
importFrom(graphics,abline,axis,axTicks,box,contour,hist,image,lines,
......@@ -90,6 +93,7 @@ importFrom(nlme,Dim,corMatrix,logDet,pdConstruct,pdFactor,pdMatrix,getGroupsForm
importMethodsFrom(Matrix,t,colMeans,colSums,chol,solve,lu,expand)
importFrom(Matrix,Diagonal,sparseMatrix,Matrix)
importFrom(methods,cbind2,as)
importFrom(stats,weighted.mean)
S3method(anova, gam)
S3method(influence, gam)
......@@ -131,6 +135,7 @@ S3method(fix.family.link,family)
S3method(fix.family.link,extended.family)
S3method(fix.family.link,general.family)
S3method(smooth.construct,ad.smooth.spec)
S3method(smooth.construct,bs.smooth.spec)
S3method(smooth.construct, cc.smooth.spec)
S3method(smooth.construct,cp.smooth.spec)
S3method(smooth.construct, cr.smooth.spec)
......@@ -150,6 +155,7 @@ S3method(smooth.construct, tensor.smooth.spec)
S3method(smooth.construct, t2.smooth.spec)
S3method(smooth.construct, ts.smooth.spec)
S3method(Predict.matrix,Bspline.smooth)
S3method(Predict.matrix,cr.smooth)
S3method(Predict.matrix,cs.smooth)
S3method(Predict.matrix,cyclic.smooth)
......
......@@ -690,6 +690,7 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
object$iter <- iter
object$wt <- w
object$y <- G$y
object$prior.weights <- G$w
rm(G);if (gc.level>0) gc()
object
} ## end bgam.fitd
......@@ -1015,6 +1016,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
object$iter <- iter
object$wt <- wt
object$y <- G$y
object$prior.weights <- G$w
rm(G);if (gc.level>0) gc()
object
} ## end bgam.fit
......@@ -1184,6 +1186,7 @@ bgam.fit2 <- function (G, mf, chunk.size, gp ,scale ,gamma,method, etastart = NU
object$wt <- w
object$R <- qrx$R
object$y <- G$y
object$prior.weights <- G$w
rm(G);gc()
object
} ## end bgam.fit2
......@@ -1596,7 +1599,7 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
}
G$smooth <- G$X <- NULL
object$prior.weights <- G$w
object$AR1.rho <- rho
if (rho!=0) { ## need to store last model matrix row, to allow update
object$yX.last <- yX.last
......@@ -2209,7 +2212,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
object$nsdf <- G$nsdf
if (G$nsdf>0) names(object$coefficients)[1:G$nsdf] <- colnamesX[1:G$nsdf]
object$offset <- G$offset
object$prior.weights <- G$w
##object$prior.weights <- G$w
object$pterms <- G$pterms
object$pred.formula <- G$pred.formula
object$smooth <- G$smooth
......@@ -2243,12 +2246,13 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
sign(object$y-object$fitted.values)
if (rho!=0) object$std.rsd <- AR.resid(object$residuals,rho,object$model$"(AR.start)")
object$deviance <- sum(object$residuals^2)
object$aic <- family$aic(object$y,1,object$fitted.values,object$prior.weights,object$deviance) +
2 * (length(object$y) - sum(AR.start))*log(1/sqrt(1-rho^2)) + ## correction for AR
dev <- object$deviance <- sum(object$residuals^2)
if (rho!=0&&family$family=="gaussian") dev <- sum(object$std.rsd^2)
object$aic <- family$aic(object$y,1,object$fitted.values,object$prior.weights,dev) -
2 * (length(object$y) - sum(sum(object$model[["(AR.start)"]])))*log(1/sqrt(1-rho^2)) + ## correction for AR
2*sum(object$edf)
if (!is.null(object$edf2)&&sum(object$edf2)>sum(object$edf1)) object$edf2 <- object$edf1
object$null.deviance <- sum(family$dev.resids(object$y,mean(object$y),object$prior.weights))
object$null.deviance <- sum(family$dev.resids(object$y,weighted.mean(object$y,object$prior.weights),object$prior.weights))
if (!is.null(object$full.sp)) {
if (length(object$full.sp)==length(object$sp)&&
all.equal(object$sp,object$full.sp)==TRUE) object$full.sp <- NULL
......
## (c) Simon N. Wood (2013, 2014) coxph model extended family.
## (c) Simon N. Wood (2013, 2014) coxph model general family.
## Released under GPL2 ...
cox.ph <- function (link = "identity") {
......@@ -20,12 +20,6 @@ cox.ph <- function (link = "identity") {
env <- new.env(parent = .GlobalEnv)
validmu <- function(mu) all(is.finite(mu))
# aic <- function(y, mu, theta=NULL, wt, dev) {
# ## this needs to call coxlpl - not really enough info
# ## use store and retrieve approach, actually this is not needed (gam.fit5 computes from likelihood)
# get(".log.partial.likelihood")
# }
## initialization is tough here... need data frame in reverse time order,
......@@ -76,7 +70,8 @@ cox.ph <- function (link = "identity") {
## = censoring)
tr <- unique(y);r <- match(y,tr);nt <- length(tr)
oo <- .C("coxpp",as.double(X%*%beta),A=as.double(X),as.integer(r),d=as.integer(wt),
h=as.double(rep(0,nt)),q=as.double(rep(0,nt)),km=as.double(rep(0,nt)),n=as.integer(nrow(X)),p=as.integer(ncol(X)),
h=as.double(rep(0,nt)),q=as.double(rep(0,nt)),km=as.double(rep(0,nt)),
n=as.integer(nrow(X)),p=as.integer(ncol(X)),
nt=as.integer(nt),PACKAGE="mgcv")
p <- ncol(X)
list(tr=tr,h=oo$h,q=oo$q,a=matrix(oo$A[p*nt],p,nt),nt=nt,r=r,km=oo$km)
......@@ -116,15 +111,18 @@ cox.ph <- function (link = "identity") {
ll <- function(y,X,coef,wt,family,deriv=0,d1b=0,d2b=0,Hp=NULL,rank=0,fh=NULL,D=NULL) {
## function defining the cox model log lik.
## Calls C code "coxlpl"
## deriv codes: 0 - eval; 1 - grad and Hessian
## 2 - d1H (diagonal only)
## 3 - d1H; 4 d2H (diag)
## deriv codes: 0 - evaluate the log likelihood
## 1 - evaluate the grad and Hessian, H, of log lik w.r.t. coefs (beta)
## 2/3 - evaluate d1H =dH/drho given db/drho in d1b
## (2 is for evaluation of diagonal only)
## 4 - given d1b and d2b evaluate trHid2H= tr(Hp^{-1}d2H/drhodrho')
## Hp is the preconditioned penalized hessian of the log lik
## which is of rank 'rank'.
## fh is a factorization of Hp - either its eigen decomposition
## or its Choleski factor
## D is the diagonal pre-conditioning matrix used to obtain Hp
## if Hr is the raw Hp then Hp = D*t(D*Hr)
##tr <- sort(unique(y),decreasing=TRUE)
tr <- unique(y)
r <- match(y,tr)
......@@ -138,8 +136,6 @@ cox.ph <- function (link = "identity") {
} else M <- d1H <- 0
if (deriv > 2) {
d2H <- rep(0,p*M*(M+1)/2)
#X <- t(forwardsolve(t(L),t(X)))
#d1b <- L %*% d1b; d2b <- L %*% d2b
if (is.list(fh)) {
ev <- fh
} else { ## need to compute eigen here
......
This diff is collapsed.
......@@ -581,7 +581,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
n.theta=as.integer(length(theta)), Mp=as.integer(Mp),Enrow=as.integer(rows.E),
rSncol=as.integer(rSncol),deriv=as.integer(deriv),
fixed.penalty = as.integer(rp$fixed.penalty),nt=as.integer(control$nthreads),
type=as.integer(gdi.type))
type=as.integer(gdi.type),dVkk=as.double(rep(0,nSp^2)))
## test code used to ensure type 0 and type 1 produce identical results, when both should work.
# oot <- .C(C_gdi2,
# X=as.double(x[good,]),E=as.double(Sr),Es=as.double(Eb),rS=as.double(unlist(rS)),
......@@ -685,7 +685,8 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
scale.est=scale,reml.scale=scale,
aic=aic.model,
rank=oo$rank.est,
K=Kmat,control=control
K=Kmat,control=control,
dVkk = matrix(oo$dVkk,nSp,nSp)
#,D1=oo$D1,D2=D2,
#ldet=oo$ldet,ldet1=oo$ldet1,ldet2=ldet2,
#bSb=oo$P,bSb1=oo$P1,bSb2=bSb2,
......@@ -761,7 +762,8 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
## get log likelihood, grad and Hessian (w.r.t. coefs - not s.p.s) ...
ll <- family$ll(y,x,coef,weights,family,deriv=1)
llf <- family$ll
ll <- llf(y,x,coef,weights,family,deriv=1)
ll0 <- ll$l - (t(coef)%*%St%*%coef)/2
rank.checked <- FALSE ## not yet checked the intrinsic rank of problem
rank <- q;drop <- NULL
......@@ -776,7 +778,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
fdg <- ll$lb*0; fdh <- ll$lbb*0
for (k in 1:length(coef)) {
coef1 <- coef;coef1[k] <- coef[k] + eps
ll.fd <- family$ll(y,x,coef1,weights,family,deriv=1)
ll.fd <- llf(y,x,coef1,weights,family,deriv=1)
fdg[k] <- (ll.fd$l-ll$l)/eps
fdh[,k] <- (ll.fd$lb-ll$lb)/eps
}
......@@ -833,15 +835,15 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
}
## try the Newton step...
coef1 <- coef + step
ll <- family$ll(y,x,coef1,weights,family,deriv=1)
ll <- llf(y,x,coef1,weights,family,deriv=1)
ll1 <- ll$l - (t(coef1)%*%St%*%coef1)/2
khalf <- 0;fac <- 2
while (ll1 < ll0 && khalf < 25) { ## step halve until it succeeds...
step <- step/fac;coef1 <- coef + step
ll <- family$ll(y,x,coef1,weights,family,deriv=0)
ll <- llf(y,x,coef1,weights,family,deriv=0)
ll1 <- ll$l - (t(coef1)%*%St%*%coef1)/2
if (ll1>=ll0) {
ll <- family$ll(y,x,coef1,weights,family,deriv=1)
ll <- llf(y,x,coef1,weights,family,deriv=1)
} else { ## abort if step has made no difference
if (max(abs(coef1-coef))==0) khalf <- 100
}
......@@ -856,10 +858,10 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
while (ll1 < ll0 && khalf < 25) { ## step cut until it succeeds...
step <- step/10;coef1 <- coef + step
ll <- family$ll(y,x,coef1,weights,family,deriv=0)
ll <- llf(y,x,coef1,weights,family,deriv=0)
ll1 <- ll$l - (t(coef1)%*%St%*%coef1)/2
if (ll1>=ll0) {
ll <- family$ll(y,x,coef1,weights,family,deriv=1)
ll <- llf(y,x,coef1,weights,family,deriv=1)
} else { ## abort if step has made no difference
if (max(abs(coef1-coef))==0) khalf <- 100
}
......@@ -878,7 +880,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
perturbed <- perturbed + 1
coef <- coef*(1+(runif(length(coef))*.02-.01)*perturbed) +
(runif(length(coef)) - 0.5 ) * mean(abs(coef))*1e-5*perturbed
ll <- family$ll(y,x,coef,weights,family,deriv=1)
ll <- llf(y,x,coef,weights,family,deriv=1)
ll0 <- ll$l - (t(coef)%*%St%*%coef)/2
} else {
rank.checked <- TRUE
......@@ -914,7 +916,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
} ## lpi adjustment done
attr(x,"lpi") <- lpi
attr(x,"drop") <- drop ## useful if family has precomputed something from x
ll <- family$ll(y,x,coef,weights,family,deriv=1)
ll <- llf(y,x,coef,weights,family,deriv=1)
ll0 <- ll$l - (t(coef)%*%St%*%coef)/2
}
}
......@@ -943,7 +945,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
fcoef <- rep(0,length(bdrop));fcoef[!bdrop] <- coef
} else fcoef <- coef
d1l <- d2l <- d1bSb <- d2bSb <- d1b <- d2b <- d1ldetH <- d2ldetH <- d1b <- d2b <- NULL
dVkk <- d1l <- d2l <- d1bSb <- d2bSb <- d1b <- d2b <- d1ldetH <- d2ldetH <- d1b <- d2b <- NULL
if (deriv>0) { ## Implicit differentiation for derivs...
......@@ -953,6 +955,9 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
if (nSp) for (i in 1:m) d1b[,i] <-
-D*(backsolve(L,forwardsolve(t(L),(D*Sib[[i]][!bdrop])[piv]))[ipiv])
## obtain the curvature check matrix...
dVkk <- crossprod(L[,ipiv]%*%(d1b/D))
if (!is.null(drop)) { ## create full version of d1b with zeros for unidentifiable
fd1b <- matrix(0,q,m)
fd1b[!bdrop,] <- d1b
......@@ -961,7 +966,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
## Now call the family again to get first derivative of Hessian w.r.t
## smoothing parameters, in list d1H...
ll <- family$ll(y,x,coef,weights,family,deriv=3,d1b=d1b)
ll <- llf(y,x,coef,weights,family,deriv=3,d1b=d1b)
d1l <- colSums(ll$lb*d1b)
......@@ -978,7 +983,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
## Now call family for last time to get trHid2H the tr(H^{-1} d^2 H / drho_i drho_j)...
llr <- family$ll(y,x,coef,weights,family,deriv=4,d1b=d1b,d2b=d2b,
llr <- llf(y,x,coef,weights,family,deriv=4,d1b=d1b,d2b=d2b,
Hp=Hp,rank=rank,fh = L,D=D)
## Now compute Hessian of log lik w.r.t. log sps using chain rule
......@@ -1095,7 +1100,12 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
#S=rp$ldetS,S1=rp$ldet1,S2=rp$ldet2,
#Hp=ldetHp,Hp1=d1ldetH,Hp2=d2ldetH,
#b2 = d2b)
H = ll$lbb,dH = ll$d1H)#,d2H=llr$d2H)
H = ll$lbb,dH = ll$d1H,dVkk=dVkk)#,d2H=llr$d2H)
## debugging code to allow components of 2nd deriv of hessian w.r.t. sp.s
## to be passed to deriv.check....
#if (!is.null(ll$ghost1)&&!is.null(ll$ghost2)) {
# ret$ghost1 <- ll$ghost1; ret$ghost2 <- ret$ghost2
#}
ret
} ## end of gam.fit5
......
......@@ -1163,6 +1163,11 @@ gammPQL <- function (fixed, random, family, data, correlation, weights,
converged <- FALSE
if (family$family %in% c("poisson","binomial")) {
control$sigma <- 1 ## set scale parameter to 1
control$apVar <- FALSE ## not available
}
for (i in 1:niter) {
if (verbose) message(gettextf("iteration %d", i))
fit <- lme(fixed=fixed,random=random,data=data,correlation=correlation,
......@@ -1609,6 +1614,7 @@ test.gamm <- function(control=nlme::lmeControl(niterEM=3,tolerance=1e-11,msTol=1
x<-runif(n)/20;z<-runif(n);
f <- test1(x,z)
y <- f + rnorm(n)*0.2
control$sigma <- NULL ## avoid failure on silly test
cat("testing covariate scale invariance ... ")
b <- gamm(y~te(x,z), control=control )
x1 <- x*100
......
......@@ -3,7 +3,7 @@
## Code offering JAGS/BUGS support for mgcv.
## In particular autogenerates the code and data to fit an mgcv
## style GAM in JAGS, and re-packages the simulation output
## in a form suitable for plotting and predcition.
## in a form suitable for plotting and prediction.
## Idea is that the code would be modified to add the sort
## of random effects structure most appropriately handled in JAGS.
......
This diff is collapsed.
......@@ -18,6 +18,82 @@ rmvn <- function(n,mu,V) {
z
} ## rmvn
sdiag <- function(A,k=0) {
## extract sub or super diagonal of matrix (k=0 is leading)
p <- ncol(A)
n <- nrow(A)
if (k>p-1||-k > n-1) return()
if (k >= 0) {
i <- 1:n
j <- (k+1):p
} else {
i <- (-k+1):n
j <- 1:p
}
if (length(i)>length(j)) i <- i[1:length(j)] else j <- j[1:length(i)]
ii <- i + (j-1) * n
A[ii]
} ## sdiag
"sdiag<-" <- function(A,k=0,value) {
p <- ncol(A)
n <- nrow(A)
if (k>p-1||-k > n-1) return()
if (k >= 0) {
i <- 1:n
j <- (k+1):p
} else {
i <- (-k+1):n
j <- 1:p
}
if (length(i)>length(j)) i <- i[1:length(j)] else j <- j[1:length(i)]
ii <- i + (j-1) * n
A[ii] <- value
A
} ## "sdiag<-"
bandchol <- function(B) {
## obtain R such that R'R = A. Where A is banded matrix contained in R.
n <- ncol(B)
k <- 0
if (n==nrow(B)) { ## square matrix. Extract the diagonals
A <- B*0
for (i in 1:n) {
b <- sdiag(B,i-1)
if (sum(b!=0)!=0) {
k <- i ## largest index of a non-zero band
A[i,1:length(b)] <- b
}
}
B <- A[1:k,]
}
oo <- .C(C_band_chol,B=as.double(B),n=as.integer(n),k=as.integer(nrow(B)),info=as.integer(0))
if (oo$info<0) stop("something wrong with inputs to LAPACK routine")
if (oo$info>0) stop("not positive definite")
B <- matrix(oo$B,nrow(B),n)
if (k>0) { ## was square on entry, so also on exit...
A <- A * 0
for (i in 1:k) sdiag(A,i-1) <- B[i,1:(n-i+1)]
B <- A
}
B
} ## bandchol
trichol <- function(ld,sd) {
## obtain chol factor R of symm tridiag matrix, A, with leading diag
## ld and sub/super diags sd. R'R = A. On exit ld is diag of R and
## sd its super diagonal.
n <- length(ld)
if (n<2) stop("don't be silly")
if (n!=length(sd)+1) stop("sd should have exactly one less entry than ld")
oo <- .C(C_tri_chol,ld=as.double(ld),sd=as.double(sd),n=as.integer(n),info=as.integer(0))
if (oo$info<0) stop("something wrong with inputs to LAPACK routine")
if (oo$info>0) stop("not positive definite")
ld <- sqrt(oo$ld)
sd <- oo$sd*ld[1:(n-1)]
list(ld=ld,sd=sd)
}
mgcv.omp <- function() {
## does open MP appear to be available?
oo <- .C(C_mgcv_omp,a=as.integer(-1))
......
......@@ -160,9 +160,11 @@ mvn <- function(d=2) {
ll <- function(y,X,coef,wt,family,deriv=0,d1b=NULL,d2b=NULL,Hp=NULL,rank=0,fh=NULL,D=NULL) {
## function defining the Multivariate Normal model log lik.
## Calls C code "mvn_ll"
## deriv codes: 0 - eval; 1 - grad and Hessian
## 2 - d1H (diagonal only - not implemented efficiently)
## 3 - d1H; 4 d2H (diag - not implemented)
## deriv codes: 0 - evaluate the log likelihood
## 1 - evaluate the grad and Hessian, H, of log lik w.r.t. coefs (beta)
## 2/3 - evaluate d1H =dH/drho given db/drho in d1b
## (2 is diagonal only - not implemented efficiently)
## 4 - given d1b and d2b evaluate trHid2H= tr(Hp^{-1}d2H/drhodrho') (not implemented)
## Hp is the preconditioned penalized hessian of the log lik
## which is of rank 'rank'.
## fh is a factorization of Hp - either its eigen decomposition
......
......@@ -909,10 +909,11 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (scheme == 1) { ## perspective plot
persp(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
zlab=P$main,ylim=P$ylim,xlim=P$xlim,theta=theta,phi=phi,...)
} else if (scheme==2) {
} else if (scheme==2||scheme==3) {
if (scheme==3) colors <- grey(0:50/50)
image(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
main=P$main,xlim=P$xlim,ylim=P$ylim,col=colors,...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=3,...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=contour.col,...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
......@@ -954,7 +955,8 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (scheme==1) {
persp(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
zlab=P$main,theta=theta,phi=phi,xlim=P$xlim,ylim=P$ylim,...)
} else if (scheme==2) {
} else if (scheme==2||scheme==3) {
if (scheme==3) colors <- grey(0:50/50)
image(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
main=P$main,xlim=P$xlim,ylim=P$ylim,col=colors,...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=contour.col,...)
......
This diff is collapsed.
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
......@@ -60,8 +60,8 @@ which \code{gam} is called.}
\item{offset}{Can be used to supply a model offset for use in fitting. Note
that this offset will always be completely ignored when predicting, unlike an offset
included in \code{formula}: this conforms to the behaviour of
\code{lm} and \code{glm}.}
included in \code{formula} (this used to conform to the behaviour of
\code{lm} and \code{glm}).}
......
\name{bandchol}
\alias{bandchol}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Choleski decomposition of a band diagonal matrix}
\description{
Computes Choleski decomposition of a (symmetric positive definite) band-diagonal matrix, \code{A}.
}