Commit 2ffccbf0 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-5

parent 3c9df280
** denotes quite substantial/important changes
*** denotes really big changes
1.8-5
* scale parameter handling in edf in logLik.gam made consistent with glm
(affects AIC).
* 'bam', 'gam' and 'gamm' modified to often produce smaller files when models
saved (and never to produce absurdly large files). Achieved by setting
environment of formula, terms etc to .GlobalEnv. Previously 'save' could
save entire contents of environment of formula/terms with fitted model
object. Note that change can cause failure in user written functions calling
gam/bam and then 'predict' without supplying all prediction variables
(fix obvious).
* A help file 'single.index' supplied illustrating how single index models
can be estimated in mgcv.
* predict.gam now only creates a "constant" attribute if the model has one.
* gam.fit4 convergence testing of coefs modified to more robust test of
gradients of penalized dev w.r.t. params, rather than change in params,
which can fail under rank deficiency.
* mgcv_qrqy was not thread safe. Not noticeable on many platforms as all
threads did exactly the same thing to the same matrix, but very noticeable
on Windows. Thread safe mgcv_qrqy0 added and used in any parallel sections.
* Allow openMP support if compiler supports it and provides pre-defined macro
_OPENMP, even if SUPPORT_OPENMP undefined. (Allows multi-threading on
Windows, for example.)
* 'eps' is now an argument to 'betar' allowing some control on how to
handle response values too close to 0 or 1. Help file expanded to
emphasise the problems with using beta regression with 0s and 1s in
the data.
* fix of bug in multi-formula contrast handling, causing failure of prediction
in some cases.
* ziP and ziplss now check for non-integer (or binary) responses and produce
an error message if these are found. Previously this was not trapped and
could lead to a segfault.
1.8-4
** JAGS/BUGS support added, enabling auto-generation of code and data
......
Package: mgcv
Version: 1.8-4
Version: 1.8-5
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
estimation
Description: Routines for GAMs and other generalized ridge regression
with multiple smoothing parameter selection by GCV, REML
or UBRE/AIC. Also GAMMs. Includes a gam() function.
Title: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness
Estimation
Description: GAMs, GAMMs and other generalized ridge regression with
multiple smoothing parameter estimation by GCV, REML or UBRE/AIC.
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
......@@ -14,7 +15,7 @@ Suggests: splines, parallel, survival, MASS
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2014-11-27 09:49:02 UTC; sw283
Packaged: 2015-03-02 21:05:39 UTC; sw283
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2014-11-27 13:16:04
Date/Publication: 2015-03-03 09:18:09
0ea61afceef0bf25b8c9d74103884936 *ChangeLog
1a88c7e031e86b1a4d2f0a3733714046 *DESCRIPTION
dcdf699f4b790dfe04b100f3c980436c *ChangeLog
cadc1146dc050c7dc00162207980193b *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
68be7e47206e4b5ee210fb6de27f8169 *NAMESPACE
27f3c91b6ebcefafbf8f370ea38d08c9 *R/bam.r
af08e19c4214ac86d453bd90e1ce9e4e *NAMESPACE
4d8460a0933930189b287f8f77bdfabf *R/bam.r
52ea7081e235df890d3e51198bcc3092 *R/coxph.r
d6f40f3d2acd3e30f3a59cbc85da6bc3 *R/efam.r
a824a6e49020786a1eba60b99939eff1 *R/fast-REML.r
065215b64c502afab215f6e97e9ec915 *R/gam.fit3.r
6586de2772276237525c646362926f8b *R/gam.fit4.r
8e0796908b72ff7d4d5f57740050a488 *R/efam.r
8cde8b6379034d9eb78caaf1b0d04cf8 *R/fast-REML.r
5df55d3cb82ac4669fc89909769d5a95 *R/gam.fit3.r
cbeb60506c33d463465efb5d1ada43f5 *R/gam.fit4.r
e63276b78a4ff2566af2e651bfc6aa9c *R/gam.sim.r
bd651b8eec80d83d6e0de5f47715a0de *R/gamlss.r
449b21801a62d0d6b753e753d4861570 *R/gamm.r
739156695988beabea60d2e259073821 *R/jagam.r
b5887ec9a16b1842918901bc5e4ecc94 *R/mgcv.r
a0d1da72334aefd546cad58d1b2062b1 *R/misc.r
c287514d201a02f060fee98aab1e93c8 *R/mvam.r
8e3b23448498fd03e668165ca3b01a60 *R/plots.r
92302835adba3710624e434ad52173bf *R/smooth.r
90f72ef91dfc2a25c262fc002f781bf7 *R/soap.r
0c5d2acdc985d2af9de6bb736d4df2ab *R/sparse.r
9a4ef9d8c7362256e314be0e65b5a96c *R/gamlss.r
28ddc8624513070a62d0cc2aba7776a2 *R/gamm.r
09240f99d77e54848bf15d540c267022 *R/jagam.r
9a451b851f1e08fb7c7c60738556392d *R/mgcv.r
2f7ab8cc9f9cf2f6cb52fc72f9f69bfb *R/misc.r
6c53ea135a7ac956fbe8a3d16c59baaa *R/mvam.r
58413439e82b60d8942d070ca72f8e14 *R/plots.r
3e87a7ef880c8aa40ce57201f6ea5458 *R/smooth.r
68348617c5d3b0e2ea805c4360a8cdd4 *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
074adce4131eced4cc71c67ad1d63c52 *inst/CITATION
20b3c68124572bfdbba2938b6cda8ada *man/Beta.Rd
5280bb00fa6595ff4a4b70d2e80114a8 *inst/po/de/LC_MESSAGES/R-mgcv.mo
c745969a1292eb3d49dfd0d0c2c997d4 *inst/po/de/LC_MESSAGES/mgcv.mo
9a49cc8d28e67874337cbe11d8632138 *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
447116b9728f9a82da5ce5d821850917 *inst/po/en@quot/LC_MESSAGES/mgcv.mo
2d98abb49e987bb8ebdf9be5cdbbb68b *inst/po/fr/LC_MESSAGES/R-mgcv.mo
643671acdb430cb0790f43559addbb0d *inst/po/fr/LC_MESSAGES/mgcv.mo
ea6f2b5c19372fda131171a03c983d6c *inst/po/ko/LC_MESSAGES/R-mgcv.mo
c1b1475e5fef49fe49929d2796ff87b6 *inst/po/ko/LC_MESSAGES/mgcv.mo
f8a930e0893c980fc2e9cfe276107e0f *inst/po/pl/LC_MESSAGES/R-mgcv.mo
715e52c0debf9848bbda15e94f5e7315 *inst/po/pl/LC_MESSAGES/mgcv.mo
251140c0e4b621fd578dae4ce8a28c8f *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
e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
......@@ -45,7 +57,7 @@ c7561a0f2e114247955e212aeabc6ae9 *man/formXtViX.Rd
d87ff6287d7e343ea24d2053f4724b35 *man/formula.gam.Rd
4da4d585b329769eb44f0c7a6e7dd554 *man/fs.test.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
a1492616b069844ea05762687473c4fd *man/gam.Rd
a54ea44bb9641389bcbd21ce0e578d8f *man/gam.Rd
fe61dd0efab9e920c17335faf3d5764c *man/gam.check.Rd
a65bc22f606e45d185bc375fbf5698a1 *man/gam.control.Rd
44db24b66ce63bc16d2c8bc3f5b42ac5 *man/gam.convergence.Rd
......@@ -68,7 +80,7 @@ a2ea1233d43fac89e0cacbc09a8d31e2 *man/in.out.Rd
2f222eeeb3d7bc42f93869bf8c2af58a *man/influence.gam.Rd
bafea2eef12fdc819f8ac1fb41d8b914 *man/initial.sp.Rd
c00bcfe2d0b44b2ea955f3934421807c *man/interpret.gam.Rd
ae18c335a5bb7192ad608a5a500c2794 *man/jagam.Rd
469e313e7a1b0520593c8e7a938111b5 *man/jagam.Rd
07d2c259b9edf164f42935170b4fccd0 *man/ldTweedie.Rd
58e73ac26b93dc9d28bb27c8699e12cf *man/linear.functional.terms.Rd
93035193b0faa32700e1421ce8c1e9f6 *man/logLik.gam.Rd
......@@ -77,7 +89,7 @@ ae18c335a5bb7192ad608a5a500c2794 *man/jagam.Rd
5169af4be5fccf9fa79b6de08e9ea035 *man/magic.post.proc.Rd
59672d8599211ce9cf44252f92110920 *man/mgcv-FAQ.Rd
01926b8c2bda6fc74b51f6ec7f5cc620 *man/mgcv-package.Rd
934b90e0845a4971b97f1d67d12d61ec *man/mgcv-parallel.Rd
92499f4c6bc65a0cffd31bd8bc993a9f *man/mgcv-parallel.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
bc9b89db7e7ff246749551c16f5f1f07 *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
......@@ -93,9 +105,9 @@ c64b44bc684cbadfe77e7f6dc833ddb4 *man/pcls.Rd
8c0f8575b427f30316b639a326193aeb *man/pdTens.Rd
1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd
edf57071572275a8443b2f0b66d44424 *man/place.knots.Rd
4aed50caabd7f058f90437dab4cdbee0 *man/plot.gam.Rd
c488918033a996ce97943bf20ff8ac54 *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
afca36f5b1a5d06a7fcab2eaaa029e7e *man/predict.bam.Rd
085dff270e5255315cb244368faf79df *man/predict.bam.Rd
900647986ab1365d437368cc28b62658 *man/predict.gam.Rd
7562cb5ce9a7fbf0cd60124cf4305315 *man/print.gam.Rd
b8ac3ed8fc05bb14c605b996404853cd *man/qq.gam.Rd
......@@ -105,6 +117,7 @@ c523210ae95cb9aaa0aaa1c37da1a4c5 *man/residuals.gam.Rd
f6f1333be2587ffef5970905b13468ea *man/rig.Rd
7258dfc0149fff020054117fd2ee6bd8 *man/s.Rd
c438eb3e41cb1f51e72283380145d210 *man/scat.Rd
2199afd400a1821b74a24e2578bc0224 *man/single.index.Rd
6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
b5a06918956fd10f23285b453c05bdb4 *man/smooth.construct.Rd
4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.smooth.spec.Rd
......@@ -130,37 +143,39 @@ f0791d830687d6155efb8a73db787401 *man/summary.gam.Rd
6eebb6ef90374ee09453d6da6449ed79 *man/tensor.prod.model.matrix.Rd
71add05e5823dcc0509f435ab4dc4627 *man/uniquecombs.Rd
a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
0e045fc94ea00c17e703ae07923f549b *man/vis.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
2c5f6815e609f2cdb56b0067a183f915 *man/ziP.Rd
455f8e0a6cd4218442cfdba368cbb940 *man/ziplss.Rd
becbe3e1f1588f7292a74a97ef07a9ae *po/R-de.po
ae8388103d8b1be39f55f426b205b576 *man/ziplss.Rd
a7c1582252be4e0013e4275ffd3aac4c *po/R-de.po
0bdfcf98961b0d52b60f806dc1dba77e *po/R-en@quot.po
9126ed91030cd1c168d669854cf1f510 *po/R-fr.po
7563cdf47066589a08ba01eacfcf0bcb *po/R-mgcv.pot
ccbf345b8ca9544e9239fa85fff897a1 *po/R-pl.po
d5a0f198090ecbfedaa6549f2918b997 *po/R-po.po
3550a794504bdfd4f75a83de1148bb40 *po/de.po
eb3707a62b0e3d406ddf5c15bc40c673 *po/R-fr.po
c7317bc4fb746b3ce2a9bcdf9d928716 *po/R-ko.po
5e7691372a93b11f3c0efa48ec9dae2f *po/R-mgcv.pot
62fa356ec68d7cc8a9ec1b3915b0db13 *po/R-pl.po
a5335d43184be0f9da9df45cfadc3e33 *po/R-po.po
ccf3140169b68ec7aff4b15e6a97e5db *po/de.po
93f72334356fe6f05a64e567efd35c8e *po/en@quot.po
1a4a267ddcb87bb83f09c291d3e97523 *po/fr.po
813514ea4e046ecb4563eb3ae8aa202a *po/mgcv.pot
bd5ec325242ca6e41056c2f1f7f05bfe *po/pl.po
749ac663240fd6eaa2b725602d47ef2a *po/po.po
fb829b82760779929951d49fe29ed2e5 *po/fr.po
dc1ef92ff4454734c3a24876e299b760 *po/ko.po
5fee1431cf41b342cde9e5eb6dd27cab *po/mgcv.pot
dfd4eec9edc7d1ab6354d47b6e2bd42f *po/pl.po
b56dac4547037ea45e2c8f9bce7aa9ef *po/po.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
56d501c682dc6699ff0803f25533e849 *src/coxph.c
114701c0d407f946ac693da55da864da *src/gdi.c
49af97195accb65adc75620183d39a4c *src/general.h
761581e96a73e1bec00a9553356eaaa6 *src/init.c
6dbd109048a3bd18b67db6230c064c21 *src/magic.c
f45aa90d4ce295f0deebc54798dad330 *src/mat.c
02122df2193f86e2a08f68fe5b8ff972 *src/matrix.c
e71c1a1624b431fbab0a4c8f151d2a97 *src/coxph.c
08e156711c686a1f1efe505d63fabef5 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
e5bf24371be5ea7f3ffb40060a648803 *src/init.c
9a5d7cb3cf93cbdfc08353fbd20a270e *src/magic.c
a9735df73ae117df1b4c1197f4748389 *src/mat.c
aae0f298e384952bb8b6b923d40520a8 *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
7a12ef0cbd8f5312cd67669a36d75ebf *src/mgcv.c
f06b1075027c76e9a7ce483c3d70c6ad *src/mgcv.h
ec2ae157a9e3bedcccc88b053d0a4e0c *src/mgcv.c
b62374abb9910e539a3280d66ac4193f *src/mgcv.h
00f8a024faef17f90ed04f01e736df71 *src/misc.c
08c1706ffeec4277c484435a0644b0e3 *src/mvn.c
90fe2ad1997db56f2293ad16c7110f88 *src/qp.c
cbe54250deb38aa5f88f8b669a4468cd *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h
15ed450da804a9dbd0833459af1f2474 *src/soap.c
b3ff61207cd26308abf386d6a0666320 *src/sparse-smooth.c
75d9bf91f5c2dc536918055b6a00d93f *src/tprs.c
294190980f5a98be73fa9a657f4d0b91 *src/tprs.c
5bd85bf0319a7b7c755cf49c91a7cd94 *src/tprs.h
......@@ -125,6 +125,8 @@ S3method(smooth.construct, mrf.smooth.spec)
S3method(smooth.construct,ps.smooth.spec)
S3method(smooth.construct, re.smooth.spec)
S3method(smooth.construct,so.smooth.spec)
S3method(smooth.construct,sw.smooth.spec)
S3method(smooth.construct,sf.smooth.spec)
S3method(smooth.construct,sos.smooth.spec)
S3method(smooth.construct, tp.smooth.spec)
S3method(smooth.construct, tensor.smooth.spec)
......
......@@ -185,13 +185,13 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
weights <- G$w
conv <- FALSE
nobs <- nrow(mf)
nvars <- ncol(G$X)
##nvars <- ncol(G$X)
offset <- G$offset
family <- G$family
G$family <- gaussian() ## needed if REML/ML used
variance <- family$variance
dev.resids <- family$dev.resids
aic <- family$aic
## aic <- family$aic
linkinv <- family$linkinv
mu.eta <- family$mu.eta
if (!is.function(variance) || !is.function(linkinv))
......@@ -211,7 +211,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
mustart <- mukeep
}
coefold <- NULL
##coefold <- NULL
eta <- if (!is.null(etastart))
etastart
else family$linkfun(mustart)
......@@ -220,7 +220,8 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
if (!(validmu(mu) && valideta(eta)))
stop("cannot find valid starting values: please specify some")
dev <- sum(dev.resids(y, mu, weights))*2 ## just to avoid converging at iter 1
boundary <- conv <- FALSE
##boundary <-
conv <- FALSE
G$coefficients <- rep(0,ncol(G$X))
class(G) <- "gam"
......@@ -506,13 +507,13 @@ bgam.fit2 <- function (G, mf, chunk.size, gp ,scale ,gamma,method, etastart = NU
weights <- G$w
conv <- FALSE
nobs <- nrow(mf)
nvars <- ncol(G$X)
##nvars <- ncol(G$X)
offset <- G$offset
family <- G$family
G$family <- gaussian() ## needed if REML/ML used
variance <- family$variance
dev.resids <- family$dev.resids
aic <- family$aic
##aic <- family$aic
linkinv <- family$linkinv
mu.eta <- family$mu.eta
if (!is.function(variance) || !is.function(linkinv))
......@@ -532,7 +533,7 @@ bgam.fit2 <- function (G, mf, chunk.size, gp ,scale ,gamma,method, etastart = NU
mustart <- mukeep
}
coefold <- NULL
##coefold <- NULL
eta <- if (!is.null(etastart))
etastart
else family$linkfun(mustart)
......@@ -541,7 +542,8 @@ bgam.fit2 <- function (G, mf, chunk.size, gp ,scale ,gamma,method, etastart = NU
if (!(validmu(mu) && valideta(eta)))
stop("cannot find valid starting values: please specify some")
dev <- sum(dev.resids(y, mu, weights))*2 ## just to avoid converging at iter 1
boundary <- conv <- FALSE
##boundary <-
conv <- FALSE
G$n <- nobs
X <- G$X
......@@ -733,7 +735,7 @@ predict.bam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
gc()
if (!is.null(cluster)&&inherits(cluster,"cluster")) {
require(parallel)
## require(parallel)
n.threads <- length(cluster)
} else n.threads <- 1
if (missing(newdata)) n <- nrow(object$model) else n <- nrow(newdata)
......@@ -814,7 +816,7 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
if (n>chunk.size) { ## then use QR accumulation approach
if (!is.null(cl)&&inherits(cl,"cluster")) {
require(parallel)
## require(parallel)
n.threads <- length(cl)
} else n.threads <- 1
......@@ -926,7 +928,7 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
# }
## now consolidate the results from the parallel threads...
R <- res[[1]]$R;f <- res[[1]]$f;dev <- res[[1]]$dev
R <- res[[1]]$R;f <- res[[1]]$f; ## dev <- res[[1]]$dev
y.norm2 <- res[[1]]$y.norm2
for (i in 2:n.threads) {
if (use.chol) {
......@@ -1195,15 +1197,15 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
H=NULL,absorb.cons=TRUE,sparse.cons=sparse.cons,select=FALSE,
idLinksBases=TRUE,scale.penalty=control$scalePenalty,
paraPen=paraPen)
## no advantage to "fREML" with no free smooths...
if (((!is.null(G$L)&&ncol(G$L) < 1)||(length(G$sp)==0))&&method=="fREML") method <- "REML"
G$var.summary <- var.summary
G$family <- family
G$terms<-terms;##G$pterms<-pterms
pvars <- all.vars(delete.response(terms))
G$pred.formula <- gp$pred.formula # if (length(pvars)>0) reformulate(pvars) else NULL
G$terms<-terms;
G$pred.formula <- gp$pred.formula
n <- nrow(mf)
......@@ -1227,8 +1229,11 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
if (G$m) for (i in 1:G$m) G$min.edf<-G$min.edf+G$smooth[[i]]$null.space.dim
G$formula<-formula
environment(G$formula)<-environment(formula)
## environment(G$formula)<-environment(formula)
environment(G$pterms) <- environment(G$terms) <- environment(G$pred.formula) <-
environment(G$formula) <- .BaseNamespaceEnv
G$conv.tol<-control$mgcv.tol # tolerence for mgcv
G$max.half<-control$mgcv.half # max step halving in bfgs optimization
......@@ -1297,7 +1302,6 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
object$family <- family
object$formula<-G$formula
#object$linear.predictors <- NA
if (method=="GCV.Cp") {
if (scale<=0) object$method <- "GCV" else object$method <- "UBRE"
} else {
......@@ -1345,6 +1349,9 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
if (length(object$full.sp)==length(object$sp)&&
all.equal(object$sp,object$full.sp)==TRUE) object$full.sp <- NULL
}
environment(object$formula) <- environment(object$pred.formula) <-
environment(object$terms) <- environment(object$pterms) <-
environment(attr(object$model,"terms")) <- .GlobalEnv
object
} ## end of bam
......@@ -1451,7 +1458,6 @@ bam.update <- function(b,data,chunk.size=10000) {
} else if (method=="fREML") { ## fast REML
um <- Sl.Xprep(b$Sl,b$qrx$R)
lambda.0 <- initial.sp(b$qrx$R,b$G$S,b$G$off)
lsp0 <- log(b$sp) ## initial s.p.
log.phi <- log(b$sig2) ## initial or fixed scale
fit <- fast.REML.fit(um$Sl,um$X,b$qrx$f,rho=lsp0,L=b$G$L,rho.0=b$G$lsp0,
......
This diff is collapsed.
......@@ -85,7 +85,7 @@ Sl.setup <- function(G) {
## overlap testing requires the block ranges
for (j in 1:m) { ## get block range for each S[[j]]
ir <- range((1:nb)[rowSums(abs(Sl[[b]]$S[[j]]))>0])
ic <- range((1:nb)[colSums(abs(Sl[[b]]$S[[j]]))>0])
#ic <- range((1:nb)[colSums(abs(Sl[[b]]$S[[j]]))>0]) ## symmetric not needed
sbStart[j] <- ir[1];sbStop[j] <- ir[2] ## individual ranges
}
split.ok <- TRUE
......
......@@ -306,7 +306,8 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
}
else family$linkfun(mustart)
#etaold <- eta
muold <- mu <- linkinv(eta)
##muold <-
mu <- linkinv(eta)
#if (!(validmu(mu) && valideta(eta)))
# stop("Can't find valid starting values: please specify some")
......@@ -490,7 +491,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
old.pdev <- pdev
coef <- coefold <- start
etaold <- eta
muold <- mu
##muold <- mu
} else {
conv <- TRUE
coef <- start
......@@ -1222,7 +1223,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
weights = weights, start = start,
offset = offset,Mp=Mp,family = family,
control = control,deriv=deriv,eps=eps,spe=spe,
Sl=Sl,...)
Sl=Sl,...) ## ignore codetools warning
}
} ## end of derivative checking
......@@ -1419,7 +1420,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
if (kk==1||score3<=score2) { ## accept step - better than last try
score2 <- score3
lsp2 <- lsp3
if (!is.null(lsp.max)) delta2 <- delta3
## if (!is.null(lsp.max)) delta2 <- delta3
}
## stop when improvement found, and shorter step is worse...
if (score2<score&&score3>score2||kk==40) ok <- FALSE
......@@ -1854,7 +1855,7 @@ gam2derivative <- function(lsp,args,...)
b<-gam.fit3(x=args$X, y=args$y, sp=lsp,Eb=args$Eb,UrS=args$UrS,
offset = args$offset,U1=args$U1,Mp=args$Mp,family = args$family,weights=args$w,deriv=1,
control=args$control,gamma=args$gamma,scale=args$scale,scoreType=args$scoreType,
null.coef=args$null.coef,...)
null.coef=args$null.coef,n.true=args$n.true,...)
if (reml) {
ret <- b$REML1
} else if (args$scoreType=="GACV") {
......@@ -1878,7 +1879,7 @@ gam2objective <- function(lsp,args,...)
b<-gam.fit3(x=args$X, y=args$y, sp=lsp,Eb=args$Eb,UrS=args$UrS,
offset = args$offset,U1=args$U1,Mp=args$Mp,family = args$family,weights=args$w,deriv=0,
control=args$control,gamma=args$gamma,scale=args$scale,scoreType=args$scoreType,
null.coef=args$null.coef,...)
null.coef=args$null.coef,n.true=args$n.true,...)
if (reml) {
ret <- b$REML
} else if (args$scoreType=="GACV") {
......
......@@ -210,7 +210,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
sp <- sp[-ind] ## log smoothing parameters
} else theta <- family$getTheta() ## fixed value
penalized <- if (length(UrS)>0) TRUE else FALSE
## penalized <- if (length(UrS)>0) TRUE else FALSE
if (scale>0) scale.known <- TRUE else {
## unknown scale parameter, trial value supplied as
......@@ -292,8 +292,8 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
else family$linkfun(mustart)
mu.eta <- family$mu.eta
Dd <- family$Dd
##mu.eta <- family$mu.eta
##Dd <- family$Dd
linkinv <- family$linkinv
valideta <- family$valideta
......@@ -445,11 +445,16 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
## convergence testing...
if (abs(pdev - old.pdev)/(0.1 + abs(pdev)) < control$epsilon) {
if (max(abs(start-coefold))>control$epsilon*max(abs(start+coefold))/2) {
## Need to check coefs converged adequately, to ensure implicit differentiation
## ok. Testing coefs unchanged is problematic under rank deficiency (not guaranteed to
## drop same parameter every iteration!)
grad <- 2 * t(x[good,])%*%(w*((x%*%start)[good]-z))+ 2*St%*%start
if (max(abs(grad)) > control$epsilon*max(abs(start+coefold))/2) {
## if (max(abs(start-coefold))>control$epsilon*max(abs(start+coefold))/2) {
old.pdev <- pdev ## not converged quite enough
coef <- coefold <- start
etaold <- eta
muold <- mu
##muold <- mu
} else { ## converged
conv <- TRUE
coef <- start
......@@ -491,7 +496,7 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
residuals[good] <- z - (eta - offset)[good]
ntot <- length(theta) + length(sp)
if (deriv>1) n2d <- ntot*(1+ntot)/2 else n2d <- 0
## if (deriv>1) n2d <- ntot*(1+ntot)/2 else n2d <- 0
rSncol <- unlist(lapply(UrS,ncol))
## Now drop any elements of dd that have been dropped in fitting...
if (sum(!good)>0) { ## drop !good from fields of dd
......@@ -651,10 +656,11 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
penalized <- if (length(Sl)>0) TRUE else FALSE
nSp <- length(lsp)
sp <- exp(lsp)
rank.tol <- .Machine$double.eps*100 ## tolerance to use for rank deficiency
##sp <- exp(lsp)
## rank.tol <- .Machine$double.eps*100 ## tolerance to use for rank deficiency
q <- ncol(x)
n <- nobs <- length(y)
##n <-
nobs <- length(y)
if (penalized) {
Eb <- attr(Sl,"E") ## balanced penalty sqrt
......@@ -967,7 +973,7 @@ 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>1) -(d2l - d2bSb/2 + rp$ldet2/2 - d2ldetH/2) else NULL
bSb <- t(coef)%*%St%*%coef
## bSb <- t(coef)%*%St%*%coef
lpi <- attr(x,"lpi")
if (is.null(lpi)) {
linear.predictors <- as.numeric(x%*%coef)
......
......@@ -79,7 +79,7 @@ gamlss.etamu <- function(l1,l2,l3=NULL,l4=NULL,ig1,g2,g3=NULL,g4=NULL,i2,i3=NULL
d1[,i] <- l1[,i]*ig1[,i]
}
n <- length(ig1[,1])
##n <- length(ig1[,1])
k <- 0
d2 <- l2
......@@ -659,7 +659,8 @@ zipll <- function(y,g,eta,deriv=0) {
l1 <- El2 <- l2 <- l3 <- l4 <- NULL
zind <- y == 0 ## the index of the zeroes
yz <- y[zind];yp <- y[!zind]
## yz <- y[zind];
yp <- y[!zind]
l <- et <- exp(eta)
l[zind] <- -et[zind] # -exp(eta[ind])
l[!zind] <- l1ee(eta[!zind]) + yp*g[!zind] - lee1(g[!zind]) - lgamma(yp+1)
......@@ -873,8 +874,8 @@ ziplss <- function(link=list("identity","identity")) {
eta1 <- X[,jj[[2]],drop=FALSE]%*%coef[jj[[2]]]
p <- family$linfo[[2]]$linkinv(eta1)
n <- length(y)
l1 <- matrix(0,n,2)
##n <- length(y)
## l1 <- matrix(0,n,2)
zl <- zipll(y,lambda,p,deriv)
if (deriv>0) {
......@@ -883,7 +884,8 @@ ziplss <- function(link=list("identity","identity")) {
g2 <- cbind(family$linfo[[1]]$d2link(lambda),family$linfo[[2]]$d2link(p))
}
l3 <- l4 <- g3 <- g4 <- 0 ## defaults
## l3 <- l4 <-
g3 <- g4 <- 0 ## defaults
if (deriv>1) {
## the third derivatives
......@@ -921,6 +923,10 @@ ziplss <- function(link=list("identity","identity")) {
## which is basically penalizing something different here.
## best we can do here is to use E only as a regularizer.
n <- rep(1, nobs)
if (all.equal(y,round(y))!=TRUE) {
stop("Non-integer response variables are not allowed with ziplss ")
}
if ((min(y)==0&&max(y)==1)) stop("Using ziplss for binary data makes no sense")
## should E be used unscaled or not?..
use.unscaled <- if (!is.null(attr(E,"use.unscaled"))) TRUE else FALSE
if (is.null(start)) {
......
......@@ -375,8 +375,9 @@ smooth2random.fs.interaction <- function(object,vnames,type=1) {
rind <- c(rind,k:(k+ncol(X)-1))
rinc <- c(rinc,rep(ncol(X),ncol(X)))
k <- k + n.lev * ncol(X)
if (type==1) { ## gamm form for use with lme
form <- as.formula(paste("~",term.name,"-1",sep=""))
if (type==1) { ## gamm form for use with lme
## env set to avoid 'save' saving whole environment to file...
form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
random[[i]] <- pdIdnot(form)
names(random)[i] <- object$fterm ## supplied factor name
attr(random[[i]],"group") <- object$fac ## factor supplied as part of term
......@@ -414,7 +415,7 @@ smooth2random.fs.interaction <- function(object,vnames,type=1) {
Xf <- matrix(0,nrow(object$X),0)
list(rand=random,trans.D=D,Xf=Xf,fixed=FALSE,rind=rind,rinc=rinc,
pen.ind=pen.ind) ## pen.ind==i is TRUE for coefs penalized by ith penalty
}
} ## smooth2random.fs.interaction
smooth2random.t2.smooth <- function(object,vnames,type=1) {
## takes a smooth object and turns it into a form suitable for estimation as a random effect
......@@ -454,7 +455,8 @@ smooth2random.t2.smooth <- function(object,vnames,type=1) {
group.name <- new.name("g",vnames)
vnames <- c(vnames,term.name,group.name)
if (type==1) { ## gamm form for lme
form <- as.formula(paste("~",term.name,"-1",sep=""))
## env set to avoid 'save' saving whole environment to file...
form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
random[[i]] <- pdIdnot(form)
names(random)[i] <- group.name
attr(random[[i]],"group") <- factor(rep(1,nrow(X)))
......@@ -472,7 +474,7 @@ smooth2random.t2.smooth <- function(object,vnames,type=1) {
} else Xf <- matrix(0,nrow(object$X),0)
list(rand=random,trans.D=diagU,Xf=Xf,fixed=FALSE,
rind=1:n.para,rinc=rep(n.para,n.para),pen.ind=pen.ind)
}
} ## smooth2random.t2.smooth
smooth2random.mgcv.smooth <- function(object,vnames,type=1) {
## takes a smooth object and turns it into a form suitable for estimation as a random effect
......@@ -511,7 +513,8 @@ smooth2random.mgcv.smooth <- function(object,vnames,type=1) {
term.name <- new.name("Xr",vnames)
if (type==1) { ## gamm form for lme
form <- as.formula(paste("~",term.name,"-1",sep=""))
## env set to avoid 'save' saving whole environment with fitted model object
form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
random <- list(pdIdnot(form))
group.name <- new.name("g",vnames)
names(random) <- group.name
......@@ -593,7 +596,7 @@ smooth2random.tensor.smooth <- function(object,vnames,type=1) {
}
term.name <- new.name("Xr",vnames)
form <- as.formula(paste("~",term.name,"-1",sep=""))
form <- as.formula(paste("~",term.name,"-1",sep=""),env=.GlobalEnv)
attr(form,"S") <- object$S ## this class needs penalty matrices to be supplied
random <- list(pdTens(form)) ## lme random effect class
......@@ -678,17 +681,16 @@ gamm.setup <- function(formula,pterms,
sm <- G$smooth[[pord[i]]]
sm$X <- G$X[,sm$first.para:sm$last.para,drop=FALSE]
rasm <- smooth2random(sm,names(data)) ## convert smooth to random effect and fixed effects
sm$fixed <- rasm$fixed
if (!is.null(sm$fac)) {
flev <- levels(sm$fac) ## grouping factor for smooth
n.lev <- length(flev)
} else n.lev <- 1
##n.lev <- length(flev)
} ##else n.lev <- 1
## now append constructed variables to model frame and random effects to main list
n.para <- 0 ## count random coefficients
rinc <- rind <- rep(0,0)
## rinc <- rind <- rep(0,0)
if (!sm$fixed) {
# kk <- 1;
for (k in 1:length(rasm$rand)) {
......@@ -1258,7 +1260,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
else mf <- gmf
rm(gmf)
if (nrow(mf)<2) stop("Not enough (non-NA) data to do anything meaningful")
Terms <- attr(mf,"terms")
##Terms <- attr(mf,"terms")
## summarize the *raw* input variables
## note can't use get_all_vars here -- buggy with matrices
......@@ -1570,8 +1572,18 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
object$weights <- object$prior.weights
if (!is.null(G$Xcentre)) object$Xcentre <- G$Xcentre ## column centering values
## set environments to global to avoid enormous saved object files
environment(attr(object$model,"terms")) <-
environment(object$terms) <- environment(object$pterms) <-
environment(object$formula) <-environment(object$pred.formula) <- .GlobalEnv
ret$gam <- object
environment(attr(ret$lme$data,"terms")) <- environment(ret$lme$terms) <- .GlobalEnv
if (!is.null(ret$lme$modelStruct$varStruct)) {
environment(attr(ret$lme$modelStruct$varStruct,"formula")) <- .GlobalEnv
}
if (!is.null(ret$lme$modelStruct$corStruct)) {
environment(attr(ret$lme$modelStruct$corStruct,"formula")) <- .GlobalEnv
}
class(ret) <- c("gamm","list")
ret
......
......@@ -26,7 +26,7 @@ write.jagslp <- function(resp,family,file,use.weights,offset=FALSE) {
cat(" for (i in 1:n) { mu[i] <- ",iltab[family$link],"} ## expected response\n",file=file,append=TRUE)
}
## code the response given mu and any scale parameter prior...
scale <- TRUE ## is scale parameter free?
#scale <- TRUE ## is scale parameter free?
cat(" for (i in 1:n) { ",file=file,append=TRUE)
if (family$family=="gaussian") {
if (use.weights) cat(resp,"[i] ~ dnorm(mu[i],tau*w[i]) } ## response \n",sep="",file=