Commit dafd9351 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Debian changes 1.8-17-1

mgcv (1.8-17-1) unstable; urgency=medium

  * New upstream release
parents f9e25031 a9b37307
......@@ -5,6 +5,34 @@ Currently deprecated and liable to be removed:
- single penalty tensor product smooths.
- p.type!=0 in summary.gam.
1.8-17
* Export gamlss.etamu, gamlss.gH and trind.generator to facilitate user
addition of new location-scale families.
* Re-ordering of initialization in gam.fit4 to avoid possible failure of
dev.resids call before initialization.
* trap in fast.REML.fit for situation in which all smoothing parameters
satisfy conditions for indefinite convergence on entry, with an
immediate warning that this probably indicates iteration divergence (of bam).
* "bs" basis modified to allow easier control of the interval over which the
spline penalty applies which in turn allows more sensible control of
extrapolation behaviour, when this is unavoidable.
* Fix in uniquecombs - revised faster code (from 1.8-13) could occasionally
generate false matches between different input combinations for integer
variables or factors. Thanks to Rohan Sadler for reporting the issue that
uncovered this.
* A very bad initial model for uninformative data could lead to a negative
fletcher estimate of the scale parameter and optimizer failure - fixed.
* "fREML" allowed in sp.vcov so that it works with bam fitted models.
* 2 occurances of 'return' replaced by (correct) return().
1.8-16
* slightly improved intial value heuristics for overlapping penalties in
......
Package: mgcv
Version: 1.8-16
Version: 1.8-17
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-11-07 10:22:00 UTC; sw283
Packaged: 2017-02-06 11:03:57 UTC; sw283
Repository: CRAN
Date/Publication: 2016-11-07 19:28:16
Date/Publication: 2017-02-08 21:30:16
a38e3e6c9f0ae98bdcf3d14127347fde *ChangeLog
e051228d6327e0b3d640616d4fa510c3 *DESCRIPTION
fdc45fe0bc09d3a9dea994040574ba9d *ChangeLog
1671dbe5a90499d24659b13d04370628 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
595f7fa74dd678ae4be1f6968e174555 *NAMESPACE
fe1f2bde9713ba97eeb59d928310f628 *NAMESPACE
1aee82e17a9feb6166dbe107659029c9 *R/bam.r
87933a9e2845f871334d96f537ee11e9 *R/coxph.r
a25f5a55d9d8b4ada56ac0df3ec99a18 *R/efam.r
4b370253beb1eda19d389132ade30212 *R/fast-REML.r
cb165a0a11de9686910fb27b0d6450ec *R/gam.fit3.r
7734ba62022328a50b771e721df98a72 *R/gam.fit4.r
97c5e0261e26e0aedd8127f9ed41f817 *R/efam.r
deaf0c6a0e39d42a9bad487772eb5904 *R/fast-REML.r
dbf376ceeec4b343d15daa90fe5a81fb *R/gam.fit3.r
9f73b1e3e9781d8ee26a46c1a2127c73 *R/gam.fit4.r
1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
66257913c556f135657b7c12b6ed733d *R/gamlss.r
a6b6a5e70c7889370d33876e08396a09 *R/gamlss.r
cceac26b3d01513f8f87eacae91ddae0 *R/gamm.r
10facb791e4cfd123d183f05660119c6 *R/jagam.r
8cebce4d68d4c45d661073756ca7b202 *R/mgcv.r
088976e7d2ee3df696d09b3fcdfdd563 *R/mgcv.r
2feca0dc9d354de7bc707c67a5891364 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r
5d07d6e4f75a64b6d04cda9e23278d70 *R/plots.r
cf3543d1d4b7fe6fb8576f681b345722 *R/smooth.r
9a4db620ad080f012b56134c1aba10bb *R/smooth.r
666d7fd36fda68b928993d5388b0d7fc *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
074adce4131eced4cc71c67ad1d63c52 *inst/CITATION
7279083a0b1eb3505c3e8105c7c5d2e7 *inst/po/de/LC_MESSAGES/R-mgcv.mo
e74e6d792e3dc69f7a20c68c16097145 *inst/CITATION
342b1840b6a569b8d8d4a3355a0d2bfa *inst/po/de/LC_MESSAGES/R-mgcv.mo
fe9d11e3087789da149e3688309df670 *inst/po/de/LC_MESSAGES/mgcv.mo
14104d65bd4a0005bbc5a3932a37ae8f *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
6b8ae16fa2affd2676d2a47c61885dd1 *inst/po/en@quot/LC_MESSAGES/mgcv.mo
a78bad32fe0b74e954a15ccd98f3a995 *inst/po/fr/LC_MESSAGES/R-mgcv.mo
72ccd181324c1acd4671f69311ed494a *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
dcb794aad0b69ab80cb6877510c3a925 *inst/po/en@quot/LC_MESSAGES/mgcv.mo
72ffab2a2903e1efa30cc81424bce43d *inst/po/fr/LC_MESSAGES/R-mgcv.mo
418bef2f2c1ed07bff6bbdb6884d2858 *inst/po/fr/LC_MESSAGES/mgcv.mo
1553db580cf9e89f3876fe888f9e4e94 *inst/po/ko/LC_MESSAGES/R-mgcv.mo
fbe9f7d7c30fbcebf97c28bf17c2ef50 *inst/po/ko/LC_MESSAGES/R-mgcv.mo
e6196a86ad3a8e42df84242a861d362a *inst/po/ko/LC_MESSAGES/mgcv.mo
6291a3775233a0315dd626e50fc89f2f *inst/po/pl/LC_MESSAGES/R-mgcv.mo
5c8c8783edc20ee9196bb17bb86c420d *inst/po/pl/LC_MESSAGES/R-mgcv.mo
07e822258166c032ff9f1a4e96025841 *inst/po/pl/LC_MESSAGES/mgcv.mo
c574fe1ca9d55a9818d308906f16d16e *man/Beta.Rd
5bf12ddc0dab9daae72271b96a15c539 *man/Predict.matrix.Rd
......@@ -73,6 +73,8 @@ b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
eb8648cc6b3c9374b899abd2f2b12f7b *man/gam2objective.Rd
717401fd7efa3b39d90418a5d1d0c216 *man/gamObject.Rd
a2593990d6a87f7b783d0435062dec02 *man/gamSim.Rd
e5d2541f32dab56972f58b0773eba50c *man/gamlss.etamu.Rd
c7f140d128d1d1d76909499900faf49e *man/gamlss.gH.Rd
351009f7345aaff40f3377e77f4ce7ad *man/gamm.Rd
ec2d1b2aa87cc3f8424e07b6dc0340d5 *man/gaulss.Rd
398a5c12285401c1d37a8edb58780bc3 *man/get.var.Rd
......@@ -83,7 +85,7 @@ faf325db054dce698286f87e09ddc202 *man/gevlss.Rd
2f222eeeb3d7bc42f93869bf8c2af58a *man/influence.gam.Rd
39b9de9dbac7d9dc5c849e1a37def675 *man/initial.sp.Rd
2a37ae59a9f9f5a0a58af45947eca524 *man/interpret.gam.Rd
557a6b7a01d8c0ca7d4db3d9cf329ce6 *man/jagam.Rd
29f94d99b9a3decd4d71618eb839a7f7 *man/jagam.Rd
87d942b17bee05bb662270b894b183e6 *man/ldTweedie.Rd
58e73ac26b93dc9d28bb27c8699e12cf *man/linear.functional.terms.Rd
93035193b0faa32700e1421ce8c1e9f6 *man/logLik.gam.Rd
......@@ -129,7 +131,7 @@ d515e51ec98d73af6166f7b31aeaba9b *man/scat.Rd
6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
8020154bd5c709d11f0e7cf043df886d *man/smooth.construct.Rd
4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.smooth.spec.Rd
9590099266d7f6fafd6bdd1e5c3f61da *man/smooth.construct.bs.smooth.spec.Rd
24927a9c5ce97da988ac8ff3299c533a *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
......@@ -141,7 +143,7 @@ b45d8e71bda4ceca8203dffea577e441 *man/smooth.construct.re.smooth.spec.Rd
0bfe981f2c3e6ea5b8d5372076ccde53 *man/smooth.construct.sos.smooth.spec.Rd
3cb4e59f915c8d64b90754eaeeb5a86f *man/smooth.construct.t2.smooth.spec.Rd
8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
3e6d88ef6a8ab21bd6f120120602dcf6 *man/smooth.construct.tp.smooth.spec.Rd
a088e69cf148d07a78ce95a69759c95c *man/smooth.construct.tp.smooth.spec.Rd
d4083ff900aa69fa07610e0af2a2987b *man/smooth.terms.Rd
ca676e32c30ddaedbce9f706327cc55b *man/smooth2random.Rd
844f9653d74441293d05a24dd3e2876a *man/smoothCon.Rd
......@@ -153,35 +155,36 @@ a0b0988dba55cca5b4b970e035e3c749 *man/t2.Rd
a27690f33b9a7bd56d9c1779c64896cc *man/te.Rd
6eebb6ef90374ee09453d6da6449ed79 *man/tensor.prod.model.matrix.Rd
f22f1cee0ff2b70628846d1d0f8e9a66 *man/trichol.Rd
94154ff18af819a7bb83919ee10db0de *man/uniquecombs.Rd
87e6b4437d00fab4fc814f4cefa3795c *man/trind.generator.Rd
96c48dd705710f639d76a0d0cc3fb128 *man/uniquecombs.Rd
a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
07a73758156dfa580c6e92edd34b0654 *man/ziP.Rd
8bc98d4cb86d851ea0970d68799522cb *man/ziplss.Rd
92cee4501729e5715f0bb95b498f8425 *po/R-de.po
f32c7932ef429aa3113007140da82221 *po/R-de.po
0bdfcf98961b0d52b60f806dc1dba77e *po/R-en@quot.po
1a3e7a9e3d2f18c8df0ccd49b66bf26d *po/R-fr.po
6855e193db0dd078d534c913c21a0398 *po/R-ko.po
461a961615a7c780d6ce5d9694b112d9 *po/R-mgcv.pot
dcca72cd51cd749f5c54478e31abb706 *po/R-pl.po
18713570046eecabea9478b9b864a2cf *po/R-fr.po
d24a793282e6f43c9b67235e46102199 *po/R-ko.po
753bc70ea9b8cb7c06f658a1eb81ce36 *po/R-mgcv.pot
1bb61f6df62fc1787864f66933757701 *po/R-pl.po
8c33d89a914170dbc9f7c5fe598d2135 *po/de.po
93f72334356fe6f05a64e567efd35c8e *po/en@quot.po
b3dfaf74ca2d76ca26eec986a14f5584 *po/fr.po
9116fc041ab458e49b3e498f8c0ac0d9 *po/ko.po
24f393ff94fa39c8e66250eb0d0e2fcd *po/mgcv.pot
cd15f571b582fdd31e597083a5fe66ff *po/mgcv.pot
ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
342aa30c8f6f1e99ffa2576a6f29d7ce *src/coxph.c
555f6948880bff0b6fa23eeb51609c1c *src/discrete.c
9721ea07a126278eacd041c814161968 *src/gdi.c
0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c
dba99f7d7cc412dd9255f6307c8c7fa7 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
acf1a2cff05d68856f2b6acee1e63cf7 *src/init.c
5d9a48e07b438a7c3b32c94fe67ae30c *src/magic.c
195fa3d343a58a362dfac0ce75b5eed7 *src/mat.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c
16b3a72d8bf608026f45de356cf37e65 *src/mat.c
0545dabf3524a110d616ea5e6373295d *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
48cde0e19d5dd54b131ba66c777c0ec2 *src/mgcv.c
3e6fb646794403a7ad4d375f286117f0 *src/mgcv.h
473bc0f430ec18b53d4649c7276c2962 *src/mgcv.h
97e3717e95a70b1470b4c3071e144d17 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c
......
......@@ -8,7 +8,9 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
fs.test,fs.boundary,gam, gam2derivative,
gam2objective,
gamm, gam.check, gam.control,gam.fit3,
gam.fit, gam.outer,gam.vcomp, gamSim ,
gam.fit,
gamlss.etamu,gamlss.gH,
gam.outer,gam.vcomp, gamSim ,
gaulss,gam.side,get.var,gevlss,
influence.gam,
in.out,inSide,interpret.gam,initial.sp,
......@@ -68,7 +70,8 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
summary.gam,sp.vcov,
spasm.construct,spasm.sp,spasm.smooth,
t2,te,ti,tensor.prod.model.matrix,tensor.prod.penalties,
trichol,Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
trichol,trind.generator,
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,
......
......@@ -439,7 +439,7 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
preinitialize <- expression({
ocat.ini <- function(R,y) {
## initialize theta from raw counts in each category
if (R<3) return
if (R<3) return()
y <- c(1:R,y) ## make sure there is *something* in each class
p <- cumsum(tabulate(y[is.finite(y)])/length(y[is.finite(y)]))
eta <- if (p[1]==0) 5 else -1 - log(p[1]/(1-p[1])) ## mean of latent
......
......@@ -941,7 +941,13 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
## idea in following is only to exclude terms with zero first and second derivative
## from optimization, as it is only these that slow things down if included...
uconv.ind <- (abs(grad) > reml.scale*conv.tol*.1)|(abs(grad2)>reml.scale*conv.tol*.1)
## if all appear converged at this stage, then there is probably something wrong,
## but reset anyway to see if situation can be recovered. If we don't do this then
## need to abort immediately, otherwise fails trying to eigen a 0 by 0 matrix
if (sum(uconv.ind)==0) {
warning("Possible divergence detected in fast.REML.fit",call.=FALSE,immediate.=TRUE)
uconv.ind <- rep(TRUE,length(grad))
}
step.failed <- FALSE
for (iter in 1:200) { ## the Newton loop
## Work only with unconverged (much quicker under indefiniteness)
......
......@@ -663,7 +663,8 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
pearson <- sum(weights*(y-mu)^2/family$variance(mu))
scale.est <- (pearson+dev.extra)/(n.true-trA)
if (control$scale.est%in%c("fletcher","Fletcher")) { ## Apply Fletcher (2012) correction
s.bar = mean(family$dvar(mu)*(y-mu)*sqrt(weights)/family$variance(mu))
## note limited to 10 times Pearson...
s.bar = max(-.9,mean(family$dvar(mu)*(y-mu)*sqrt(weights)/family$variance(mu)))
if (is.finite(s.bar)) scale.est <- scale.est/(1+s.bar)
}
} else { ## use the deviance estimator
......@@ -1927,25 +1928,30 @@ bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
### Derivative testing code. Not usually called and not part of BFGS...
ok <- check.derivs
while (ok) { ## derivative testing
deriv <- 1
#deriv <- 1
ok <- FALSE ## set to TRUE to re-run (e.g. with different eps)
deriv.check(x=X, y=y, sp=L%*%lsp+lsp0, Eb=Eb,UrS=UrS,
offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=deriv,
offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=1,
control=control,gamma=gamma,scale=scale,
printWarn=FALSE,mustart=mustart,start=start,
scoreType=scoreType,eps=eps,null.coef=null.coef,Sl=Sl,...)
fdH <- b$dH
fdb.dr <- b$db.drho*0
## deal with fact that deriv might be 0...
bb <- if (deriv==1) b else 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=1,
control=control,gamma=gamma,scale=scale,printWarn=FALSE,start=prev$start,
mustart=prev$mustart,scoreType=scoreType,null.coef=null.coef,
pearson.extra=pearson.extra,dev.extra=dev.extra,n.true=n.true,Sl=Sl,...)
fdH <- bb$dH
fdb.dr <- bb$db.drho*0
for (j in 1:length(lsp)) { ## check dH and db.drho
lsp1 <- lsp;lsp1[j] <- lsp[j] + eps
ba <- 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=deriv,
offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=1,
control=control,gamma=gamma,scale=scale,printWarn=FALSE,start=prev$start,
mustart=prev$mustart,scoreType=scoreType,null.coef=null.coef,
pearson.extra=pearson.extra,dev.extra=dev.extra,n.true=n.true,Sl=Sl,...)
fdH[[j]] <- (ba$H - b$H)/eps
fdb.dr[,j] <- (ba$coefficients - b$coefficients)/eps
fdH[[j]] <- (ba$H - bb$H)/eps
fdb.dr[,j] <- (ba$coefficients - bb$coefficients)/eps
}
}
### end of derivative testing. BFGS code resumes...
......
......@@ -278,23 +278,34 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
## need an initial `null deviance' to test for initial divergence...
## if (!is.null(start)) null.coef <- start - can be on edge of feasible - not good
null.eta <- as.numeric(x%*%null.coef + as.numeric(offset))
old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef
if (!is.null(start)) { ## check it's at least better than null.coef
pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start
if (pdev>old.pdev) start <- mustart <- etastart <- NULL
}
#old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef
#if (!is.null(start)) { ## check it's at least better than null.coef
# pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start
# if (pdev>old.pdev) start <- mustart <- etastart <- NULL
#}
## call the families initialization code...
if (is.null(mustart)) {
eval(family$initialize)
mukeep <- NULL
} else {
mukeep <- mustart
eval(family$initialize)
mustart <- mukeep
#mustart <- mukeep
}
old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef
if (!is.null(start)) { ## check it's at least better than null.coef
pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start
if (pdev>old.pdev) start <- mukeep <- etastart <- NULL
}
if (!is.null(mukeep)) mustart <- mukeep
## and now finalize initialization of mu and eta...
eta <- if (!is.null(etastart)) etastart
......
## (c) Simon N. Wood (2013,2014) distributed under GPL2
## (c) Simon N. Wood (2013-2016) distributed under GPL2
## Code for the gamlss families.
## idea is that there are standard functions converting
## derivatives w.r.t. mu to derivatives w.r.t. eta, given
......@@ -186,8 +186,8 @@ gamlss.gH0 <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=
## fh is a factorization of the penalized hessian, while D contains the corresponding
## Diagonal pre-conditioning weights.
## deriv: 0 - just grad and Hess
## 1 - diagonal of first deriv of Hess
## 2 - first deriv of Hess
## 1 - diagonal of first deriv of Hess wrt sps
## 2 - first deriv of Hess wrt sps
## 3 - everything.
K <- length(jj)
p <- ncol(X);n <- nrow(X)
......@@ -304,7 +304,11 @@ gamlss.gH0 <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=
}
} ## if deriv>2
list(lb=lb,lbb=lbb,d1H=d1H,d2H=d2H,trHid2H=trHid2H)
list(lb=lb, ## grad w.r.t. coefs
lbb=lbb, ## hess w.r.t. coefs, H
d1H=d1H, ## grad of H wrt sps
## d2H=d2H,
trHid2H=trHid2H) ## tr(H^{-1}d^2H/dspsp)
} ## end of gamlss.gH0
......
......@@ -44,9 +44,19 @@ rig <- function(n,mean,scale) {
## inverse guassian deviates generated by algorithm 5.7 of
## Gentle, 2003. scale = 1/lambda.
if (length(n)>1) n <- length(n)
y <- rnorm(n)^2
mu2 <- 0*y + mean^2 ## y is there to ensure mu2 is a vector
x <- mean + 0.5*scale*(mu2*y - mean*sqrt(4*mean*y/scale + mu2*y^2))
x <- y <- rnorm(n)^2
mys <- mean*scale*y
mu <- 0*y + mean ## y is there to ensure mu is a vector
mu2 <- mu^2;
ind <- mys < .Machine$double.eps^-.5 ## cut off for tail computation
x[ind] <- mu[ind]*(1 + 0.5*(mys[ind] - sqrt(mys[ind]*4+mys[ind]^2)))
x[!ind] <- mu[!ind]/mys[!ind] ## tail term (derived from Taylor of sqrt(1+eps) etc)
#my <- mean*y; sc <- 0*y + scale
#ind <- my > 1 ## cancellation error can be severe, without splitting
#x[!ind] <- mu[!ind]*(1 + 0.5*sc[!ind]*(my[!ind] - sqrt(4*my[!ind]/sc[!ind] + my[!ind]^2)))
## really the sqrt in the next term should be expanded beyond first order and then
## worked on - otherwise too many exact zeros?
#x[ind] <- pmax(0,mu[ind]*(1+my[ind]*.5*sc[ind]*(1-sqrt(1+ 4/(sc[ind]*my[ind])))))
ind <- runif(n) > mean/(mean+x)
x[ind] <- mu2[ind]/x[ind]
x ## E(x) = mean; var(x) = scale*mean^3
......@@ -1580,7 +1590,7 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,start=NU
#}
## make sure optimizer appropriate for available derivatives
if (!is.null(G$family$available.derivs)) {
if (G$family$available.derivs==1) optimizer <- c("outer","bfgs")
if (G$family$available.deriv==1 && optimizer[1]!="efs") optimizer <- c("outer","bfgs")
if (G$family$available.derivs==0) optimizer <- "efs"
}
}
......@@ -3863,7 +3873,7 @@ cooks.distance.gam <- function(model,...)
sp.vcov <- function(x) {
## 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")&&!is.null(x$outer.info$hess)) {
if (x$method%in%c("ML","P-ML","REML","P-REML","fREML")&&!is.null(x$outer.info$hess)) {
return(solve(x$outer.info$hess))
} else return(NULL)
}
......@@ -3873,7 +3883,7 @@ gam.vcomp <- function(x,rescale=TRUE,conf.lev=.95) {
## in a fitted `gam' object.
if (!inherits(x,"gam")) stop("requires an object of class gam")
if (!is.null(x$reml.scale)&&is.finite(x$reml.scale)) scale <- x$reml.scale else scale <- x$sig2
if (length(x$sp)==0) return
if (length(x$sp)==0) return()
if (rescale) { ## undo any rescaling of S[[i]] that may have been done
m <- length(x$smooth)
if (is.null(x$paraPen)) {
......
......@@ -154,6 +154,87 @@ mono.con<-function(x,up=TRUE,lower=NA,upper=NA)
uniquecombs <- function(x,ordered=FALSE) {
## takes matrix x and counts up unique rows
## `unique' now does this in R
if (is.null(x)) stop("x is null")
if (is.null(nrow(x))||is.null(ncol(x))) x <- data.frame(x)
recheck <- FALSE
if (inherits(x,"data.frame")) {
xoo <- xo <- x
## reset character, logical and factor to numeric, to guarantee that text versions of labels
## are unique iff rows are unique (otherwise labels containing "*" could in principle
## fool it).
is.char <- rep(FALSE,length(x))
for (i in 1:length(x)) {
if (is.character(xo[[i]])) {
is.char[i] <- TRUE
xo[[i]] <- as.factor(xo[[i]])
}
if (is.factor(xo[[i]])||is.logical(xo[[i]])) x[[i]] <- as.numeric(xo[[i]])
if (!is.numeric(x[[i]])) recheck <- TRUE ## input contains unknown type cols
}
#x <- data.matrix(xo) ## ensure all data are numeric
} else xo <- NULL
if (ncol(x)==1) { ## faster to use R
xu <- if (ordered) sort(unique(x[,1])) else unique(x[,1])
ind <- match(x[,1],xu)
if (is.null(xo)) x <- matrix(xu,ncol=1,nrow=length(xu)) else {
x <- data.frame(xu)
names(x) <- names(xo)
}
} else { ## no R equivalent that directly yields indices
if (ordered) {
chloc <- Sys.getlocale("LC_CTYPE")
Sys.setlocale("LC_CTYPE","C")
}
## txt <- paste("paste0(",paste("x[,",1:ncol(x),"]",sep="",collapse=","),")",sep="")
## ... this can produce duplicate labels e.g. x[,1] = c(1,11), x[,2] = c(12,2)...
## solution is to insert separator not present in representation of a number (any
## factor codes are already converted to numeric by data.matrix call above.)
txt <- paste("paste0(",paste("x[,",1:ncol(x),"]",sep="",collapse=",\"*\","),")",sep="")
xt <- eval(parse(text=txt)) ## text representation of rows
dup <- duplicated(xt) ## identify duplicates
xtu <- xt[!dup] ## unique text rows
x <- x[!dup,] ## unique rows in original format
#ordered <- FALSE
if (ordered) { ## return unique in same order regardless of entry order
## ordering of character based labels is locale dependent
## so that e.g. running the same code interactively and via
## R CMD check can give different answers.
coloc <- Sys.getlocale("LC_COLLATE")
Sys.setlocale("LC_COLLATE","C")
ii <- order(xtu)
Sys.setlocale("LC_COLLATE",coloc)
Sys.setlocale("LC_CTYPE",chloc)
xtu <- xtu[ii]
x <- x[ii,]
}
ind <- match(xt,xtu) ## index each row to the unique duplicate deleted set
}
if (!is.null(xo)) { ## original was a data.frame
x <- as.data.frame(x)
names(x) <- names(xo)
for (i in 1:ncol(xo)) {
if (is.factor(xo[,i])) { ## may need to reset factors to factors
xoi <- levels(xo[,i])
x[,i] <- if (is.ordered(xo[,i])) ordered(x[,i],levels=1:length(xoi),labels=xoi) else
factor(x[,i],levels=1:length(xoi),labels=xoi)
contrasts(x[,i]) <- contrasts(xo[,i])
}
if (is.char[i]) x[,i] <- as.character(x[,i])
if (is.logical(xo[,i])) x[,i] <- as.logical(x[,i])
}
}
if (recheck) {
if (all.equal(xoo,x[ind,],check.attributes=FALSE)!=TRUE) warning("uniquecombs has not worked properly")
}
attr(x,"index") <- ind
x
} ## uniquecombs
uniquecombs0 <- function(x,ordered=FALSE) {
## takes matrix x and counts up unique rows
## `unique' now does this in R
if (is.null(x)) stop("x is null")
if (is.null(nrow(x))||is.null(ncol(x))) x <- data.frame(x)
......@@ -170,7 +251,11 @@ uniquecombs <- function(x,ordered=FALSE) {
chloc <- Sys.getlocale("LC_CTYPE")
Sys.setlocale("LC_CTYPE","C")
}
txt <- paste("paste0(",paste("x[,",1:ncol(x),"]",sep="",collapse=","),")",sep="")
## txt <- paste("paste0(",paste("x[,",1:ncol(x),"]",sep="",collapse=","),")",sep="")
## ... this can produce duplicate labels e.g. x[,1] = c(1,11), x[,2] = c(12,2)...
## solution is to insert separator not present in representation of a number (any
## factor codes are already converted to numeric by data.matrix call above.)
txt <- paste("paste0(",paste("x[,",1:ncol(x),"]",sep="",collapse=",\":\","),")",sep="")
xt <- eval(parse(text=txt)) ## text representation of rows
dup <- duplicated(xt) ## identify duplicates
xtu <- xt[!dup] ## unique text rows
......@@ -203,7 +288,7 @@ uniquecombs <- function(x,ordered=FALSE) {
}
attr(x,"index") <- ind
x
} ## uniquecombs
} ## uniquecombs0
cSplineDes <- function (x, knots, ord = 4,derivs=0)
{ ## cyclic version of spline design...
......@@ -626,15 +711,15 @@ tensor.prod.penalties <- function(S)
smooth.construct.tensor.smooth.spec <- function(object,data,knots)
smooth.construct.tensor.smooth.spec <- function(object,data,knots) {
## the constructor for a tensor product basis object
{ inter <- object$inter ## signal generation of a pure interaction
inter <- object$inter ## signal generation of a pure interaction
m <- length(object$margin) # number of marginal bases
if (inter) {
object$mc <- if (is.null(object$mc)) rep(TRUE,m) else as.logical(object$mc)
if (inter) { ## interaction term so at least some marginals subject to constraint
object$mc <- if (is.null(object$mc)) rep(TRUE,m) else as.logical(object$mc) ## which marginals to constrain
object$sparse.cons <- if (is.null(object$sparse.cons)) rep(0,m) else object$sparse.cons
} else {
object$mc <- rep(FALSE,m)
object$mc <- rep(FALSE,m) ## all marginals unconstrained
}
Xm <- list();Sm<-list();nr<-r<-d<-array(0,m)
C <- NULL
......@@ -675,8 +760,8 @@ smooth.construct.tensor.smooth.spec <- function(object,data,knots)
km <- which(mono)
g <- list(); for (i in 1:length(km)) g[[i]] <- object$margin[[km[i]]]$g.index
for (i in 1:length(object$margin)) {
d <- ncol(object$margin[[i]]$X)
for (j in length(km)) if (i!=km[j]) g[[j]] <- if (i > km[j]) rep(g[[j]],each=d) else rep(g[[j]],d)
dx <- ncol(object$margin[[i]]$X)
for (j in length(km)) if (i!=km[j]) g[[j]] <- if (i > km[j]) rep(g[[j]],each=dx) else rep(g[[j]],dx)
}
object$g.index <- as.logical(rowSums(matrix(unlist(g),length(g[[1]]),length(g))))
}
......@@ -1771,12 +1856,21 @@ smooth.construct.bs.smooth.spec <- function(object,data,knots) {
if (length(k)==2) {
xl <- min(k);xu <- max(k);
if (xl>min(x)||xu<max(x)) stop("knot range does not include data")
}
if (is.null(k)||length(k)==2) {
}
if (!is.null(k)&&length(k)==4&&length(k)<nk+2*m[1]) {
## 4 knots supplied: lower prediction limit, lower data limit,
## upper data limit, upper prediction limit
k <- sort(k)
dx <- (k[4]-k[1])/(nk-1)
ko <- c(k[1]-dx*m[1],k[4]+dx*m[1]) ## limits for outer knots
k <- c(seq(ko[1],k[1],length=m[1]+1),
seq(k[2],k[3],length=max(0,nk-2)),
seq(k[4],ko[2],length=m[1]+1))
} else if (is.null(k)||length(k)==2) {
xr <- xu - xl # data limits and range
xl <- xl-xr*0.001;xu <- xu+xr*0.001;dx <- (xu-xl)/(nk-1)
k <- seq(xl-dx*(m[1]),xu+dx*(m[1]),length=nk+2*m[1])
k <- seq(xl-dx*m[1],xu+dx*m[1],length=nk+2*m[1])
} else {
if (length(k)!=nk+2*m[1])
stop(paste("there should be ",nk+2*m[1]," supplied knots"))
......@@ -3311,7 +3405,7 @@ ExtractData <- function(object,data,knots) {
knt[[object$term[i]]] <- get.var(object$term[i],knots)
}
names(dat) <- object$term;m <- length(object$term)
names(dat) <- object$term; m <- length(object$term)
if (!is.null(attr(dat[[1]],"matrix"))) { ## strip down to unique covariate combinations
n <- length(dat[[1]])
X <- matrix(unlist(dat),n,m)
......
mgcv (1.8-17-1) unstable; urgency=medium
* New upstream release
-- Dirk Eddelbuettel <edd@debian.org> Wed, 08 Feb 2017 20:10:07 -0600
mgcv (1.8-16-1) unstable; urgency=medium
* New upstream release
......
citHeader("2011 for generalized additive model method;
2004 for strictly additive GCV based model method and basics of gamm;
2006 for overview; 2003 for thin plate regression splines;
2000 is the original method, now superceded.")
2016 for jagam.")
citEntry(
entry="Article",
bibentry(
bibtype="Article",
title="Fast stable restricted maximum likelihood and marginal
likelihood estimation of semiparametric generalized linear models",
journal="Journal of the Royal Statistical Society (B)",
......@@ -17,8 +17,8 @@ likelihood estimation of semiparametric generalized linear models",
and marginal likelihood estimation of semiparametric generalized linear
models. Journal of the Royal Statistical Society (B) 73(1):3-36" )
citEntry(
entry="Article",
bibentry(
bibtype="Article",
title= "Stable and efficient multiple smoothing parameter estimation for
generalized additive models",
journal="Journal of the American Statistical Association",
......@@ -32,8 +32,8 @@ citEntry(
American Statistical Association. 99:673-686."
)
citEntry(
entry="Book",
bibentry(
bibtype="Book",
title="Generalized Additive Models: An Introduction with R",
year="2006",
author="S.N Wood",
......@@ -42,8 +42,8 @@ textVersion="Wood, S.N. (2006) Generalized Additive Models: An
Introduction with R. Chapman and Hall/CRC. "
)
citEntry(
entry="Article",
bibentry(
bibtype="Article",
title="Thin-plate regression splines",
journal="Journal of the Royal Statistical Society (B)",
volume= "65",
......@@ -55,16 +55,19 @@ citEntry(
Journal of the Royal Statistical Society (B) 65(1):95-114." )
citEntry(
entry="Article",
title="Modelling and smoothing parameter estimation with
multiple quadratic penalties",
journal="Journal of the Royal Statistical Society (B)",
volume= "62",
number="2",
pages="413-428",
year="2000",
author="S. N. Wood",
textVersion="Wood, S.N. (2000) Modelling and smoothing parameter
estimation with multiple quadratic penalties. Journal of the Royal
Statistical Society (B) 62(2):413-428. " )
bibentry(bibtype = "Article",
title = "Just Another Gibbs Additive Modeler: Interfacing {JAGS} and {mgcv}",
author = "S. N. Wood",
journal = "Journal of Statistical Software",
year = "2016",
volume = "75",
number = "7",
pages = "1--15",
doi = "10.18637/jss.v075.i07",
textVersion =
paste("Simon N. Wood (2016).",
"Just Another Gibbs Additive Modeler: Interfacing JAGS and mgcv.",
"Journal of Statistical Software, 75(7), 1-15.",
"doi:10.18637/jss.v075.i07")
)
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
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
\name{gamlss.etamu}
\alias{gamlss.etamu}
\title{Transform derivatives wrt mu to derivatives wrt linear predictor}
\usage{
gamlss.etamu(l1, l2, l3 = NULL, l4 = NULL, ig1, g2, g3 = NULL,
g4 = NULL, i2, i3 = NULL, i4 = NULL, deriv = 0)