Commit 0159f649 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Debian changes 1.8-8-1

mgcv (1.8-8-1) unstable; urgency=low

  * New upstream release
  
  * debian/control: Set Build-Depends: to current R version 
parents 8909ce92 464de942
** denotes quite substantial/important changes
*** denotes really big changes
1.8-8
** New "gp" smooth class (see ?gp.smooth) implemeting the Matern
covariance based Gaussian process model of Kamman and Wand (2003), and a
variety of other simple GP smoothers.
* some smooth plot methods now accept 'colors' and 'contour.col' argument to
set color palette in image plots and contour line colors.
* predict.gam and predict.bam now accept an 'exclude' argument allowing
terms (e.g. random effects) to be zeroed for prediction. For efficiency,
smooth terms not in 'terms' or in 'exclude' are no longer evaluated, and
are instead set to zero or not returned. See ?predict.gam.
* ocat saturated likelihood definition changed to zero, leading to better
comprability of deviance between model fits (thanks to Herwig Friedl).
* null.deviance calculation for extended families modified to make more sense
when `mu' is the mean of a latent variable, rather than the response itself.
* bam now returns standarized residuals 'std.rsd' if `rho!=0'.
* bam(...,discrete=TRUE) can now handle 'fs' terms.
* bam(...,discrete=TRUE) now accepts 'by' variables. Thanks to Zheyaun Li
for debugging on this.
* bam now works with drop.unused.levels == TRUE when random effects should
have more levels than those that exist in data. (Thanks Alec Leh)
* bam chunk.size logic error fix - error could be triggered if chunk.size
reset automaticlly to be larger than data size.
* uniqucombs can now accept a data frame with some or all factor columns,
as well as purely numeric marices.
* discrete.mf modified to avoid discretizing a covariate more than once,
and to halt if a model requires the same covariate to be discretized
two different ways (e.g. s(x) + s(x,z)). This affects only
bam(...,discrete=TRUE).
* Some changes to ziP and ziplss families to improve numerical robustness,
and to ziP help file to suggest appropriate checking. Thanks to Keren
Raiter, for reporting problems.
* numerical robustness of extended gam methods (gam.fit4) improved for cases
with many zero or near zero iterative weights. Handling of zero weights
modified to avoid divide by (near) zero problems. Also tests for poor
scaling of sqrt(abs(w))*z and substitutes computations based on w*z if
detected. Also 'newton' routine now step halves if REML score not finite!
* Sl.setup (used by bam) modification to allow more efficient handling of terms
with multiple diagonal penalties with no non-zero elements in common, but
possibly with non zero elements `interleaved' between penalties.
1.8-7
** 'gam' default scale parameter changed to modified Pearson estimator
......
Package: mgcv
Version: 1.8-7
Version: 1.8-8
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: 2015-07-22 17:38:22 UTC; sw283
Packaged: 2015-10-24 06:15:25 UTC; sw283
Repository: CRAN
Date/Publication: 2015-07-23 07:05:08
Date/Publication: 2015-10-24 09:17:33
f876a44fccfab3bcf4e8c6a8916f480f *ChangeLog
68b8a0864224d96ff69d755507862f24 *DESCRIPTION
bcb7f62d1e1297b1d685a095b00220aa *ChangeLog
f133d9c24128ba8a977d1674d59f2a23 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
053c3c177572e1570e63329445eacc14 *NAMESPACE
e99642616897453fa1cab582bc73c84f *R/bam.r
b1b86938209b64b673e096da331769eb *NAMESPACE
9cf2059cc5ab5fd2f8a9f4439fec53e5 *R/bam.r
52ea7081e235df890d3e51198bcc3092 *R/coxph.r
209a7a1d10cec589213ffed9efc1d5cc *R/efam.r
9131483dd63ef2b80e8159601cc03922 *R/fast-REML.r
88dbdd15671f0a924bdd9f5f42a89a34 *R/gam.fit3.r
a6208b52c8af1f8779dd0b4b482c3fa7 *R/gam.fit4.r
09790e96aac078861d04d160d852c6d7 *R/efam.r
5a542999fd494e3872c877c6e1dc2706 *R/fast-REML.r
57f8b12b97304227870a4e43a33b5d75 *R/gam.fit3.r
472d6e64f3f403aad4f6550f921a4bba *R/gam.fit4.r
e63276b78a4ff2566af2e651bfc6aa9c *R/gam.sim.r
633c048e36083646c7d409e99cf7b037 *R/gamlss.r
29f3e7f145dabb21b25b6942066fab3a *R/gamlss.r
ff163969e9ad7a38857907bf38a39ec0 *R/gamm.r
3b0d5cac4a59ef1a8cb325b249699476 *R/jagam.r
d5ce48c9ea98c2f428cfb00bcc9fb5ea *R/mgcv.r
479a34c52f6ec8a04e7c66cfcf4bdf58 *R/mgcv.r
59a92561754c8156506af4004f594ac9 *R/misc.r
66c24aed2f8fc9f6bce321794f8aff87 *R/mvam.r
f031621e351c06b12c328e4033e2f097 *R/plots.r
966c267722f1027c53aa5621614333ae *R/smooth.r
68348617c5d3b0e2ea805c4360a8cdd4 *R/soap.r
76192b5b36d8828af6e50013823732b4 *R/plots.r
6069c2b79e5678d8c6ff35b65c6d180c *R/smooth.r
9d47fd13872418b318323699de89dfb4 *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
......@@ -39,7 +39,7 @@ e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
f12468625253dd3603907de233762fd6 *man/Rrank.Rd
f6cadda5999c35800fd65c08d6812f7b *man/Tweedie.Rd
80f8763baa4987579e2aa56073a9e94e *man/anova.gam.Rd
e49e38bc1741829ae997f8b0b4fcd832 *man/bam.Rd
b57775363a39cd704e37f3b3d06962f6 *man/bam.Rd
71bde8b8caa24a36529ce7e0ac3165d8 *man/bam.update.Rd
a2beb811b1093c5e82ef32d7de1f7d32 *man/cSplineDes.Rd
9b4d616d1b6c4a46ca77d16cded3f806 *man/choose.k.Rd
......@@ -86,13 +86,13 @@ b440c31b77c368b484d9665b2b3e89fb *man/jagam.Rd
19a3f7da03e76d7d2fc0a4e544564c72 *man/magic.Rd
5169af4be5fccf9fa79b6de08e9ea035 *man/magic.post.proc.Rd
59672d8599211ce9cf44252f92110920 *man/mgcv-FAQ.Rd
01926b8c2bda6fc74b51f6ec7f5cc620 *man/mgcv-package.Rd
92499f4c6bc65a0cffd31bd8bc993a9f *man/mgcv-parallel.Rd
38c75ee7bb718c277da5bc649e6886e7 *man/mgcv-package.Rd
9742c9bd1d42aaebe27e2b4218d12672 *man/mgcv-parallel.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
bc9b89db7e7ff246749551c16f5f1f07 *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
d0c2b73407680aac50a17753346d814f *man/mvn.Rd
9e95f6b9eec476188de830beb9ff3e0b *man/negbin.Rd
61e6a0a7fa986b4b6278c2d3c44c9c71 *man/negbin.Rd
8b766a6ad848b0f1ca469e381ded0169 *man/new.name.Rd
7d8f62e182f7c428cc9d46ddd4d97d43 *man/notExp.Rd
445571cf0b115bf19a2d3474e70f59e7 *man/notExp2.Rd
......@@ -103,10 +103,10 @@ c64b44bc684cbadfe77e7f6dc833ddb4 *man/pcls.Rd
8c0f8575b427f30316b639a326193aeb *man/pdTens.Rd
1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd
edf57071572275a8443b2f0b66d44424 *man/place.knots.Rd
04a9af3e294bd14ad590aacb6d26332d *man/plot.gam.Rd
b903ebcf31703db156e033fdfa527d73 *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
085dff270e5255315cb244368faf79df *man/predict.bam.Rd
900647986ab1365d437368cc28b62658 *man/predict.gam.Rd
fe39dcd63a4a01bdac0604389d59bb41 *man/predict.bam.Rd
3d811dd86252174c0108e2b1b75aaf9a *man/predict.gam.Rd
7562cb5ce9a7fbf0cd60124cf4305315 *man/print.gam.Rd
b8ac3ed8fc05bb14c605b996404853cd *man/qq.gam.Rd
f77ca1471881d2f93c74864d076c0a0e *man/rTweedie.Rd
......@@ -116,40 +116,41 @@ f6f1333be2587ffef5970905b13468ea *man/rig.Rd
9f6f46f5c5da080bc82f9aa4685d364a *man/rmvn.Rd
7258dfc0149fff020054117fd2ee6bd8 *man/s.Rd
c438eb3e41cb1f51e72283380145d210 *man/scat.Rd
2199afd400a1821b74a24e2578bc0224 *man/single.index.Rd
9641bb5e2573b5fcbdc303fcf160f3e6 *man/single.index.Rd
6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
b5a06918956fd10f23285b453c05bdb4 *man/smooth.construct.Rd
2318b61ce2e4770d880abfeaba194577 *man/smooth.construct.Rd
4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.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
3a0545da6f80f00e0b7347f8083a10b3 *man/smooth.construct.gp.smooth.spec.Rd
4aaa84b520992fbc32b0c37f7f63c1dd *man/smooth.construct.mrf.smooth.spec.Rd
abe15377f471a2d8957a59c19eeef0bb *man/smooth.construct.ps.smooth.spec.Rd
6aaff8575f68775ed21930893ea9e03d *man/smooth.construct.re.smooth.spec.Rd
a8d5eb4772641567fee15714f01a4727 *man/smooth.construct.so.smooth.spec.Rd
9c35c35b9583281eaa55b97c9b9fe770 *man/smooth.construct.sos.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
75a042d74460dd0c6e0c0a95c1f5d3f7 *man/smooth.terms.Rd
8120512d52c741f8c27c2689f69b63d3 *man/smoothCon.Rd
d4083ff900aa69fa07610e0af2a2987b *man/smooth.terms.Rd
8fd3089f8fd1379736ae93a67f4febd6 *man/smoothCon.Rd
b55a396da77559dac553613146633f97 *man/sp.vcov.Rd
83bd8e097711bf5bd0fff09822743d43 *man/spasm.construct.Rd
b9394812e5398ec95787c65c1325a027 *man/step.gam.Rd
f0791d830687d6155efb8a73db787401 *man/summary.gam.Rd
47d22ea4fe93bc0fa0f47bd262979941 *man/t2.Rd
fd71e0a218ad840e6ce4b873f9fa083e *man/t2.Rd
63e2050f9a01b12096d640e87916b78a *man/te.Rd
6eebb6ef90374ee09453d6da6449ed79 *man/tensor.prod.model.matrix.Rd
71add05e5823dcc0509f435ab4dc4627 *man/uniquecombs.Rd
42542181aa314eed22f05907a8546735 *man/uniquecombs.Rd
a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
2c5f6815e609f2cdb56b0067a183f915 *man/ziP.Rd
1e514afeea034edfbd6ba563561a71de *man/ziP.Rd
ae8388103d8b1be39f55f426b205b576 *man/ziplss.Rd
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
be84cc1bdb81bb411322266c34a0bf1d *po/R-mgcv.pot
7b07899266c3acf3d2a625850d7cd6ef *po/R-pl.po
382c94188dbc193fca9628287b66d1af *po/de.po
93f72334356fe6f05a64e567efd35c8e *po/en@quot.po
......@@ -160,15 +161,15 @@ dfd4eec9edc7d1ab6354d47b6e2bd42f *po/pl.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
e71c1a1624b431fbab0a4c8f151d2a97 *src/coxph.c
48f0c9de0da60d48cf4b306c2fdd039a *src/discrete.c
e4236680abfb9c89c7bf09b441c755c2 *src/gdi.c
5a8e93ae5f4fae722c1304cecc9c94b3 *src/gdi.c
2436f9b328e80370ce2203dbf1dd813c *src/general.h
890105488631271ad8b184aa0524b59f *src/init.c
d95aa7d646833d1e56f2414f32637bb3 *src/init.c
9a5d7cb3cf93cbdfc08353fbd20a270e *src/magic.c
55b92c3ab6a037f9742542c2af8cad56 *src/mat.c
53cde45313485d2d15e6377461dd3f46 *src/mat.c
aae0f298e384952bb8b6b923d40520a8 *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
ec2ae157a9e3bedcccc88b053d0a4e0c *src/mgcv.c
39b22a8fa2128fd6b52b3601841b0393 *src/mgcv.h
d437a2cf1548e79c8dc910657fb1cabb *src/mgcv.h
6b90633f745d3b3314ec5a784bde99d0 *src/misc.c
08c1706ffeec4277c484435a0644b0e3 *src/mvn.c
cbe54250deb38aa5f88f8b669a4468cd *src/qp.c
......
......@@ -37,6 +37,7 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
Predict.matrix.pspline.smooth,
Predict.matrix.random.effect,
Predict.matrix.t2.smooth,
Predict.matrix.gp.smooth,
qq.gam,
residuals.gam,rig,rTweedie,rmvn,
Rrank,s,scat,
......@@ -54,6 +55,7 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
smooth.construct.ps.smooth.spec,
smooth.construct.re.smooth.spec,
smooth.construct.mrf.smooth.spec,
smooth.construct.gp.smooth.spec,
smooth.construct.sos.smooth.spec,
smooth.construct.so.smooth.spec,
smooth.construct.sf.smooth.spec,
......@@ -136,6 +138,7 @@ S3method(smooth.construct, cs.smooth.spec)
S3method(smooth.construct,ds.smooth.spec)
S3method(smooth.construct, fs.smooth.spec)
S3method(smooth.construct, mrf.smooth.spec)
S3method(smooth.construct, gp.smooth.spec)
S3method(smooth.construct,ps.smooth.spec)
S3method(smooth.construct, re.smooth.spec)
S3method(smooth.construct,so.smooth.spec)
......@@ -164,6 +167,7 @@ S3method(Predict.matrix,soap.film)
S3method(Predict.matrix,sf)
S3method(Predict.matrix,sw)
S3method(Predict.matrix,sos.smooth)
S3method(Predict.matrix,gp.smooth)
S3method(spasm.construct,cus)
......
This diff is collapsed.
This diff is collapsed.
......@@ -81,19 +81,27 @@ Sl.setup <- function(G) {
## easily here means no overlap in penalties
Sl[[b]]$S <- G$smooth[[i]]$S
nb <- nrow(Sl[[b]]$S[[1]])
sbStart <- sbStop <- rep(NA,m)
sbdiag <- sbStart <- sbStop <- rep(NA,m)
ut <- upper.tri(Sl[[b]]$S[[1]],diag=FALSE)
## overlap testing requires the block ranges
for (j in 1:m) { ## get block range for each S[[j]]
sbdiag[j] <- sum(abs(Sl[[b]]$S[[j]][ut]))==0 ## is penalty diagonal??
ir <- range((1:nb)[rowSums(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
for (j in 1:m) { ## test for overlap
itot <- rep(FALSE,nb)
for (k in 1:m) if (j!=k) itot[sbStart[k]:sbStop[k]] <- TRUE
if (sum(itot[sbStart[j]:sbStop[j]])>0) { ## no good, some overlap detected
split.ok <- FALSE; break
if (all(sbdiag)) { ## it's all diagonal - can allow interleaving
for (k in 1:m) if (j!=k) itot[diag(Sl[[b]]$S[[k]])!=0] <- TRUE
if (sum(itot[diag(Sl[[b]]$S[[j]])!=0])>0) { ## no good, some overlap detected
split.ok <- FALSE; break
}
} else { ## not diagonal - really need on overlapping blocks
for (k in 1:m) if (j!=k) itot[sbStart[k]:sbStop[k]] <- TRUE
if (sum(itot[sbStart[j]:sbStop[j]])>0) { ## no good, some overlap detected
split.ok <- FALSE; break
}
}
}
if (split.ok) { ## can split this block into m separate singleton blocks
......@@ -108,9 +116,8 @@ Sl.setup <- function(G) {
}
} else { ## not possible to split
Sl[[b]]$S <- G$smooth[[i]]$S
b <- b + 1 ## next block!!
} ## additive block finished
b <- b + 1 ## next block!!
} ## additive block finished
}
......
......@@ -386,20 +386,22 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
if (sum(good)<ncol(x)) stop("Not enough informative observations.")
if (control$trace) t1 <- proc.time()
oo <- .C(C_pls_fit1,y=as.double(z),X=as.double(x[good,]),w=as.double(w),
oo <- .C(C_pls_fit1,y=as.double(z),X=as.double(x[good,]),w=as.double(w),wy=as.double(w*z),
E=as.double(Sr),Es=as.double(Eb),n=as.integer(sum(good)),
q=as.integer(ncol(x)),rE=as.integer(rows.E),eta=as.double(z),
penalty=as.double(1),rank.tol=as.double(rank.tol),nt=as.integer(control$nthreads))
penalty=as.double(1),rank.tol=as.double(rank.tol),nt=as.integer(control$nthreads),
use.wy=as.integer(0))
if (control$trace) tc <- tc + sum((proc.time()-t1)[c(1,4)])
if (!fisher&&oo$n<0) { ## likelihood indefinite - switch to Fisher for this step
z <- (eta - offset)[good] + (yg - mug)/mevg
w <- (weg * mevg^2)/var.mug
if (control$trace) t1 <- proc.time()
oo <- .C(C_pls_fit1,y=as.double(z),X=as.double(x[good,]),w=as.double(w),
oo <- .C(C_pls_fit1,y=as.double(z),X=as.double(x[good,]),w=as.double(w),wy=as.double(w*z),
E=as.double(Sr),Es=as.double(Eb),n=as.integer(sum(good)),
q=as.integer(ncol(x)),rE=as.integer(rows.E),eta=as.double(z),
penalty=as.double(1),rank.tol=as.double(rank.tol),nt=as.integer(control$nthreads))
penalty=as.double(1),rank.tol=as.double(rank.tol),nt=as.integer(control$nthreads),
use.wy=as.integer(0))
if (control$trace) tc <- tc + sum((proc.time()-t1)[c(1,4)])
}
......@@ -1436,7 +1438,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
## accept if improvement, else step halve
ii <- 0 ## step halving counter
##sc.extra <- 1e-4*sum(grad*Nstep) ## -ve sufficient decrease
if (score1<score && pdef) { ## immediately accept step if it worked and positive definite
if (is.finite(score1) && score1<score && pdef) { ## immediately accept step if it worked and positive definite
old.score <- score
mustart <- b$fitted.values
etastart <- b$linear.predictors
......@@ -1460,10 +1462,10 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
grad <- rho$rho1*grad
}
} else if (score1>=score) { ## initial step failed to improve score, try step halving ...
} else if (!is.finite(score1) || score1>=score) { ## initial step failed to improve score, try step halving ...
step <- Nstep ## start with the (pseudo) Newton direction
##sc.extra <- 1e-4*sum(grad*step) ## -ve sufficient decrease
while (score1>=score && ii < maxHalf) {
while ((!is.finite(score1) || score1>=score) && ii < maxHalf) {
if (ii==3&&i<10) { ## Newton really not working - switch to SD, but keeping step length
s.length <- min(sum(step^2)^.5,maxSstep)
step <- Sstep*s.length/sum(Sstep^2)^.5 ## use steepest descent direction
......@@ -1489,7 +1491,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
score1 <- b1$UBRE
} else score1 <- b1$GCV
##sc.extra <- 1e-4*sum(grad*Nstep) ## -ve sufficient decrease
if (score1 < score) { ## accept
if (is.finite(score1) && score1 < score) { ## accept
if (pdef||!sd.unused) { ## then accept and compute derivatives
b <- 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=2,
......@@ -1523,7 +1525,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
}
score1 <- score - abs(score) - 1 ## make sure that score1 < score
} # end of if (score1<= score ) # accept
if (score1>=score) ii <- ii + 1
if (!is.finite(score1) || score1>=score) ii <- ii + 1
} ## end while (score1>score && ii < maxHalf)
if (!pdef&&sd.unused&&ii<maxHalf) score1 <- score2
} ## end of step halving branch
......@@ -1535,7 +1537,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
if (!pdef&&sd.unused) {
step <- Sstep*2
kk <- 0
kk <- 0;score2 <- NA
ok <- TRUE
while (ok) { ## step length loop for steepest....
step <- step/2;kk <- kk+1
......@@ -1556,18 +1558,18 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
} else if (scoreType=="UBRE") {
score3 <- b1$UBRE
} else score3 <- b1$GCV
if (kk==1||score3<=score2) { ## accept step - better than last try
if (!is.finite(score2)||(is.finite(score3)&&score3<=score2)) { ## accept step - better than last try
score2 <- score3
lsp2 <- lsp3
## 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
if ((is.finite(score2)&&is.finite(score3)&&score2<score&&score3>score2)||kk==40) ok <- FALSE
} ## while (ok) ## step length control loop
## now pick the step that led to the biggest decrease
if (score2<score1) {
if (is.finite(score2) && score2<score1) {
lsp1 <- lsp2
if (!is.null(lsp.max)) delta1 <- delta
score1 <- score2
......@@ -1621,9 +1623,13 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
if (ii==maxHalf) converged <- TRUE ## step failure
if (converged) break
} ## end of iteration loop
if (ii==maxHalf) ct <- "step failed"
else if (i==200) ct <- "iteration limit reached"
else ct <- "full convergence"
if (ii==maxHalf) {
ct <- "step failed"
warning("Fitting terminated with step failure - check results carefully")
} else if (i==200) {
ct <- "iteration limit reached"
warning("Iteration limit reached without full convergence - check carefully")
} else ct <- "full convergence"
list(score=score,lsp=lsp,lsp.full=L%*%lsp+lsp0,grad=grad,hess=hess,iter=i,conv =ct,score.hist = score.hist[!is.na(score.hist)],object=b)
} ## newton
......
This diff is collapsed.
......@@ -730,6 +730,7 @@ ldg <- function(g,deriv=4) {
## so l'' = alpha*(b-eg)...
b <- egi*(1+egi/6)/2
l2[ind] <- a[ind]*(b-egi)
l2[ii] <- -exp(g[ii])
l3 <- l4 <- NULL
## in a similar vein l3 can be robustified...
if (deriv>1) {
......@@ -745,7 +746,13 @@ ldg <- function(g,deriv=4) {
l4[ii] <- -exp(g[ii])
}
list(l1=-a,l2=l2,l3=l3,l4=l4)
l1=-a
ghi <- log(.Machine$double.xmax)/5
ii <- g > ghi
if (sum(ii)) {
l1[ii] <- l2[ii] <- l3[ii] <- l4[ii] <- -exp(ghi)
}
list(l1=l1,l2=l2,l3=l3,l4=l4)
} ## ldg
lde <- function(eta,deriv=4) {
......@@ -832,8 +839,8 @@ zipll <- function(y,g,eta,deriv=0) {
## order gggg,ggge,ggee,geee,eeee
l4 <- matrix(0,n,5)
l4[!zind,1] <- lg$l4[!zind] ## l_gggg, y>0
l4[!zind,4] <- le$l4[!zind] ## l_eeee, y>0
l4[zind,4] <- l[zind] ## l_eeee, y=0
l4[!zind,5] <- le$l4[!zind] ## l_eeee, y>0
l4[zind,5] <- l[zind] ## l_eeee, y=0
}
list(l=l,l1=l1,l2=l2,l3=l3,l4=l4,El2=El2)
} ## zipll
......
......@@ -372,7 +372,7 @@ augment.smX <- function(sm,nobs,np) {
## identifiability constraint purposes.
ns <- length(sm$S) ## number of penalty matrices
if (ns==0) { ## nothing to do
return(rbind(sm$X),matrix(0,np,np))
return(rbind(sm$X,matrix(0,np,ncol(sm$X))))
}
ind <- colMeans(abs(sm$S[[1]]))!=0
sqrmaX <- mean(abs(sm$X[,ind]))^2
......@@ -503,10 +503,10 @@ gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5,with.pen=FALSE)
}
} ## penalty matrices finished
## Now we need to establish null space rank for the term
m <- length(sm[[i]]$S)
if (m>0) {
mi <- length(sm[[i]]$S)
if (mi>0) {
St <- sm[[i]]$S[[1]]/norm(sm[[i]]$S[[1]],type="F")
if (m>1) for (j in 1:m) St <- St +
if (mi>1) for (j in 1:mi) St <- St +
sm[[i]]$S[[j]]/norm(sm[[i]]$S[[j]],type="F")
es <- eigen(St,symmetric=TRUE,only.values=TRUE)
sm[[i]]$null.space.dim <- sum(es$values<max(es$values)*.Machine$double.eps^.75)
......@@ -872,7 +872,7 @@ gam.setup <- function(formula,pterms,
data=stop("No data supplied to gam.setup"),knots=NULL,sp=NULL,
min.sp=NULL,H=NULL,absorb.cons=TRUE,sparse.cons=0,select=FALSE,idLinksBases=TRUE,
scale.penalty=TRUE,paraPen=NULL,gamm.call=FALSE,drop.intercept=FALSE,
diagonal.penalty=FALSE)
diagonal.penalty=FALSE,apply.by=TRUE)
## set up the model matrix, penalty matrices and auxilliary information about the smoothing bases
## needed for a gam fit.
## elements of returned object:
......@@ -1007,12 +1007,14 @@ gam.setup <- function(formula,pterms,
id <- split$smooth.spec[[i]]$id
if (is.null(id)||!idLinksBases) { ## regular evaluation
sml <- smoothCon(split$smooth.spec[[i]],data,knots,absorb.cons,scale.penalty=scale.penalty,
null.space.penalty=select,sparse.cons=sparse.cons,diagonal.penalty=diagonal.penalty)
null.space.penalty=select,sparse.cons=sparse.cons,
diagonal.penalty=diagonal.penalty,apply.by=apply.by)
} else { ## it's a smooth with an id, so basis setup data differs from model matrix data
names(id.list[[id]]$data) <- split$smooth.spec[[i]]$term ## give basis data suitable names
sml <- smoothCon(split$smooth.spec[[i]],id.list[[id]]$data,knots,
absorb.cons,n=nrow(data),dataX=data,scale.penalty=scale.penalty,
null.space.penalty=select,sparse.cons=sparse.cons,diagonal.penalty=diagonal.penalty)
null.space.penalty=select,sparse.cons=sparse.cons,
diagonal.penalty=diagonal.penalty,apply.by=apply.by)
}
for (j in 1:length(sml)) {
newm <- newm + 1
......@@ -1024,7 +1026,15 @@ gam.setup <- function(formula,pterms,
## at this stage, it is neccessary to impose any side conditions required
## for identifiability
if (m>0) sm <- gam.side(sm,X,tol=.Machine$double.eps^.5)
if (m>0) {
sm <- gam.side(sm,X,tol=.Machine$double.eps^.5)
if (!apply.by) for (i in 1:length(sm)) { ## restore any by-free model matrices
if (!is.null(sm[[i]]$X0)) { ## there is a by-free matrix to restore
ind <- attr(sm[[i]],"del.index") ## columns, if any to delete
sm[[i]]$X <- if (is.null(ind)) sm[[i]]$X0 else sm[[i]]$X0[,-ind,drop=FALSE]
}
}
}
## The matrix, L, mapping the underlying log smoothing parameters to the
## log of the smoothing parameter multiplying the S[[i]] must be
......@@ -2501,7 +2511,7 @@ model.matrix.gam <- function(object,...)
predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclude=NULL,
block.size=NULL,newdata.guaranteed=FALSE,na.action=na.pass,
unconditional=FALSE,...) {
......@@ -2771,11 +2781,14 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
if (!is.null(drop.ind)) X <- X[,-drop.ind]
if (n.smooth) for (k in 1:n.smooth) { ## loop through smooths
Xfrag <- PredictMat(object$smooth[[k]],data)
X[,object$smooth[[k]]$first.para:object$smooth[[k]]$last.para] <- Xfrag
Xfrag.off <- attr(Xfrag,"offset") ## any term specific offsets?
if (!is.null(Xfrag.off)) { Xoff[,k] <- Xfrag.off; any.soff <- TRUE }
if (type=="terms"||type=="iterms") ColNames[n.pterms+k] <- object$smooth[[k]]$label
klab <- object$smooth[[k]]$label
if ((is.null(terms)||(klab%in%terms))&&(is.null(exclude)||!(klab%in%exclude))) {
Xfrag <- PredictMat(object$smooth[[k]],data)
X[,object$smooth[[k]]$first.para:object$smooth[[k]]$last.para] <- Xfrag
Xfrag.off <- attr(Xfrag,"offset") ## any term specific offsets?
if (!is.null(Xfrag.off)) { Xoff[,k] <- Xfrag.off; any.soff <- TRUE }
}
if (type=="terms"||type=="iterms") ColNames[n.pterms+k] <- klab
} ## smooths done
......@@ -2913,13 +2926,27 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
rm(X)
} ## end of prediction block loop
if ((type=="terms"||type=="iterms")&&!is.null(terms)) { # return only terms requested via `terms'
if (sum(!(terms %in%colnames(fit))))
warning("non-existent terms requested - ignoring")
else {
fit <- fit[,terms,drop=FALSE]
if (se.fit) {
se <- se[,terms,drop=FALSE]
if ((type=="terms"||type=="iterms")&&(!is.null(terms)||!is.null(exclude))) { # return only terms requested via `terms'
cnames <- colnames(fit)
if (!is.null(terms)) {
if (sum(!(terms %in%cnames)))
warning("non-existent terms requested - ignoring")
else {
fit <- fit[,terms,drop=FALSE]
if (se.fit) {
se <- se[,terms,drop=FALSE]
}
}
}
if (!is.null(exclude)) {
if (sum(!(exclude %in%cnames)))
warning("non-existent exclude terms requested - ignoring")
else {
exclude <- which(cnames%in%exclude) ## convert to numeric column index
fit <- fit[,-exclude,drop=FALSE]
if (se.fit) {
se <- se[,-exclude,drop=FALSE]
}
}
}
}
......@@ -3639,8 +3666,9 @@ summary.gam <- function (object, dispersion = NULL, freq = FALSE, p.type=0, ...)
chi.sq[i] <- t(p)%*%V%*%p
df[i] <- attr(V, "rank")
} else { ## Better founded alternatives...
Xt <- X[,start:stop,drop=FALSE]
if (object$smooth[[i]]$null.space.dim==0&&!is.null(object$R)) { ## random effect or fully penalized term
Xt <- X[,start:stop,drop=FALSE]
fx <- if (inherits(object$smooth[[i]],"tensor.smooth")) all(object$smooth[[i]]$fx) else object$smooth[[i]]$fixed
if (!fx&&object$smooth[[i]]$null.space.dim==0&&!is.null(object$R)) { ## random effect or fully penalized term
res <- reTest(object,i)
} else { ## Inverted Nychka interval statistics
df[i] <- min(ncol(Xt),edf1[i])
......
......@@ -438,13 +438,15 @@ plot.sos.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,
pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,...) {
shift=0,trans=I,by.resids=FALSE,scheme=0,colors=heat.colors(100),
contour.col=3,...) {
## plot method function for sos.smooth terms
if (scheme>1) return(plot.mgcv.smooth(x,P=P,data=data,label=label,se1.mult=se1.mult,se2.mult=se2.mult,
partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,
pers=pers,theta=theta,phi=phi,jit=jit,xlab=xlab,ylab=ylab,main=main,
ylim=ylim,xlim=xlim,too.far=too.far,shade=shade,shade.col=shade.col,
shift=shift,trans=trans,by.resids=by.resids,scheme=scheme-2,...))
shift=shift,trans=trans,by.resids=by.resids,scheme=scheme-2,
colors=colors,contour.col=contour.col,...))
## convert location of pole in plotting grid to radians
phi <- phi*pi/180
theta <- theta*pi/180
......@@ -504,7 +506,7 @@ plot.sos.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (scheme == 0) {
col <- 1# "lightgrey
zz[P$ind] <- P$fit
image(P$xm,P$ym,matrix(zz,m,m),col=heat.colors(100),axes=FALSE,xlab="",ylab="",...)
image(P$xm,P$ym,matrix(zz,m,m),col=colors,axes=FALSE,xlab="",ylab="",...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
......@@ -514,7 +516,7 @@ plot.sos.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
zz[P$ind] <- P$lo
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,levels=c(-8:9*20),col=col,...)
zz[P$ind] <- P$fit
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,col=3,...)
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,col=contour.col,...)
} else if (scheme == 1) {
col <- 1
zz[P$ind] <- P$fit
......@@ -694,7 +696,8 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,
pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,...) {
shift=0,trans=I,by.resids=FALSE,scheme=0,colors=heat.colors(50),
contour.col=3,...) {
## default plot method for smooth objects `x' inheriting from "mgcv.smooth"
## `x' is a smooth object, usually part of a `gam' fit. It has an attribute
## 'coefficients' containg the coefs for the smooth, but usually these
......@@ -908,7 +911,7 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
zlab=P$main,ylim=P$ylim,xlim=P$xlim,theta=theta,phi=phi,...)
} else if (scheme==2) {
image(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
main=P$main,xlim=P$xlim,ylim=P$ylim,col=heat.colors(50),...)
main=P$main,xlim=P$xlim,ylim=P$ylim,col=colors,...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=3,...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
......@@ -953,8 +956,8 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
zlab=P$main,theta=theta,phi=phi,xlim=P$xlim,ylim=P$ylim,...)
} else if (scheme==2) {
image(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
main=P$main,xlim=P$xlim,ylim=P$ylim,col=heat.colors(50),...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=3,...)
main=P$main,xlim=P$xlim,ylim=P$ylim,col=colors,...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=contour.col,...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
......
This diff is collapsed.
......@@ -747,7 +747,8 @@ plot.soap.film <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,
pers=FALSE,theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,...) {
shift=0,trans=I,by.resids=FALSE,scheme=0,colors=heat.colors(100),
contour.col=1,...) {
## plot method function for soap.smooth terms
if (scheme==3) {
if (is.null(P)) outline <- FALSE else outline <- TRUE
......@@ -758,7 +759,8 @@ plot.soap.film <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,
pers=pers,theta=theta,phi=phi,jit=jit,xlab=xlab,ylab=ylab,main=main,
ylim=ylim,xlim=xlim,too.far=too.far,shade=shade,shade.col=shade.col,
shift=shift,trans=trans,by.resids=by.resids,scheme=scheme,...)
shift=shift,trans=trans,by.resids=by.resids,scheme=scheme,colors=colors,
contour.col=contour.col,...)
if (outline) { if (is.null(names(P$bnd))) {