Commit ee00362c authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-7

parent 4079654a
** denotes quite substantial/important changes
*** denotes really big changes
1.8-7
** 'gam' default scale parameter changed to modified Pearson estimator
developed by Fletcher 2012 Biometrika 99(1), 230-237. See ?gam.scale.
** 'bam' now has a 'discrete' argument to allow discretization of covariates
for more efficient computation, with substantially more parallelization
(via 'nthreads'). Still somewhat experimental.
* Slightly more accurate smoothing parameter uncertainty correction. Changes
edf2 used for AIC (under RE/ML), and hence may change AIC values.
* jagam prior variance on fixed effects is now set with reference to data and
model during initialization step.
* bam could lose offset for small datasets in gaussian additive case. fixed.
* gam.side now setup to include penalties in computations if fewer data than
coefs (an exceedingly specialist topic).
* p-value computation for smooth terms modified to avoid an ambiguity in the
choice of test statistic that could lead to p-value changing somewhat between
platforms.
* gamm now warns if attempt is made to use extended family.
* step fail logic improved for "fREML" optimization in 'bam'.
* fix of openMP error in mgcv_pbsi, which could cause a problem in
multi-threaded bam computation (failure to declare a variable as private).
* Smoothing parameter uncertainty corrected AIC calculations had an
indexing problem in Sl.postproc, which could result in failure of bam with
linked smooths.
* mroot patched for fact that chol(...,pivot=TRUE) does not operate as
documented on rank deficient matrices: trailing block of triangular factor
has to be zeroed for pivoted crossprod of factor to equal original matrix.
* bam(...,sparse=TRUE) deprecated as no examples found where it is really
worthwhile (but let me know if this is a problem).
* marginal model matrices in tensor product smooths now stored in
re-parameterized form, if re-parameterization happened (shouldn't change
anything!).
* initial.spg could fail if response vector had dim attributes and extended
family used. fixed.
1.8-6
* Generalization of list formula handling to allow linear predictors to
......
Package: mgcv
Version: 1.8-6
Version: 1.8-7
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
......@@ -15,7 +15,7 @@ Suggests: splines, parallel, survival, MASS
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2015-03-30 13:02:31 UTC; sw283
NeedsCompilation: yes
Packaged: 2015-07-22 17:38:22 UTC; sw283
Repository: CRAN
Date/Publication: 2015-03-31 12:46:38
Date/Publication: 2015-07-23 07:05:08
2f31dc4982b7fa7c7f08ff3fad94b88d *ChangeLog
2244e729b443669c81a2732d005dd044 *DESCRIPTION
f876a44fccfab3bcf4e8c6a8916f480f *ChangeLog
68b8a0864224d96ff69d755507862f24 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
4d14a8a6cdbe2e25634829667d3251ec *NAMESPACE
bb90db9c74f22d202abc4b0fc562448e *R/bam.r
053c3c177572e1570e63329445eacc14 *NAMESPACE
e99642616897453fa1cab582bc73c84f *R/bam.r
52ea7081e235df890d3e51198bcc3092 *R/coxph.r
a53c7592bb4c10a961f58e6537c87355 *R/efam.r
8cde8b6379034d9eb78caaf1b0d04cf8 *R/fast-REML.r
5b6cfb632772fb514428730b22066480 *R/gam.fit3.r
fc965ed3d5edf72c6f0335347c383dbb *R/gam.fit4.r
209a7a1d10cec589213ffed9efc1d5cc *R/efam.r
9131483dd63ef2b80e8159601cc03922 *R/fast-REML.r
88dbdd15671f0a924bdd9f5f42a89a34 *R/gam.fit3.r
a6208b52c8af1f8779dd0b4b482c3fa7 *R/gam.fit4.r
e63276b78a4ff2566af2e651bfc6aa9c *R/gam.sim.r
633c048e36083646c7d409e99cf7b037 *R/gamlss.r
28ddc8624513070a62d0cc2aba7776a2 *R/gamm.r
09240f99d77e54848bf15d540c267022 *R/jagam.r
70eb07e6f5c4d61f6519863b8eebdfc0 *R/mgcv.r
bf5f3e43220b1e6db57c80caf66f1e54 *R/misc.r
ff163969e9ad7a38857907bf38a39ec0 *R/gamm.r
3b0d5cac4a59ef1a8cb325b249699476 *R/jagam.r
d5ce48c9ea98c2f428cfb00bcc9fb5ea *R/mgcv.r
59a92561754c8156506af4004f594ac9 *R/misc.r
66c24aed2f8fc9f6bce321794f8aff87 *R/mvam.r
f031621e351c06b12c328e4033e2f097 *R/plots.r
3e87a7ef880c8aa40ce57201f6ea5458 *R/smooth.r
966c267722f1027c53aa5621614333ae *R/smooth.r
68348617c5d3b0e2ea805c4360a8cdd4 *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
......@@ -32,8 +32,6 @@ c745969a1292eb3d49dfd0d0c2c997d4 *inst/po/de/LC_MESSAGES/mgcv.mo
c1b1475e5fef49fe49929d2796ff87b6 *inst/po/ko/LC_MESSAGES/mgcv.mo
cd7e6d1282796c089c320fbff388047f *inst/po/pl/LC_MESSAGES/R-mgcv.mo
715e52c0debf9848bbda15e94f5e7315 *inst/po/pl/LC_MESSAGES/mgcv.mo
6224353bb8556c6ae4106352e244a3ab *inst/po/po/LC_MESSAGES/R-mgcv.mo
4e7e8faae00111b7db61daab04358202 *inst/po/po/LC_MESSAGES/mgcv.mo
c574fe1ca9d55a9818d308906f16d16e *man/Beta.Rd
5bf12ddc0dab9daae72271b96a15c539 *man/Predict.matrix.Rd
c45c0f78f753461b33a295883461e732 *man/Predict.matrix.cr.smooth.Rd
......@@ -41,7 +39,7 @@ e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
f12468625253dd3603907de233762fd6 *man/Rrank.Rd
f6cadda5999c35800fd65c08d6812f7b *man/Tweedie.Rd
80f8763baa4987579e2aa56073a9e94e *man/anova.gam.Rd
de9b8ca4369a4fc7ba69c68c1cdb98c2 *man/bam.Rd
e49e38bc1741829ae997f8b0b4fcd832 *man/bam.Rd
71bde8b8caa24a36529ce7e0ac3165d8 *man/bam.update.Rd
a2beb811b1093c5e82ef32d7de1f7d32 *man/cSplineDes.Rd
9b4d616d1b6c4a46ca77d16cded3f806 *man/choose.k.Rd
......@@ -59,13 +57,13 @@ c6fd48d86959402982c75566876baa16 *man/formula.gam.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
86678ba1579e9deb8ac3337af60a9c20 *man/gam.Rd
fe61dd0efab9e920c17335faf3d5764c *man/gam.check.Rd
a65bc22f606e45d185bc375fbf5698a1 *man/gam.control.Rd
ca99cf1e23fd828e07b49a3588d35d65 *man/gam.control.Rd
44db24b66ce63bc16d2c8bc3f5b42ac5 *man/gam.convergence.Rd
58ab3b3d6f4fd0d008d73c3c4e6d3305 *man/gam.fit.Rd
dcf10ab3cc3102f7fb36d3ddf44013f5 *man/gam.fit3.Rd
8ba3991b5932b0775b452d20c9ff4d54 *man/gam.models.Rd
e969287d1a5c281faa7eb6cfce31a7c5 *man/gam.outer.Rd
39e2d1fa5374f3b5cff1d761b73f0daa *man/gam.scale.Rd
c17814cea1b11e5ca374e72d6e1cbd98 *man/gam.scale.Rd
96676186808802344a99f9d3170bf775 *man/gam.selection.Rd
956bf1a6ac1361dd0403c28153b03a9b *man/gam.side.Rd
b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
......@@ -78,9 +76,9 @@ df4a6db696749986fd5e20586fc9b718 *man/gaulss.Rd
a2ea1233d43fac89e0cacbc09a8d31e2 *man/in.out.Rd
9c461959be1272edcb98ee7e20fdc317 *man/inSide.Rd
2f222eeeb3d7bc42f93869bf8c2af58a *man/influence.gam.Rd
bafea2eef12fdc819f8ac1fb41d8b914 *man/initial.sp.Rd
39b9de9dbac7d9dc5c849e1a37def675 *man/initial.sp.Rd
c00bcfe2d0b44b2ea955f3934421807c *man/interpret.gam.Rd
469e313e7a1b0520593c8e7a938111b5 *man/jagam.Rd
b440c31b77c368b484d9665b2b3e89fb *man/jagam.Rd
07d2c259b9edf164f42935170b4fccd0 *man/ldTweedie.Rd
58e73ac26b93dc9d28bb27c8699e12cf *man/linear.functional.terms.Rd
93035193b0faa32700e1421ce8c1e9f6 *man/logLik.gam.Rd
......@@ -147,32 +145,31 @@ a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
2c5f6815e609f2cdb56b0067a183f915 *man/ziP.Rd
ae8388103d8b1be39f55f426b205b576 *man/ziplss.Rd
4550c3e06d76e6e8a9b22fb8bbd76eb3 *po/R-de.po
7bd0744ad8ea562d7a624e066ef3390c *po/R-de.po
0bdfcf98961b0d52b60f806dc1dba77e *po/R-en@quot.po
4e65e93fef4d034a399f90421e8f323a *po/R-fr.po
73cdaf7a5a69f0b7cbfe411cd0c468b6 *po/R-ko.po
7eb472ce4c2d1dc30b4dd1091c0e88df *po/R-mgcv.pot
7b07899266c3acf3d2a625850d7cd6ef *po/R-pl.po
5b91cecd0b9e52154185d380a273c623 *po/R-po.po
ccf3140169b68ec7aff4b15e6a97e5db *po/de.po
382c94188dbc193fca9628287b66d1af *po/de.po
93f72334356fe6f05a64e567efd35c8e *po/en@quot.po
fb829b82760779929951d49fe29ed2e5 *po/fr.po
dc1ef92ff4454734c3a24876e299b760 *po/ko.po
8ad4757e026d1841c8f43eb97072c06e *po/mgcv.pot
dfd4eec9edc7d1ab6354d47b6e2bd42f *po/pl.po
b56dac4547037ea45e2c8f9bce7aa9ef *po/po.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
e71c1a1624b431fbab0a4c8f151d2a97 *src/coxph.c
08e156711c686a1f1efe505d63fabef5 *src/gdi.c
48f0c9de0da60d48cf4b306c2fdd039a *src/discrete.c
e4236680abfb9c89c7bf09b441c755c2 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
e5bf24371be5ea7f3ffb40060a648803 *src/init.c
890105488631271ad8b184aa0524b59f *src/init.c
9a5d7cb3cf93cbdfc08353fbd20a270e *src/magic.c
a9735df73ae117df1b4c1197f4748389 *src/mat.c
55b92c3ab6a037f9742542c2af8cad56 *src/mat.c
aae0f298e384952bb8b6b923d40520a8 *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
ec2ae157a9e3bedcccc88b053d0a4e0c *src/mgcv.c
b62374abb9910e539a3280d66ac4193f *src/mgcv.h
00f8a024faef17f90ed04f01e736df71 *src/misc.c
39b22a8fa2128fd6b52b3601841b0393 *src/mgcv.h
6b90633f745d3b3314ec5a784bde99d0 *src/misc.c
08c1706ffeec4277c484435a0644b0e3 *src/mvn.c
cbe54250deb38aa5f88f8b669a4468cd *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h
......
......@@ -65,15 +65,29 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
t2,te,ti,tensor.prod.model.matrix,tensor.prod.penalties,
Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
importFrom(grDevices,cm.colors,gray,heat.colors,terrain.colors,topo.colors)
importFrom(graphics,axis,box,contour,hist,lines,mtext, par, persp,plot,points,
polygon,strheight,strwidth,text)
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,
mtext, par, persp,plot,points,
polygon,rect,strheight,strwidth,text,title)
importFrom(stats,.checkMFClasses,.getXlevels,anova,approx,as.formula,
binomial,coef,cooks.distance,cor,cov,
delete.response,dbeta,dgamma,dnorm,dpois,fitted,formula,gaussian,glm,
influence,logLik,lm,mad,
make.link,median,model.frame,model.offset,model.matrix,nlm,
na.pass,napredict,na.omit,naresid,optim,pchisq,pnorm,pt,pf,
power,predict,printCoefmat,quantile,
qbeta,qbinom,qchisq,qnbinom,qgamma,qnorm,qpois,qqline,qqnorm,qqplot,
reformulate,residuals,
rbeta,rbinom,rgamma,rnbinom,rnorm,rpois,runif,sd,
termplot,terms.formula,terms,uniroot,var,vcov,weights)
importFrom(utils,object.size)
importFrom(stats,anova,influence,cooks.distance,logLik,vcov,residuals,predict,model.matrix)
importFrom(nlme,Dim,corMatrix,logDet,pdConstruct,pdFactor,pdMatrix,getGroupsFormula,lme,varFixed,lmeControl)
importMethodsFrom(Matrix,t,colMeans,colSums,chol,solve,lu,expand)
importFrom(Matrix,Diagonal,sparseMatrix,Matrix)
importFrom(methods,cbind2)
importFrom(methods,cbind2,as)
S3method(anova, gam)
S3method(influence, gam)
......
This diff is collapsed.
......@@ -96,11 +96,11 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
validmu <- function(mu) all(is.finite(mu))
dev.resids <- function(y, mu, wt,theta=NULL) {
F <- function(x) {
h <- ind <- x > 0; h[ind] <- 1/(exp(-x[ind]) + 1)
x <- exp(x[!ind]); h[!ind] <- (x/(1+x))
h
}
#F <- function(x) {
# h <- ind <- x > 0; h[ind] <- 1/(exp(-x[ind]) + 1)
# x <- exp(x[!ind]); h[!ind] <- (x/(1+x))
# h
#}
Fdiff <- function(a,b) {
## cancellation resistent F(b)-F(a), b>a
h <- rep(1,length(b)); h[b>0] <- -1; eb <- exp(b*h)
......@@ -144,11 +144,11 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
Dd <- function(y, mu, theta, wt=NULL, level=0) {
## derivatives of the deviance...
F <- function(x) { ## e^(x)/(1+e^x) without overflow
h <- ind <- x > 0; h[ind] <- 1/(exp(-x[ind]) + 1)
x <- exp(x[!ind]); h[!ind] <- (x/(1+x))
h
}
# F <- function(x) { ## e^(x)/(1+e^x) without overflow
# h <- ind <- x > 0; h[ind] <- 1/(exp(-x[ind]) + 1)
# x <- exp(x[!ind]); h[!ind] <- (x/(1+x))
# h
# }
Fdiff <- function(a,b) {
## cancellation resistent F(b)-F(a), b>a
h <- rep(1,length(b)); h[b>0] <- -1; eb <- exp(b*h)
......@@ -958,7 +958,7 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
## derivatives of the -2*loglik...
## ltheta <- theta
theta <- exp(theta)
onemu <- 1 - mu; oney <- 1 - y
onemu <- 1 - mu; ## oney <- 1 - y
muth <- mu*theta; ## yth <- y*theta
onemuth <- onemu*theta ## (1-mu)*theta
psi0.th <- digamma(theta)
......@@ -1595,18 +1595,18 @@ ziP <- function (theta = NULL, link = "identity") {
})
fv <- function(lp,theta=NULL) {
## optional function to give fitted values...
if (is.null(theta)) theta <- get(".Theta")
th1 <- theta[1]; th2 <- exp(theta[2]);
eta <- th1 + th2*lp
p <- 1 - exp(-exp(eta))
fv <- lambda <- exp(lp)
ind <- lp < log(.Machine$double.eps)/2
fv[!ind] <- p[!ind] * lambda[!ind]/(1-exp(-lambda[!ind]))
fv[ind] <- p[ind]
fv
} ## fv
# fv <- function(lp,theta=NULL) {
# ## optional function to give fitted values...
# if (is.null(theta)) theta <- get(".Theta")
# th1 <- theta[1]; th2 <- exp(theta[2]);
# eta <- th1 + th2*lp
# p <- 1 - exp(-exp(eta))
# fv <- lambda <- exp(lp)
# ind <- lp < log(.Machine$double.eps)/2
# fv[!ind] <- p[!ind] * lambda[!ind]/(1-exp(-lambda[!ind]))
# fv[ind] <- p[ind]
# fv
# } ## fv
rd <- function(mu,wt,scale) {
## simulate data given fitted latent variable in mu
......
This diff is collapsed.
This diff is collapsed.
......@@ -2,7 +2,6 @@
## Routines for gam estimation beyond exponential family.
dDeta <- function(y,mu,wt,theta,fam,deriv=0) {
## What is available directly from the family are derivatives of the
## deviance and link w.r.t. mu. This routine converts these to the
......@@ -535,7 +534,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
Det2=as.double(dd$Deta2),Dth2=as.double(dd$Dth2),Det.th=as.double(dd$Detath),
Det2.th=as.double(dd$Deta2th),Det3=as.double(dd$Deta3),Det.th2 = as.double(dd$Detath2),
Det4 = as.double(dd$Deta4),Det3.th=as.double(dd$Deta3th), Deta2.th2=as.double(dd$Deta2th2),
beta=as.double(coef),b1=as.double(rep(0,ntot*ncol(x))),
beta=as.double(coef),b1=as.double(rep(0,ntot*ncol(x))),w1=rep(0,ntot*length(z)),
D1=as.double(rep(0,ntot)),D2=as.double(rep(0,ntot^2)),
P=as.double(0),P1=as.double(rep(0,ntot)),P2 = as.double(rep(0,ntot^2)),
ldet=as.double(1-2*(scoreType=="ML")),ldet1 = as.double(rep(0,ntot)),
......@@ -551,6 +550,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
rV <- T %*% rV
## derivatives of coefs w.r.t. sps etc...
db.drho <- if (deriv) T %*% matrix(oo$b1,ncol(x),ntot) else NULL
dw.drho <- if (deriv) matrix(oo$w1,length(z),ntot) else NULL
Kmat <- matrix(0,nrow(x),ncol(x))
Kmat[good,] <- oo$X ## rV%*%t(K)%*%(sqrt(wf)*X) = F; diag(F) is edf array
......@@ -607,6 +607,11 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
ww <- wt <- rep.int(0, nobs)
wt[good] <- wf
ww[good] <- w
if (deriv && nrow(dw.drho)!=nrow(x)) {
w1 <- dw.drho
dw.drho <- matrix(0,nrow(x),ncol(w1))
dw.drho[good,] <- w1
}
aic.model <- family$aic(y, mu, theta, weights, dev) # note: incomplete 2*edf needs to be added
......@@ -619,7 +624,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
df.null = nulldf, y = y, converged = conv,
boundary = boundary,
REML=REML,REML1=REML1,REML2=REML2,
rV=rV,db.drho=db.drho,
rV=rV,db.drho=db.drho,dw.drho=dw.drho,
scale.est=scale,reml.scale=scale,
aic=aic.model,
rank=oo$rank.est,
......@@ -700,7 +705,7 @@ 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)
ll0 <- ll$l - t(coef)%*%St%*%coef/2
ll0 <- ll$l - (t(coef)%*%St%*%coef)/2
rank.checked <- FALSE ## not yet checked the intrinsic rank of problem
rank <- q;drop <- NULL
eigen.fix <- FALSE
......@@ -773,16 +778,36 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
coef1 <- coef + step
ll <- family$ll(y,x,coef1,weights,family,deriv=1)
ll1 <- ll$l - (t(coef1)%*%St%*%coef1)/2
khalf <- 0
while (ll1 < ll0 && khalf < 50) { ## step halve until it succeeds...
step <- step/2;coef1 <- coef + step
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)
ll1 <- ll$l - (t(coef1)%*%St%*%coef1)/2
if (ll1>=ll0) {
ll <- family$ll(y,x,coef1,weights,family,deriv=1)
} else { ## abort if step has made no difference
if (max(abs(coef1-coef))==0) khalf <- 100
}
khalf <- khalf + 1
if (khalf>5) fac <- 5
} ## end step halve
if (ll1 < ll0) { ## switch to steepest descent...
step <- -.5*drop(grad)*mean(abs(coef))/mean(abs(grad))
khalf <- 0
}
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)
ll1 <- ll$l - (t(coef1)%*%St%*%coef1)/2
if (ll1>=ll0) {
ll <- family$ll(y,x,coef1,weights,family,deriv=1)
} else { ## abort if step has made no difference
if (max(abs(coef1-coef))==0) khalf <- 100
}
khalf <- khalf + 1
}
if (ll1 >= ll0||iter==control$maxit) { ## step ok. Accept and test
coef <- coef + step
......@@ -1017,7 +1042,7 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
ret
} ## end of gam.fit5
gam.fit5.post.proc <- function(object,Sl,L) {
gam.fit5.post.proc <- function(object,Sl,L,S,off) {
## object is object returned by gam.fit5, Sl is penalty object, L maps working sp
## vector to full sp vector
## Computes:
......@@ -1089,13 +1114,27 @@ gam.fit5.post.proc <- function(object,Sl,L) {
}
## compute the smoothing parameter uncertainty correction...
if (!is.null(object$outer.info$hess)) {
ev <- eigen(object$outer.info$hess,symmetric=TRUE)
ind <- ev$values <= 0
ev$values[ind] <- 0;ev$values[!ind] <- 1/sqrt(ev$values[!ind])
if (!is.null(object$outer.info$hess)) {
if (!is.null(L)) object$db.drho <- object$db.drho%*%L ## transform to derivs w.r.t. working
Vc <- crossprod((ev$values*t(ev$vectors))%*%t(object$db.drho))
Vc <- Vb + Vc ## Bayesian cov matrix with sp uncertainty
ev <- eigen(object$outer.info$hess,symmetric=TRUE)
d <- ev$values;ind <- d <= 0
d[ind] <- 0;d[!ind] <- 1/sqrt(d[!ind])
Vc <- crossprod((d*t(ev$vectors))%*%t(object$db.drho))
#dpv <- rep(0,ncol(object$outer.info$hess));M <- length(off)
#dpv[1:M] <- 1/100 ## prior precision (1/var) on log smoothing parameters
#Vr <- chol2inv(chol(object$outer.info$hess + diag(dpv,ncol=length(dpv))))[1:M,1:M]
#Vc <- object$db.drho%*%Vr%*%t(object$db.drho)
#dpv[1:M] <- 1/10 ## prior precision (1/var) on log smoothing parameters
#Vr <- chol2inv(chol(object$outer.info$hess + diag(dpv,ncol=length(dpv))))[1:M,1:M]
#M <- length(off)
d <- ev$values; d[ind] <- 0;
d <- d + 1/50 #d[1:M] <- d[1:M] + 1/50
d <- 1/sqrt(d)
Vr <- crossprod(d*t(ev$vectors))
#Vc2 <- Vb.corr(R,L,S,off,dw=NULL,w=NULL,log(object$sp),Vr)
Vc <- Vb + Vc #+ Vc2 ## Bayesian cov matrix with sp uncertainty
## reverse the various re-parameterizations...
} else Vc <- Vb
Vc <- Sl.repara(object$rp,Vc,inverse=TRUE)
......@@ -1112,8 +1151,14 @@ gam.fit5.post.proc <- function(object,Sl,L) {
## model. This is larger than edf2 should be, because of bias correction variability,
## but is bounded in a way that is not *guaranteed* for edf2. Note that
## justification only applies to sum(edf1/2) not elementwise
if (!is.null(object$outer.info$hess)) {
## second correction term is easier computed in original parameterization...
Vc2 <- Vb.corr(R,L,S,off,dw=NULL,w=NULL,log(object$sp),Vr)
Vc <- Vc + Vc2
}
edf1 <- 2*edf - rowSums(t(F)*F)
edf2 <- diag(Vc%*%crossprod(R))
#edf2 <- diag(Vc%*%crossprod(R))
edf2 <- rowSums(Vc*crossprod(R))
if (sum(edf2)>sum(edf1)) edf2 <- edf1
## note hat not possible here...
list(Vc=Vc,Vb=Vb,Ve=Ve,edf=edf,edf1=edf1,edf2=edf2,F=F,R=R)
......
......@@ -1195,9 +1195,8 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
# parts of the smooth terms are treated as random effects. The onesided formula random defines additional
# random terms. correlation describes the correlation structure. This routine is basically an interface
# between the basis constructors provided in mgcv and the gammPQL routine used to estimate the model.
# NOTE: need to fill out the gam object properly
{
## if (!require("nlme")) stop("gamm() requires package nlme to be installed")
{ if (inherits(family,"extended.family")) warning("family are not designed for use with gamm!")
control <- do.call("lmeControl",control)
# check that random is a named list
if (!is.null(random))
......
......@@ -75,7 +75,14 @@ jini <- function(G,lambda) {
X[(nobs+1):(nobs+jj),uoff[i]:(uoff[i]+ncol(S)-1)] <- S
nobs <- nobs + jj
}
qr.coef(qr(X),z)
## we need some idea of initial coeffs and some idea of
## associated standard error...
qrx <- qr(X,LAPACK=TRUE)
rp <- qrx$pivot;rp[rp] <- 1:ncol(X)
Ri <- backsolve(qr.R(qrx),diag(1,nrow=ncol(X)))[rp,]
beta <- qr.coef(qrx,z)
se <- sqrt(rowSums(Ri^2))*sqrt(sum((z-X%*%beta)^2)/nrow(X))
list(beta=beta,se=se)
} ## jini
jagam <- function(formula,family=gaussian,data=list(),file,weights=NULL,na.action,
......@@ -154,11 +161,22 @@ sp.prior = "gamma",diagonalize=FALSE) {
if (use.weights) jags.stuff$w <- weights
if (family$family == "binomial") jags.stuff$y <- G$y*weights ## JAGS not expecting observed prob!!
## get initial values, for use by JAGS, and to guess suitable values for
## uninformative priors...
lambda <- initial.spg(G$X,G$y,G$w,family,G$S,G$off,G$L) ## initial sp values
jags.ini <- list()
lam <- if (is.null(G$L)) lambda else G$L%*%lambda
jin <- jini(G,lam)
jags.ini$b <- jin$beta
prior.tau <- signif(0.01/(abs(jin$beta) + jin$se)^2,2)
## set the fixed effect priors...
if (G$nsdf>0) {
cat(" ## Parameteric effect priors CHECK tau is appropriate!\n",file=file,append=TRUE)
cat(" for (i in 1:",G$nsdf,") { b[i] ~ dnorm(0,0.001) }\n",file=file,append=TRUE,sep="")
ptau <- min(prior.tau[1:G$nsdf])
cat(" ## Parametric effect priors CHECK tau=1/",signif(1/sqrt(ptau),2),"^2 is appropriate!\n",file=file,append=TRUE,sep="")
cat(" for (i in 1:",G$nsdf,") { b[i] ~ dnorm(0,",ptau,") }\n",file=file,append=TRUE,sep="")
}
## Work through smooths.
......@@ -191,9 +209,10 @@ sp.prior = "gamma",diagonalize=FALSE) {
if (seperable) {
b0 <- G$smooth[[i]]$first.para
if (M==0) {
cat(" ## Note fixed vague prior, CHECK tau...\n",file=file,append=TRUE,sep="")
cat(" ## Note fixed vague prior, CHECK tau = 1/",signif(1/sqrt(ptau),2),"^2...\n",file=file,append=TRUE,sep="")
b1 <- G$smooth[[i]]$last.para
cat(" for (i in ",b0,":",b1,") { b[i] ~ dnorm(0, 1e-6) }\n",file=file,append=TRUE,sep="")
ptau <- min(prior.tau[b0:b1])
cat(" for (i in ",b0,":",b1,") { b[i] ~ dnorm(0,",ptau,") }\n",file=file,append=TRUE,sep="")
} else for (j in 1:M) {
D <- diag(G$smooth[[i]]$S[[j]]) > 0
b1 <- sum(as.numeric(D)) + b0 - 1
......@@ -222,10 +241,7 @@ sp.prior = "gamma",diagonalize=FALSE) {
} ## smoothing penalties finished
## Write the smoothing parameter prior code, using L if it exists.
lambda <- initial.spg(G$X,G$y,G$w,family,G$S,G$off,G$L) ## initial sp values
jags.ini <- list()
lam <- if (is.null(G$L)) lambda else G$L%*%lambda
jags.ini$b <- jini(G,lam)
cat(" ## smoothing parameter priors CHECK...\n",file=file,append=TRUE,sep="")
if (is.null(G$L)) {
if (sp.prior=="log.uniform") {
......
This diff is collapsed.
......@@ -52,6 +52,78 @@ mvn.ll <- function(y,X,beta,dbeta=NULL) {
list(l=oo$ll,lb=oo$lb,lbb=matrix(oo$lbb,nb,nb),dH=dH)
}
## discretized covariate routines...
XWXd <- function(X,w,k,ts,dt,v,qc,nthreads=1,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1) {
## Form X'WX given weights in w and X in compressed form in list X.
## each element of X is a (marginal) model submatrix. Full version
## is given by X[[i]][k[,i],]. list X relates to length(ds) separate
## terms. ith term starts at matrix ts[i] and has dt[i] marginal matrices.
## Terms with several marginals are tensor products and may have
## constraints (if qc[i]>1), stored as a householder vector in v[[i]].
## check ts and k index start (assumed 1 here)
## if drop is non-NULL it contains index of rows/cols to drop from result
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
n <- length(w);pt <- 0;
for (i in 1:nt) pt <- pt + prod(p[ts[i]:(ts[i]+dt[i]-1)]) - as.numeric(qc[i]>0)
oo <- .C(C_XWXd,XWX =as.double(rep(0,pt^2)),X= as.double(unlist(X)),w=as.double(w),
k=as.integer(k-1),m=as.integer(m),p=as.integer(p), n=as.integer(n),
ns=as.integer(nx), ts=as.integer(ts-1), as.integer(dt), nt=as.integer(nt),
v = as.double(unlist(v)),qc=as.integer(qc),nthreads=as.integer(nthreads),
ar.stop=as.integer(ar.stop-1),ar.row=as.integer(ar.row-1),ar.weights=as.double(ar.w))
if (is.null(drop)) matrix(oo$XWX,pt,pt) else matrix(oo$XWX,pt,pt)[-drop,-drop]
} ## XWXd
XWyd <- function(X,w,y,k,ts,dt,v,qc,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1) {
## X'Wy...
m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
n <- length(w);pt <- 0;
for (i in 1:nt) pt <- pt + prod(p[ts[i]:(ts[i]+dt[i]-1)]) - as.numeric(qc[i]>0)
oo <- .C(C_XWyd,XWy=rep(0,pt),y=as.double(y),X=as.double(unlist(X)),w=as.double(w),k=as.integer(k-1),
m=as.integer(m),p=as.integer(p),n=as.integer(n), nx=as.integer(nx), ts=as.integer(ts-1),
dt=as.integer(dt),nt=as.integer(nt),v=as.double(unlist(v)),qc=as.integer(qc),
ar.stop=as.integer(ar.stop-1),ar.row=as.integer(ar.row-1),ar.weights=as.double(ar.w))
if (is.null(drop)) oo$XWy else oo$XWy[-drop]
} ## XWyd
Xbd <- function(X,beta,k,ts,dt,v,qc,drop=NULL) {
## note that drop may contain the index of columns of X to drop before multiplying by beta.
## equivalently we can insert zero elements into beta in the appropriate places.
n <- nrow(k);m <- unlist(lapply(X,nrow));p <- unlist(lapply(X,ncol))
nx <- length(X);nt <- length(ts)
if (!is.null(drop)) {
b <- rep(0,length(beta)+length(drop))
b[-drop] <- beta
beta <- b
}
oo <- .C(C_Xbd,f=as.double(rep(0,n)),beta=as.double(beta),X=as.double(unlist(X)),k=as.integer(k-1),
m=as.integer(m),p=as.integer(p), n=as.integer(n), nx=as.integer(nx), ts=as.integer(ts-1),
as.integer(dt), as.integer(nt),as.double(unlist(v)),as.integer(qc))
oo$f
} ## Xbd
dchol <- function(dA,R) {
## if dA contains matrix dA/dx where R is chol factor s.t. R'R = A
## then this routine returns dR/dx...
p <- ncol(R)
oo <- .C(C_dchol,dA=as.double(dA),R=as.double(R),dR=as.double(R*0),p=as.integer(ncol(R)))
return(matrix(oo$dR,p,p))
} ## dchol
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.
p <- ncol(dR[[1]])
M <- if (trans) ncol(Vr) else -ncol(Vr) ## sign signals transpose or not to C code
if (abs(M)!=length(dR)) stop("internal error in vcorr, please report to simon.wood@r-project.org")
oo <- .C(C_vcorr,dR=as.double(unlist(dR)),Vr=as.double(Vr),Vb=as.double(rep(0,p*p)),
p=as.integer(p),M=as.integer(M))
return(matrix(oo$Vb,p,p))
} ## vcorr
pinv <- function(X,svd=FALSE) {
## a pseudoinverse for n by p, n>p matrices
qrx <- qr(X,tol=0,LAPACK=TRUE)
......@@ -93,7 +165,7 @@ pqr2 <- function(x,nt=1,nb=30) {
pbsi <- function(R,nt=1,copy=TRUE) {
## parallel back substitution inversion of upper triangular R
## library(mgcv); n <- 5000;p<-4000;x <- matrix(runif(n*p),n,p)
## library(mgcv); n <- 500;p<-400;x <- matrix(runif(n*p),n,p)
## qrx <- mgcv:::pqr2(x,2);R <- qr.R(qrx)
## system.time(Ri <- mgcv:::pbsi(R,2))
## system.time(Ri2 <- backsolve(R,diag(p)));range(Ri-Ri2)
......@@ -117,6 +189,20 @@ pchol <- function(A,nt=1,nb=30) {
A
}
pforwardsolve <- function(R,B,nt=1) {
## parallel forward solve via simple col splitting...
if (!is.matrix(B)) B <- as.matrix(B)
.Call(C_mgcv_Rpforwardsolve,R,B,nt)
}
pcrossprod <- function(A,trans=FALSE,nt=1,nb=30) {
## parallel cross prod A'A or AA' if trans==TRUE...
if (!is.matrix(A)) A <- as.matrix(A)
if (trans) A <- t(A)
.Call(C_mgcv_Rpcross,A,nt,nb)
}
pRRt <- function(R,nt=1) {
## parallel RR' for upper triangular R
## following creates index of lower triangular elements...
......
......@@ -508,7 +508,7 @@ tensor.prod.model.matrix1 <- function(X) {
# X is a list of model matrices, from which a tensor product model matrix is to be produced.
# e.g. ith row is basically X[[1]][i,]%x%X[[2]][i,]%x%X[[3]][i,], but this routine works
# column-wise, for efficiency
# old version, which is rather slow becuase of using cbind.
# old version, which is rather slow because of using cbind.
m <- length(X)
X1 <- X[[m]]
n <- nrow(X1)
......@@ -603,8 +603,8 @@ smooth.construct.tensor.smooth.spec <- function(object,data,knots)
}
XP <- list()
if (object$np) # reparameterize
for (i in 1:m)
{ if (object$margin[[i]]$dim==1) {
for (i in 1:m) {
if (object$margin[[i]]$dim==1) {
# only do classes not already optimal (or otherwise excluded)
if (!inherits(object$margin[[i]],c("cs.smooth","cr.smooth","cyclic.smooth","random.effect"))) {
x <- get.var(object$margin[[i]]$term,data)
......@@ -625,7 +625,7 @@ smooth.construct.tensor.smooth.spec <- function(object,data,knots)
warning("reparameterization unstable for margin: not done")