Commit a2c1c4ac authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Debian changes 1.8-23-1

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

  * New upstream release
 
  * debian/control: Set Build-Depends: to current R version
  * debian/control: Set Standards-Version: to current version 
parents 6e1e032b d8806868
......@@ -3,6 +3,58 @@
Currently deprecated and liable to be removed:
- gam performance iteration (1.8-19, Sept 2017)
- argument 'pers' in plot.gam (1.8-23, Nov 2018)
1.8-23
* default plot methods added for smooths of 3 and 4 variables.
* The `gamma' control parameter for gam and bam can now be used with RE/ML
smoothness selection, not just GCV/AIC. Essentially smoothing parameters
are chosen as if the sample size was n/gamma instead of n.
* The "bs" basis now allows multiple penalties of different orders on the same
spline. e.g. s(x,bs="bs",m=c(3,2,0)). See ?b.spline.
* bam(...,discrete=TRUE) can now use the smooth 'id' mechanism to link
smoothing parameters, but note the method constraint that the linked bases
are not forced to be identical in this case (unlike other fitting methods).
* summary.gam now allows random effects tests to be skipped (in some large
models the test is costly and uninteresting).
* 'interpret.gam0' modified so that masked 's', 'te', etc from other packages
does not cause failure.
* coxph fix of prediction bug introduced with stratified model (thanks
Giampiero Marra)
* bam(...,discrete=TRUE) fix to handle nested matrix arguments to smooths.
* bam(...,discrete=TRUE) fix to by variable handling with fs and re smooths
which could fail during re-representation as tensor smooths (for
discretization purposes).
* bam extended family extension had introduced a bug in null deviance
computation for Gaussian additive case when using methods other than fREML
or GCV.Cp. Fixed.
* bam(...,discrete=TRUE) now defaults to discrete=FALSE if there are no
smooths, rather than failing.
* bam was reporting wrong family under some smoothing parameter selection
methods (not default).
* null deviance computation improved for extended families. Previous version
used an approximation valid for most families, and corrected the rest - now
replaced with exact computations for all cases.
* scat initialization tweaked to avoid -ve def problems at start.
* paraPen handling in bam was broken - fixed.
* slight adjustment to sp initialization for extended families - use observed
information in weights if possible.
1.8-22
......
Package: mgcv
Version: 1.8-22
Version: 1.8-23
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with Automatic Smoothness
......@@ -18,6 +18,6 @@ LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2017-09-18 10:38:41 UTC; sw283
Packaged: 2018-01-08 15:15:43 UTC; sw283
Repository: CRAN
Date/Publication: 2017-09-19 00:27:56 UTC
Date/Publication: 2018-01-15 00:23:36 UTC
bc8bc8be9cec065f0c7ec0cf670516f2 *ChangeLog
ad2a852c227dd80d86b24784d12dba42 *DESCRIPTION
6e3790e8df181acf2ecf469373521368 *ChangeLog
a01208f1c73866153328f2ce262e28d9 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
9fa89dc9361930dca536d4e011e4e605 *NAMESPACE
b051474f30e8e779f2440e8ecd5bbd51 *R/bam.r
7b419683b0948cf6da009f614078fe90 *R/coxph.r
777a0d67a1f7fa14bf87bc312064061b *R/efam.r
dfdb821247da3e780de0d4599b88735d *R/fast-REML.r
0a9cc97524994d1563b93769284266c4 *R/gam.fit3.r
f0fd93b43b36abf05e6666b4a71c20fa *R/gam.fit4.r
20feb98661e2c2d10bf2d3337af57cb3 *NAMESPACE
256ca3ddb2f799583b97d24ac136e579 *R/bam.r
28ef756f5a8c6bc10da15ad22c15eb1b *R/coxph.r
85f72f767420468e5d0a9b8608128607 *R/efam.r
276847b54ce19aea6c10d5922c31ea8c *R/fast-REML.r
f52b79e696006be12556df1140d6e294 *R/gam.fit3.r
374afe5d0d7441db1d823d239548a1b1 *R/gam.fit4.r
1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
b99d9de0c46a230e081b9bbbbbb3c617 *R/gamlss.r
c934ba653cb1a904132b3f15a01e41c5 *R/gamm.r
f2368e22ccf048cf3d1e210e970013dd *R/gamlss.r
857d79084bc67dee97d412267d5cac44 *R/gamm.r
10facb791e4cfd123d183f05660119c6 *R/jagam.r
b9084037d058ef5d787e2818fbc29bae *R/mgcv.r
1c500e2c478843eb16c967a2472e48f5 *R/mgcv.r
2feca0dc9d354de7bc707c67a5891364 *R/misc.r
03772998ab05887f2eeae10dd6efe983 *R/mvam.r
ca13f63d1112d31258eb3c68464fa653 *R/plots.r
bdc529fe735f22e20aa0cde6708cde25 *R/smooth.r
666d7fd36fda68b928993d5388b0d7fc *R/soap.r
aeacfacc81c1726bdaee38e7f173c844 *R/sparse.r
2d0b7c4deaa2ac5d97dbc00fcf8ef395 *R/plots.r
694f20ac853e6cfa817bade3c6b1f518 *R/smooth.r
7398607a9ba7e85276bcf09df0d9344b *R/soap.r
bde1774ce7903cabc57b3196f8872ea8 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
794d1e804e26577500008820b9196f24 *inst/CITATION
e70536903ca327ade1496529ab817b60 *inst/CITATION
d0b290e7efd8b4c47a53dcde665a4d42 *inst/po/de/LC_MESSAGES/R-mgcv.mo
fe9d11e3087789da149e3688309df670 *inst/po/de/LC_MESSAGES/mgcv.mo
78e6bcd08b905ed2adb2ffac48034e06 *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
......@@ -42,7 +42,7 @@ fd0cfd64be579f9fbecdbb7f2b8ec1eb *man/Sl.initial.repara.Rd
60670020425f8749b81a8d8c3f168880 *man/Sl.setup.Rd
69ae63833438a3af2963e39482f1d72f *man/Tweedie.Rd
8087ab00d10b44c99c37f49bf90e19cd *man/anova.gam.Rd
6180a4e9ea206e1a350b993dade0a869 *man/bam.Rd
4e6114b2f161710eb862d6b02e1c6870 *man/bam.Rd
ab5e37c3bf8803de63b63c3bdc5909cd *man/bam.update.Rd
cf5f1ee0aab639c7c4b9b357434f15b2 *man/bandchol.Rd
745cbf31eb14fc1c5916fc634c74d998 *man/bug.reports.mgcv.Rd
......@@ -62,8 +62,8 @@ e75719779e18c723ee1fd17e44e7901b *man/formXtViX.Rd
88888e966394c9f9792d9537341d053c *man/formula.gam.Rd
4da4d585b329769eb44f0c7a6e7dd554 *man/fs.test.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
e2d118b71a00f65912912f629945eb53 *man/gam.Rd
430b6d60be7818b659e2ffe5ec29cd7d *man/gam.check.Rd
1d56a971aa67782fa267f5c6f85b57ec *man/gam.Rd
42d669d0f18eba79234b2e849d724863 *man/gam.check.Rd
8931cd75ddec14d91ec94dec6ba69362 *man/gam.control.Rd
afd2fdc49ac57b4b94ee233fa5da1d64 *man/gam.convergence.Rd
1cf5145859af2263f4e3459f40e1ab23 *man/gam.fit.Rd
......@@ -121,17 +121,17 @@ ee9352ba4c531a8def16deddcab9a9fd *man/pdIdnot.Rd
8bc429d92aa9f58c4c43f2852e1f8123 *man/pdTens.Rd
1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd
931c3aefb8b5e42aa230cfedd281bed1 *man/place.knots.Rd
74b9891dcf56eae77cbd11b04fa0623a *man/plot.gam.Rd
ad49b46ad1751b09ec805b0f06f0e890 *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
1a9d83c9fc67e5f0fc85d66d3112f4ef *man/predict.bam.Rd
6ea139c7a291c7dbf897a3f39afc97c0 *man/predict.gam.Rd
41b1b0883f2419c8e2d2dd913fc7d677 *man/print.gam.Rd
6d0ce4e574fabceffdbedd46c91364cb *man/qq.gam.Rd
22b7dcbc8ff4096365fa98ce56b957c9 *man/rTweedie.Rd
1a533ae9f698437faa076021194c9b1d *man/random.effects.Rd
420e3d95f4c38c73303536ee33bfc98b *man/random.effects.Rd
c523210ae95cb9aaa0aaa1c37da1a4c5 *man/residuals.gam.Rd
3c747a8066bcc28ae706ccf74f903d3e *man/rig.Rd
9f6f46f5c5da080bc82f9aa4685d364a *man/rmvn.Rd
4faef2a628f3c2c43c8a88552ff5a7df *man/rmvn.Rd
1c0b5434ba69062941aa72b262d32ff7 *man/s.Rd
c45d0a4edfa63ff362bd34195b3059ca *man/scat.Rd
898e7cc2def2ee234475e68d0b904b29 *man/sdiag.Rd
......@@ -139,11 +139,11 @@ d54f4042e212fca7704cf8428bdaea38 *man/single.index.Rd
6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
8020154bd5c709d11f0e7cf043df886d *man/smooth.construct.Rd
4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.smooth.spec.Rd
24927a9c5ce97da988ac8ff3299c533a *man/smooth.construct.bs.smooth.spec.Rd
fc7fba34b89fdf29f66744d1fdadda79 *man/smooth.construct.bs.smooth.spec.Rd
2f0463d1aca0b8798da6e681bd4c6e54 *man/smooth.construct.cr.smooth.spec.Rd
f5e6d0f5122f61c336827b3615482157 *man/smooth.construct.ds.smooth.spec.Rd
db75c958cbfb561914a3291ab58b9786 *man/smooth.construct.fs.smooth.spec.Rd
7eba0bd0ecf412e7db5dab49eba724cd *man/smooth.construct.gp.smooth.spec.Rd
92591aadf25c362bed2b07da4adbd8be *man/smooth.construct.gp.smooth.spec.Rd
ba28e87c767445d260afa5e685821d73 *man/smooth.construct.mrf.smooth.spec.Rd
2523b6cefa306210c00f3477853b7f07 *man/smooth.construct.ps.smooth.spec.Rd
b45d8e71bda4ceca8203dffea577e441 *man/smooth.construct.re.smooth.spec.Rd
......@@ -158,7 +158,7 @@ ae5e27524e37d57505754639455f18a5 *man/smooth.terms.Rd
b55a396da77559dac553613146633f97 *man/sp.vcov.Rd
83bd8e097711bf5bd0fff09822743d43 *man/spasm.construct.Rd
b9394812e5398ec95787c65c1325a027 *man/step.gam.Rd
a559de7ef43b63d58ed67faf5316e002 *man/summary.gam.Rd
62a42d898c2f1ccd7a64aef33d07b3a1 *man/summary.gam.Rd
a0b0988dba55cca5b4b970e035e3c749 *man/t2.Rd
d1358e5c7f1f9d9a56072a77787803d2 *man/te.Rd
6eebb6ef90374ee09453d6da6449ed79 *man/tensor.prod.model.matrix.Rd
......@@ -187,19 +187,19 @@ ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po
0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c
f6c4ba80ce1b71be7f4a44fac1b08c28 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
a287f92033406194afe0301ee076fda8 *src/init.c
3fc05fa48bda8202d05926014324404d *src/init.c
a9151b5236852eef9e6947590bfcf88a *src/magic.c
654ff83187dc0f7ef4e085f3348f70d2 *src/mat.c
e4cef7f1753153fbab242d1c4d4f7e7f *src/matrix.c
de37b0972199b796654405efc007f25b *src/matrix.h
8df4b96961491d76989b50856237ee2d *src/mgcv.c
b6278c9a4a32bdc16487a96e6a0045b6 *src/mgcv.h
4a5fe2fc83908d1b868b5dabcf7158ab *src/mgcv.h
97e3717e95a70b1470b4c3071e144d17 *src/misc.c
465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
563938b7bb6504ab10df5376c4360220 *src/qp.c
073a4b5b0bc6e869c5b35478c47facf1 *src/qp.h
d5673b88f6f3d85c62a1337f49abba24 *src/soap.c
dcac8c02b5f9be28d13efc834fc88d55 *src/sparse-smooth.c
7ad1c265233960a282682a6c90e27057 *src/sparse-smooth.c
fe0444bece322bc229e46b3d1c150779 *src/tprs.c
5bd85bf0319a7b7c755cf49c91a7cd94 *src/tprs.h
38e593a85a6fd0bb4fbed836f3361406 *tests/bam.R
......
......@@ -73,7 +73,8 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
totalPenaltySpace,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(grDevices,cm.colors,dev.interactive,devAskNewPage,gray,grey,
heat.colors,terrain.colors,topo.colors,axisTicks)
importFrom(graphics,abline,axis,axTicks,box,contour,hist,image,lines,
mtext, par, persp,plot,points,
polygon,rect,strheight,strwidth,text,title)
......
This diff is collapsed.
......@@ -155,7 +155,7 @@ cox.ph <- function (link = "identity") {
strat <- NULL
}
if (sum(is.na(y))>0) stop("NA times supplied for cox.ph prediction")
X <- X[ii,];y <- y[ii];
X <- X[ii,,drop=FALSE];y <- y[ii];
if (is.null(strat)) {
n <- nrow(X)
oo <- .C("coxpred",as.double(X),t=as.double(y),as.double(beta),as.double(Vb),
......
## (c) Simon N. Wood (ocat, tw, nb, ziP) & Natalya Pya (scat, beta),
## 2013-2017. Released under GPL2.
estimate.theta <- function(theta,family,y,mu,scale=1,wt=1,tol=1e-7) {
estimate.theta <- function(theta,family,y,mu,scale=1,wt=1,tol=1e-7,attachH=FALSE) {
## Simple Newton iteration to estimate theta for an extended family,
## given y and mu. To be iterated with estimation of mu given theta.
## If used within a PIRLS loop then divergence testing of coef update
......@@ -79,9 +79,32 @@ estimate.theta <- function(theta,family,y,mu,scale=1,wt=1,tol=1e-7) {
if (sum(abs(nll$g) > tol*abs(nll$nll))==0) break
} ## main Newton loop
if (step.failed) warning("step failure in theta estimation")
if (attachH) attr(theta,"H") <- nll$H
theta
} ## estimate.theta
find.null.dev <- function(family,y,eta,offset,weights) {
## obtain the null deviance given y, best fit mu and
## prior weights
fnull <- function(gamma,family,y,wt,offset) {
## evaluate deviance for single parameter model
mu <- family$linkinv(gamma+offset)
sum(family$dev.resids(y,mu, wt))
}
mu <- family$linkinv(eta-offset)
mum <- mean(mu*weights)/mean(weights) ## initial value
eta <- family$linkfun(mum) ## work on l.p. scale
deta <- abs(eta)*.1 + 1 ## search interval half width
ok <- FALSE
while (!ok) {
search.int <- c(eta-deta,eta+deta)
op <- optimize(fnull,interval=search.int,family=family,y=y,wt = weights,offset=offset)
if (op$minimum > search.int[1] && op$minimum < search.int[2]) ok <- TRUE else deta <- deta*2
}
op$objective
} ## find.null.dev
## extended families for mgcv, standard components.
## family - name of family character string
## link - name of link character string
......@@ -171,8 +194,11 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
}
theta
}
postproc <- function(family,...) {
postproc <- function(family,y,prior.weights,fitted,linear.predictors,offset,intercept) {
posr <- list()
## null.deviance needs to be corrected...
posr$null.deviance <- find.null.dev(family,y,eta=linear.predictors,offset,prior.weights)
posr$family <-
paste("Ordered Categorical(",paste(round(family$getTheta(TRUE),2),collapse=","),")",sep="")
posr
......@@ -448,6 +474,8 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
ls <- function(y,w,theta,scale) {
## the log saturated likelihood function.
#! actually only first line used since re-def as 0
#vec <- !is.null(attr(theta,"vec.grad"))
#lsth1 <- if (vec) matrix(0,length(y),R-2) else rep(0,R-2)
return(list(ls=0,lsth1=rep(0,R-2),lsth2=matrix(0,R-2,R-2)))
F <- function(x) {
h <- ind <- x > 0; h[ind] <- 1/(exp(-x[ind]) + 1)
......@@ -762,6 +790,7 @@ nb <- function (theta = NULL, link = "log") {
ls <- function(y,w,theta,scale) {
## the log saturated likelihood function.
Theta <- exp(theta)
#vec <- !is.null(attr(theta,"vec.grad")) ## lsth by component?
ylogy <- y;ind <- y>0;ylogy[ind] <- y[ind]*log(y[ind])
term <- (y + Theta) * log(y + Theta) - ylogy +
lgamma(y + 1) - Theta * log(Theta) + lgamma(Theta) -
......@@ -773,6 +802,7 @@ nb <- function (theta = NULL, link = "log") {
psi0.yth <- digamma(yth)
psi0.th <- digamma(Theta)
term <- Theta * (lyth - psi0.yth + psi0.th-theta)
#lsth <- if (vec) -term*w else -sum(term*w)
lsth <- -sum(term*w)
## second deriv wrt theta...
psi1.yth <- trigamma(yth)
......@@ -790,8 +820,9 @@ nb <- function (theta = NULL, link = "log") {
mustart <- y + (y == 0)/6
})
postproc <- function(family,...) {
postproc <- function(family,y,prior.weights,fitted,linear.predictors,offset,intercept){
posr <- list()
posr$null.deviance <- find.null.dev(family,y,eta=linear.predictors,offset,prior.weights)
posr$family <-
paste("Negative Binomial(",round(family$getTheta(TRUE),3),")",sep="")
posr
......@@ -920,6 +951,8 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) {
mup1 <- mu^(-p-1)
r$Dmu3 <- -2 * wt * mup1*p*(y/mu*(p+1) + 1-p)
r$Dmu2th <- 2 * wt * (mup1*y*(1-p*logmu)-(logmu*(1-p)+1)/mup )*dpth1
r$EDmu3 <- -2*wt*p*mup1
r$EDmu2th <- -2*wt*logmu/mup*dpth1
}
if (level>1) { ## whole damn lot
mup2 <- mup1/mu
......@@ -953,8 +986,12 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) {
ls <- function(y, w, theta, scale) {
## evaluate saturated log likelihood + derivs w.r.t. working params and log(scale)
a <- get(".a");b <- get(".b")
LS <- colSums(w * ldTweedie(y, y, rho=log(scale), theta=theta,a=a,b=b))
lsth1 <- c(LS[4],LS[2])
#vec <- !is.null(attr(theta,"vec.grad"))
LS <- w * ldTweedie(y, y, rho=log(scale), theta=theta,a=a,b=b)
#if (vec) lsth1 <- LS[,c(4,2)]
LS <- colSums(LS)
#if (!vec) lsth1 <- c(LS[4],LS[2])
lsth1 <- c(LS[4],LS[2])
lsth2 <- matrix(c(LS[5],LS[6],LS[6],LS[3]),2,2)
list(ls=LS[1],lsth1=lsth1,lsth2=lsth2)
}
......@@ -965,8 +1002,9 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) {
mustart <- y + (y == 0)*.1
})
postproc <- function(family,...) {
postproc <- function(family,y,prior.weights,fitted,linear.predictors,offset,intercept) {
posr <- list()
posr$null.deviance <- find.null.dev(family,y,eta=linear.predictors,offset,prior.weights)
posr$family <-
paste("Tweedie(p=",round(family$getTheta(TRUE),3),")",sep="")
posr
......@@ -1025,7 +1063,7 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
n.theta <- 0 ## signal that there are no theta parameters to estimate
} else iniTheta <- log(-theta) ## initial theta supplied
} else iniTheta <- 0 ## inital log theta value
env <- new.env(parent = .GlobalEnv)
assign(".Theta", iniTheta, envir = env)
assign(".betarEps",eps, envir = env)
......@@ -1369,12 +1407,13 @@ scat <- function (theta = NULL, link = "identity",min.df = 3) {
term <- 2*nu1/sig^2/(nu+3)
n <- length(y)
oo$EDmu2 <- rep(term,n)
if (level>0) { ## quantities needed for first derivatives
nu1nusig2a <- nu1/nusig2a
nu2nu <- nu2/nu
fym <- f*ym; ff1 <- f*f1; f1ym <- f1*ym; fymf1 <- fym*f1
ymsig2a <- ym/sig2a
oo$Dmu2th <- oo$Dmuth <- oo$Dth <- matrix(0,n,2)
oo$EDmu2th <- oo$Dmu2th <- oo$Dmuth <- oo$Dth <- matrix(0,n,2)
oo$Dth[,1] <- 1 * wt * nu2 * (log(a) - fym/nu)
oo$Dth[,2] <- -2 * wt * fym
oo$Dmuth[,1] <- 2 * wt *(f - ymsig2a - fymf1)*nu2nu
......@@ -1382,6 +1421,8 @@ scat <- function (theta = NULL, link = "identity",min.df = 3) {
oo$Dmu3 <- 4 * wt * f * (3/nusig2a - 4*f1^2)
oo$Dmu2th[,1] <- 2* wt * (-nu1nusig2a + 1/sig2a + 5*ff1- 2*f1ym/sig2a - 4*fymf1*f1)*nu2nu
oo$Dmu2th[,2] <- 4*wt*(-nu1nusig2a + ff1*5 - 4*ff1*f1ym)
oo$EDmu3 <- rep(0,n)
oo$EDmu2th <- cbind(4/(sig^2*(nu+3)^2)*exp(theta[1]),-2*oo$EDmu2)
}
if (level>1) { ## whole lot
## nu1nu2 <- nu1*nu2;
......@@ -1432,18 +1473,17 @@ scat <- function (theta = NULL, link = "identity",min.df = 3) {
## the log saturated likelihood function.
## (Note these are correct but do not correspond to NP notes)
if (length(w)==1) w <- rep(w,length(y))
#vec <- !is.null(attr(theta,"vec.grad"))
min.df <- get(".min.df")
nu <- exp(theta[1])+min.df; sig <- exp(theta[2]); nu2 <- nu-min.df;
nu2nu <- nu2/nu; nu12 <- (nu+1)/2
term <- lgamma(nu12) - lgamma(nu/2) - log(sig*(pi*nu)^.5)
ls <- sum(term*w)
## first derivative wrt theta...
lsth <- rep(0,2)
lsth2 <- matrix(0,2,2) ## rep(0, 3)
term <- nu2 * digamma(nu12)/2- nu2 * digamma(nu/2)/2 - 0.5*nu2nu
lsth[1] <- sum(w*term)
lsth[2] <- sum(-1*w)
#lsth <- if (vec) cbind(w*term,-1*w) else c(sum(w*term),sum(-w))
lsth <- c(sum(w*term),sum(-w))
## second deriv...
term <- nu2^2 * trigamma(nu12)/4 + nu2 * digamma(nu12)/2 -
nu2^2 * trigamma(nu/2)/4 - nu2 * digamma(nu/2)/2 + 0.5*(nu2nu)^2 - 0.5*nu2nu
......@@ -1457,7 +1497,10 @@ scat <- function (theta = NULL, link = "identity",min.df = 3) {
preinitialize <- function(y,family) {
## initialize theta from raw observations..
if (family$n.theta>0) {
Theta <- c(-1, log(0.2*var(y)^.5))
## low df and low variance promotes indefiniteness.
## Better to start with moderate df and fairly high
## variance...
Theta <- c(1.5, log(0.8*sd(y)))
return(list(Theta=Theta))
} ## otherwise fixed theta supplied
}
......@@ -1468,10 +1511,12 @@ scat <- function (theta = NULL, link = "identity",min.df = 3) {
mustart <- y + (y == 0)*.1
})
postproc <- function(family,...) {
postproc <- function(family,y,prior.weights,fitted,linear.predictors,offset,intercept) {
posr <- list()
posr$family <-
paste("Scaled t(",paste(round(family$getTheta(TRUE),3),collapse=","),")",sep="")
posr$null.deviance <- find.null.dev(family,y,eta=linear.predictors,offset,prior.weights)
th <- round(family$getTheta(TRUE),3)
if (th[1]>999) th[1] <- Inf
posr$family <- paste("Scaled t(",paste(th,collapse=","),")",sep="")
posr
}
......@@ -1683,7 +1728,9 @@ ziP <- function (theta = NULL, link = "identity",b=0) {
ls <- function(y,w,theta,scale) {
## the log saturated likelihood function.
## ls is defined as zero for REML/ML expression as deviance is defined as -2*log.lik
## ls is defined as zero for REML/ML expression as deviance is defined as -2*log.lik
#vec <- !is.null(attr(theta,"vec.grad"))
#lsth1 <- if (vec) matrix(0,length(y),2) else c(0,0)
list(ls=0,## saturated log likelihood
lsth1=c(0,0), ## first deriv vector w.r.t theta - last element relates to scale
lsth2=matrix(0,2,2)) ##Hessian w.r.t. theta
......
......@@ -51,11 +51,13 @@ Sl.setup <- function(G) {
Sl[[b]]$start <- off[ind[1]]
Sl[[b]]$stop <- Sl[[b]]$start + nr - 1
Sl[[b]]$lambda <- rep(1,length(ind)) ## dummy at this stage
Sl[[b]]$repara <- FALSE
} else { ## singleton
Sl[[b]] <- list(start=off[ind], stop=off[ind]+nrow(G$S[[ind]])-1,
rank=G$rank[ind],S=list(G$S[[ind]]))
Sl[[b]]$S <- list(G$S[[ind]])
Sl[[b]]$lambda <- 1 ## dummy at this stage
Sl[[b]]$repara <- TRUE
} ## finished singleton
b <- b + 1
} ## finished this block
......@@ -211,7 +213,7 @@ Sl.setup <- function(G) {
St <- Sl[[b]]$S[[1]]/S.norm
lambda <- c(lambda,1/S.norm)
for (j in 2:length(Sl[[b]]$S)) {
S.norm <- norm(Sl[[b]]$S[[1]])
S.norm <- norm(Sl[[b]]$S[[j]])
St <- St + Sl[[b]]$S[[j]]/S.norm
lambda <- c(lambda,1/S.norm)
}
......@@ -234,7 +236,7 @@ Sl.Sb <- function(Sl,rho,beta) {
## multi-S blocks. Logic is identical to Sl.addS.
k <- 1
a <- beta * 0
for (b in 1:length(Sl)) {
if (length(Sl)>0) for (b in 1:length(Sl)) {
ind <- (Sl[[b]]$start:Sl[[b]]$stop)[Sl[[b]]$ind]
if (length(Sl[[b]]$S)==1) { ## singleton - multiple of identity
a[ind] <- a[ind] + beta[ind] * exp(rho[k])
......@@ -385,7 +387,7 @@ ldetS <- function(Sl,rho,fixed,np,root=FALSE,repara=TRUE,nt=1) {
d2.ldS <- matrix(0,n.deriv,n.deriv)
rp <- list() ## reparameterization list
if (root) E <- matrix(0,np,np) else E <- NULL
for (b in 1:length(Sl)) { ## work through blocks
if (length(Sl)>0) for (b in 1:length(Sl)) { ## work through blocks
if (length(Sl[[b]]$S)==1) { ## singleton
ldS <- ldS + rho[k.sp] * Sl[[b]]$rank
if (!fixed[k.sp]) {
......@@ -470,7 +472,7 @@ Sl.addS <- function(Sl,A,rho) {
## and should have already been applied to A using Sl.initial.repara
k <- 1
A <- A*1 ## force a copy to be made so that A not modified in calling env!!
for (b in 1:length(Sl)) {
if (length(Sl)>0) for (b in 1:length(Sl)) {
ind <- (Sl[[b]]$start:Sl[[b]]$stop)[Sl[[b]]$ind]
if (length(Sl[[b]]$S)==1) { ## singleton
B <- exp(rho[k]);diag <- -1
......@@ -543,7 +545,7 @@ Sl.mult <- function(Sl,A,k = 0,full=TRUE) {
## If k>0 then the routine forms S_k%*%A, zero padded
## if full==TRUE, but in smallest number of rows form otherwise.
nb <- length(Sl) ## number of blocks
if (nb==0) return(A*0)
Amat <- is.matrix(A)
if (k<=0) { ## apply whole penalty
B <- A*0
......@@ -635,7 +637,7 @@ Sl.termMult <- function(Sl,A,full=FALSE,nt=1) {
SA <- list()
k <- 0 ## component counter
nb <- length(Sl) ## number of blocks
for (b in 1:nb) { ## block loop
if (nb>0) for (b in 1:nb) { ## block loop
if (length(Sl[[b]]$S)==1) { ## singleton
k <- k + 1
ind <- Sl[[b]]$start:Sl[[b]]$stop
......@@ -811,7 +813,8 @@ Sl.iftChol <- function(Sl,XX,R,d,beta,piv,nt=1) {
} ## end Sl.iftChol
Sl.fitChol <- function(Sl,XX,f,rho,yy=0,L=NULL,rho0=0,log.phi=0,phi.fixed=TRUE,nobs=0,Mp=0,nt=1,tol=0) {
Sl.fitChol <- function(Sl,XX,f,rho,yy=0,L=NULL,rho0=0,log.phi=0,phi.fixed=TRUE,
nobs=0,Mp=0,nt=1,tol=0,gamma=1) {
## given X'WX in XX and f=X'Wy solves the penalized least squares problem
## with penalty defined by Sl and rho, and evaluates a REML Newton step, the REML
## gradiant and the the estimated coefs bhat. If phi.fixed=FALSE then we need
......@@ -855,16 +858,16 @@ Sl.fitChol <- function(Sl,XX,f,rho,yy=0,L=NULL,rho0=0,log.phi=0,phi.fixed=TRUE,n
phi <- exp(log.phi)
reml1 <- (dXXS$d1[!fixed] - ldS$ldet1 +
(dift$rss1[!fixed] + dift$bSb1[!fixed])/phi)/2
(dift$rss1[!fixed] + dift$bSb1[!fixed])/(phi*gamma))/2
reml2 <- (dXXS$d2[!fixed,!fixed] - ldS$ldet2 +
(dift$rss2[!fixed,!fixed] + dift$bSb2[!fixed,!fixed])/phi)/2
(dift$rss2[!fixed,!fixed] + dift$bSb2[!fixed,!fixed])/(phi*gamma))/2
if (!phi.fixed) {
n <- length(reml1)
rss.bSb <- yy - sum(beta*f) ## use identity ||y-Xb|| + b'Sb = y'y - b'X'y (b is minimizer)
reml1[n+1] <- (-rss.bSb/phi + nobs - Mp)/2
d <- c(-(dift$rss1[!fixed] + dift$bSb1[!fixed]),rss.bSb)/(2*phi)
reml1[n+1] <- (-rss.bSb/(phi*gamma) + nobs/gamma - Mp)/2
d <- c(-(dift$rss1[!fixed] + dift$bSb1[!fixed]),rss.bSb)/(2*phi*gamma)
reml2 <- rbind(cbind(reml2,d[1:n]),d)
if (!is.null(L)) L <- rbind(cbind(L,rep(0,nrow(L))),c(rep(0,ncol(L)),1))
}
......@@ -898,7 +901,7 @@ Sl.fitChol <- function(Sl,XX,f,rho,yy=0,L=NULL,rho0=0,log.phi=0,phi.fixed=TRUE,n
hess=hess,ldetS=ldS$ldetS,ldetXXS=ldetXXS)
} ## Sl.fitChol
Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NULL,Mp=0,nt=1) {
Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NULL,Mp=0,nt=1,gamma=1) {
## fits penalized regression model with model matrix X and
## initialised block diagonal penalty Sl to data in y, given
## log smoothing parameters rho.
......@@ -930,19 +933,19 @@ Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NU
dXXS <- d.detXXS(ldS$Sl,PP,nt=nt) ## derivs of log|X'X+S|
## all ingredients are now in place to form REML score and
## its derivatives....
reml <- (rss.bSb/phi + (nobs-Mp)*log(2*pi*phi) +
reml <- (rss.bSb/(phi*gamma) + (nobs/gamma-Mp)*log(2*pi*phi) + Mp*log(gamma) +
ldetXXS - ldS$ldetS)/2
reml1 <- (dXXS$d1[!fixed] - ldS$ldet1 + # dift$bSb1[!fixed]/phi)/2
(dift$rss1[!fixed] + dift$bSb1[!fixed])/phi)/2
(dift$rss1[!fixed] + dift$bSb1[!fixed])/(phi*gamma))/2
reml2 <- (dXXS$d2[!fixed,!fixed] - ldS$ldet2 + #dift$bSb2[!fixed,!fixed]/phi)/2
(dift$rss2[!fixed,!fixed] + dift$bSb2[!fixed,!fixed])/phi)/2
(dift$rss2[!fixed,!fixed] + dift$bSb2[!fixed,!fixed])/(phi*gamma))/2
## finally add in derivatives w.r.t. log.phi
if (!phi.fixed) {
n <- length(reml1)
reml1[n+1] <- (-rss.bSb/phi + nobs - Mp)/2
reml1[n+1] <- (-rss.bSb/(phi*gamma) + nobs/gamma - Mp)/2
#d <- c(-(dift$bSb1[!fixed]),rss.bSb)/(2*phi)
d <- c(-(dift$rss1[!fixed] + dift$bSb1[!fixed]),rss.bSb)/(2*phi)
d <- c(-(dift$rss1[!fixed] + dift$bSb1[!fixed]),rss.bSb)/(2*phi*gamma)
reml2 <- rbind(cbind(reml2,d[1:n]),d)
}
## following are de-bugging lines for testing derivatives of components...
......@@ -955,7 +958,7 @@ Sl.fit <- function(Sl,X,y,rho,fixed,log.phi=0,phi.fixed=TRUE,rss.extra=0,nobs=NU
} ## Sl.fit
fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
rss.extra=0,nobs=NULL,Mp=0,conv.tol=.Machine$double.eps^.5,nt=1) {
rss.extra=0,nobs=NULL,Mp=0,conv.tol=.Machine$double.eps^.5,nt=1,gamma=gamma) {
## estimates log smoothing parameters rho, by optimizing fast REML
## using Newton's method. On input Sl is a block diagonal penalty
## structure produced by Sl.setup, while X is a model matrix
......@@ -988,7 +991,7 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
fixed <- rep(FALSE,nrow(L))
best <- Sl.fit(Sl,X,y,L%*%rho+rho.0,fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt)
best <- Sl.fit(Sl,X,y,L%*%rho+rho.0,fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt,gamma=gamma)
## get a typical scale for the reml score...
reml.scale <- abs(best$reml) + best$rss/best$nobs
......@@ -1038,7 +1041,7 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
step[uconv.ind] <- uc.step ## step includes converged
## try out the step...
rho1 <- L%*%(rho + step)+rho.0; if (!phi.fixed) log.phi <- rho1[nr+1]
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt)
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt,gamma=gamma)
k <- 0
not.moved <- 0 ## count number of consecutive steps of essentially no change from best
while (trial$reml>best$reml) { ## step half until improvement or failure
......@@ -1052,7 +1055,7 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
}
step <- step/2;k <- k + 1
rho1 <- L%*%(rho + step)+rho.0; if (!phi.fixed) log.phi <- rho1[nr+1]
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt)
trial <- Sl.fit(Sl,X,y,rho1[1:nr],fixed,log.phi,phi.fixed,rss.extra,nobs,Mp,nt=nt,gamma=gamma)
}
if (step.failed) break ## can get no further
#if ((k==35 && trial$reml>best$reml)||(sum(rho != rho + step)==0)) { ## step has failed
......
......@@ -106,12 +106,12 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
if (inherits(family,"extended.family")) { ## then actually gam.fit4/5 is needed
if (inherits(family,"general.family")) {
return(gam.fit5(x,y,sp,Sl=Sl,weights=weights,offset=offset,deriv=deriv,
family=family,control=control,Mp=Mp,start=start))
family=family,control=control,Mp=Mp,start=start,gamma=gamma))
} else
return(gam.fit4(x, y, sp, Eb,UrS=UrS,
weights = weights, start = start, etastart = etastart,
mustart = mustart, offset = offset,U1=U1, Mp=Mp, family = family,
control = control, deriv=deriv,
control = control, deriv=deriv,gamma=gamma,
scale=scale,scoreType=scoreType,null.coef=null.coef,...))
}
......@@ -610,20 +610,21 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
ls <- family$ls(y,weights,n,scale)*n.true/nobs ## saturated likelihood and derivatives
Dp <- dev + oo$conv.tol + dev.extra
REML <- Dp/(2*scale) - ls[1] + oo$rank.tol/2 - rp$det/2 - remlInd*Mp/2*log(2*pi*scale)
REML <- (Dp/(2*scale) - ls[1])/gamma + oo$rank.tol/2 - rp$det/2 -
remlInd*(Mp/2*(log(2*pi*scale)-log(gamma)))
attr(REML,"Dp") <- Dp/(2*scale)
if (deriv) {
REML1 <- oo$D1/(2*scale) + oo$trA1/2 - rp$det1/2
if (deriv==2) REML2 <- (matrix(oo$D2,nSp,nSp)/scale + matrix(oo$trA2,nSp,nSp) - rp$det2)/2
REML1 <- oo$D1/(2*scale*gamma) + oo$trA1/2 - rp$det1/2
if (deriv==2) REML2 <- (matrix(oo$D2,nSp,nSp)/(scale*gamma) + matrix(oo$trA2,nSp,nSp) - rp$det2)/2
if (sum(!is.finite(REML2))) {
stop("Non finite derivatives. Try decreasing fit tolerance! See `epsilon' in `gam.contol'")
}
}
if (!scale.known&&deriv) { ## need derivatives wrt log scale, too
##ls <- family$ls(y,weights,n,scale) ## saturated likelihood and derivatives
dlr.dlphi <- -Dp/(2 *scale) - ls[2]*scale - Mp/2*remlInd
d2lr.d2lphi <- Dp/(2*scale) - ls[3]*scale^2 - ls[2]*scale
d2lr.dspphi <- -oo$D1/(2*scale)
dlr.dlphi <- (-Dp/(2 *scale) - ls[2]*scale)/gamma - Mp/2*remlInd
d2lr.d2lphi <- (Dp/(2*scale) - ls[3]*scale^2 - ls[2]*scale)/gamma
d2lr.dspphi <- -oo$D1/(2*scale*gamma)
REML1 <- c(REML1,dlr.dlphi)
if (deriv==2) {
REML2 <- rbind(REML2,as.numeric(d2lr.dspphi))
......@@ -643,7 +644,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
K <- oo$rank.tol/2 - rp$det/2
REML <- Dp/(2*phi) - ls[1] + K - Mp/2*log(2*pi*phi)*remlInd
REML <- (Dp/(2*phi) - ls[1]) + K - Mp/2*(log(2*pi*phi))*remlInd
attr(REML,"Dp") <- Dp/(2*phi)
if (deriv) {
phi1 <- oo$P1; Dp1 <- oo$D1; K1 <- oo$trA1/2 - rp$det1/2;
......@@ -1446,12 +1447,11 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
delta1 <- delta + step
lsp1 <- rt(delta1,lsp1.max)$rho ## transform to log sp space
} else lsp1 <- lsp + step
b1<-gam.fit3(x=X, y=y, sp=L%*%lsp1+lsp0,Eb=Eb,UrS=UrS,
offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=0,
control=control,gamma=gamma,scale=scale,
printWarn=FALSE,start=start,mustart=mustart,scoreType=scoreType,
null.coef=null.coef,pearson.extra=pearson.extra,
dev.extra=dev.extra,n.true=n.true,Sl=Sl,...)
b1<-gam.fit3(x=X, y=y, sp=L%*%lsp1+lsp0,Eb=Eb,UrS=UrS,offset = offset,U1=U1,Mp=Mp,
family = family,weights=weights,deriv=0,control=control,gamma=gamma,
scale=scale,printWarn=FALSE,start=start,mustart=mustart,scoreType=scoreType,
null.coef=null.coef,pearson.extra=pearson.extra,dev.extra=dev.extra,
n.true=n.true,Sl=Sl,...)
pred.change <- sum(grad*step) + 0.5*t(step) %*% hess %*% step ## Taylor prediction of change
if (reml) {
score1 <- b1$REML
......@@ -1461,7 +1461,9 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
score1 <- b1$UBRE
} else score1 <- b1$GCV