Commit f31b3bf5 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.7-20

parent b568af2d
Package: mgcv
Version: 1.7-19
Version: 1.7-20
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
......@@ -14,6 +14,6 @@ Suggests: nlme (>= 3.1-64), splines, Matrix, parallel
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2012-07-18 06:24:25 UTC; sw283
Packaged: 2012-08-24 09:38:21 UTC; sw283
Repository: CRAN
Date/Publication: 2012-07-20 05:00:31
Date/Publication: 2012-08-24 19:00:39
882792073c8f32da732e6450ae7d1c12 *DESCRIPTION
7b0afac092bc0b2f4a344d67a1c7eaf0 *DESCRIPTION
50152051f123a389d421aa3130dce252 *NAMESPACE
7b5db03a3a80878eb8007f8c0848583a *R/bam.r
e59ee1842b4546cb64ebe39e6a7f00ee *R/fast-REML.r
2073e50d30aacc57abac1db8b3133467 *R/gam.fit3.r
902657a0ee2dedc3fdfa501bf3b37c5b *R/gam.sim.r
3527e64fa9d32baba5644c69c7d1eaf6 *R/gamm.r
41fc0ffd73b6f2ee4e18071134216299 *R/mgcv.r
958a9a13b359c6e6f23f340a05d98899 *R/plots.r
77f6910e5b8b1aeb997006ef7ad20b5c *R/smooth.r
7abad0459a93e7e451e018d0109612a8 *R/mgcv.r
1a5634ab06b9f0bd68fd86e489e9c864 *R/plots.r
d16c7de3c76cfbc1dfc418f49bc94eb5 *R/smooth.r
fb66d6c18398411a99ffcb788b854f13 *R/sparse.r
f97ab579ec183613a62c338e34e31e61 *changeLog
7a549e348bd9f2eeb51dccb1dc5de6d4 *changeLog
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
88d77139cc983317b6acd8c5f1252ab9 *gnugpl2.txt
......@@ -17,7 +17,7 @@ f693920e12f8a6f2b6cab93648628150 *index
9388a0979dab4fe1988b777ae4dcfb0a *inst/CITATION
5bf12ddc0dab9daae72271b96a15c539 *man/Predict.matrix.Rd
c51c9b8c9c73f81895176ded39b91394 *man/Predict.matrix.cr.smooth.Rd
612ab6354541ebe38a242634d73b66ba *man/Tweedie.Rd
3e5e30b44947d3ddf00fcd55be219038 *man/Tweedie.Rd
0d24940b2e0a3e9fa7b3dc0224d049ab *man/anova.gam.Rd
585d6b6a2e7816df83d5f9193195a233 *man/bam.Rd
b385d6d5419d0d6aefe03af1a79d5c4e *man/bam.update.Rd
......@@ -86,12 +86,12 @@ f77ca1471881d2f93c74864d076c0a0e *man/rTweedie.Rd
7258dfc0149fff020054117fd2ee6bd8 *man/s.Rd
690dddb35fd5986a6eeb39dd79fa33f9 *man/slanczos.Rd
3026df475af053acd2b742af875f626c *man/smooth.construct.Rd
10a456bc3d7151f60be24a7e3bba3d30 *man/smooth.construct.ad.smooth.spec.Rd
9f05449c0114748608f99693ec25cf3f *man/smooth.construct.cr.smooth.spec.Rd
747f1bdf9251d6cf95653aec67b02599 *man/smooth.construct.ad.smooth.spec.Rd
c07b9d5733913e8397a91bd2526f0d23 *man/smooth.construct.cr.smooth.spec.Rd
131c3f2131138ba5d6f6bcde4be5ac31 *man/smooth.construct.ds.smooth.spec.Rd
1bb6748d4d2934e48f0572bc5114ffcb *man/smooth.construct.fs.smooth.spec.Rd
772e6c18d25cfaf3d9c194031ee042fe *man/smooth.construct.mrf.smooth.spec.Rd
78688702b6396f24692b74c554483428 *man/smooth.construct.ps.smooth.spec.Rd
65f84caf46cc1abbfa8ef8d9f3b7bbe2 *man/smooth.construct.ps.smooth.spec.Rd
d202c6718fb1138fdd99e6102250aedf *man/smooth.construct.re.smooth.spec.Rd
53d7986dd7a54c1edb13ce84dbfe34a2 *man/smooth.construct.sos.smooth.spec.Rd
8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
......
......@@ -1406,7 +1406,7 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,...) {
## correct null deviance if there's an offset [Why not correct calc in gam.fit/3???]....
if (G$intercept&&any(G$offset!=0)) object$null.deviance <-
glm(G$y~offset(G$offset),family=object$family,weights=object$prior.weights)$deviance
glm(object$y~offset(G$offset),family=object$family,weights=object$prior.weights)$deviance
object$method <- criterion
......@@ -2542,10 +2542,8 @@ recov <- function(b,re=rep(0,0),m=0) {
## dropped.
if (!inherits(b,"gam")) stop("recov works with fitted gam objects only")
if (is.null(b$full.sp)) sp <- b$sp else sp <- b$full.sp
llr <- TRUE
if (length(re)<1) {
if (m>0) {
if (llr) { ## llr
## annoyingly, need total penalty
np <- length(coef(b))
k <- 1;S1 <- matrix(0,np,np)
......@@ -2564,21 +2562,8 @@ recov <- function(b,re=rep(0,0),m=0) {
LRB <- cbind(LRB[,-ii],LRB[,ii])
ii <- (ncol(LRB)-length(ii)+1):ncol(LRB)
Rm <- qr.R(qr(LRB,tol=0,LAPACK=FALSE))[ii,ii] ## unpivoted QR
} else {
er <- eigen(crossprod(b$R))
ii <- er$values>max(er$values)*.Machine$double.eps^0.8
er$values[!ii] <- 0
er$values[ii] <- 1/sqrt(er$values[ii])
ind <- b$smooth[[m]]$first.para:b$smooth[[m]]$last.para
Vu <- crossprod(er$values*t(er$vector))[ind,ind]
er <- eigen(Vu,symmetric=TRUE)
ii <- er$values>max(er$values)*.Machine$double.eps^0.8
er$values[!ii] <- 0
er$values[ii] <- 1/sqrt(er$values[ii])
Rm <- er$values*t(er$vectors)
}
} else Rm <- NULL
return(list(Ve=b$Ve,Rm=Rm))
return(list(Ve=(t(b$Ve)+b$Ve)*.5,Rm=Rm))
}
if (m%in%re) stop("m can't be in re")
......@@ -2615,25 +2600,27 @@ recov <- function(b,re=rep(0,0),m=0) {
}
## pseudoinvert S2
if (nrow(S2)==1) {
S2[1,1] <- 1/S2[1,1]
S2[1,1] <- 1/sqrt(S2[1,1])
} else if (max(abs(diag(diag(S2))-S2))==0) {
ds2 <- diag(S2)
ind <- ds2 > max(ds2)*.Machine$double.eps^.8
ds2[ind] <- 1/ds2[ind];ds2[!ind] <- 0
diag(S2) <- ds2
diag(S2) <- sqrt(ds2)
} else {
ev <- eigen(S2,symmetric=TRUE)
ev <- eigen((S2+t(S2))/2,symmetric=TRUE)
ind <- ev$values > max(ev$values)*.Machine$double.eps^.8
ev$values[ind] <- 1/ev$values[ind];ev$values[!ind] <- 0
S2 <- ev$vectors%*%(ev$values*t(ev$vectors))
## S2 <- ev$vectors%*%(ev$values*t(ev$vectors))
S2 <- sqrt(ev$values)*t(ev$vectors)
}
## choleski of cov matrix....
L <- chol(diag(p)+R2%*%S2%*%t(R2)) ## L'L = I + R2 S2^- R2'
## L <- chol(diag(p)+R2%*%S2%*%t(R2)) ## L'L = I + R2 S2^- R2'
L <- chol(diag(p) + crossprod(S2%*%t(R2)))
## now we need the square root of the unpenalized
## cov matrix for m
if (m>0) {
if (llr) { ## llr version
## llr version
LRB <- rbind(L%*%R1,t(mroot(S1)))
ii <- map[b$smooth[[m]]$first.para:b$smooth[[m]]$last.para]
## ii is cols of LRB related to smooth m, which need
......@@ -2641,19 +2628,6 @@ recov <- function(b,re=rep(0,0),m=0) {
LRB <- cbind(LRB[,-ii],LRB[,ii])
ii <- (ncol(LRB)-length(ii)+1):ncol(LRB) ## need to pick up final block
Rm <- qr.R(qr(LRB,tol=0,LAPACK=FALSE))[ii,ii] ## unpivoted QR
} else { ## original inverse unpenalized cov version
er <- eigen(crossprod(L%*%R1),symmetric=TRUE)
ii <- er$values>max(er$values)*.Machine$double.eps^0.8
er$values[!ii] <- 0
er$values[ii] <- 1/sqrt(er$values[ii])
ind <- map[b$smooth[[m]]$first.para:b$smooth[[m]]$last.para]
Vu <- crossprod(er$values*t(er$vector))[ind,ind]
er <- eigen(Vu,symmetric=TRUE)
ii <- er$values>max(er$values)*.Machine$double.eps^0.8
er$values[!ii] <- 0
er$values[ii] <- 1/sqrt(er$values[ii])
Rm <- er$values*t(er$vectors)
}
} else Rm <- NULL
list(Ve= crossprod(L%*%b$R%*%b$Vp)/b$sig2, ## Frequentist cov matrix
......@@ -2664,7 +2638,7 @@ recov <- function(b,re=rep(0,0),m=0) {
reTest <- function(b,m) {
## Test the mth smooth for equality to zero, using f'f as statistic,
## Test the mth smooth for equality to zero
## and accounting for all random effects in model
## find indices of random effects other than m
......@@ -2675,24 +2649,9 @@ reTest <- function(b,m) {
Ve <- rc$Ve
ind <- b$smooth[[m]]$first.para:b$smooth[[m]]$last.para
B <- mroot(Ve[ind,ind]) ## BB'=Ve
if (FALSE) {
## (R'R)^- is unpenalized cov matrix
ev <- eigen(crossprod(b$R),symmetric=TRUE)
ii <- ev$values>max(ev$values)*.Machine$double.eps^.8
ev$values[ii] <- 1/ev$values[ii]
ev$values[!ii] <- 0
Vu <- crossprod(sqrt(ev$values)*t(ev$vectors))
ev <- eigen(Vu[ind,ind])
ii <- ev$values>max(ev$values)*.Machine$double.eps^.8
ev$values[ii] <- 1/ev$values[ii]
ev$values[!ii] <- 0
Rm <- sqrt(ev$values)*t(ev$vectors) ## Rm'Rm is inv unp cov
} else {
#Rm <- b$R[,ind]
#Rm <- t(mroot(b$smooth[[m]]$S[[1]])) ## CHECK: te etc.!!
#Rm <- t(mroot(b$Ve[ind,ind]/b$sig2))
Rm <- rc$Rm
}
b.hat <- coef(b)[ind]
d <- Rm%*%b.hat
stat <- sum(d^2)/b$sig2
......
......@@ -574,7 +574,7 @@ polys.plot <- function(pc,z=NULL,scheme="heat",lab="",...) {
xlim[1] <- xlim[1] - .1 * (xlim[2]-xlim[1]) ## allow space for scale
n.col <- 100
if (scheme=="heat") scheme <- heat.colors(n.col) else
if (scheme=="heat") scheme <- heat.colors(n.col+1) else
scheme <- gray(0:n.col/n.col)
zlim <- range(pretty(z))
......@@ -588,7 +588,7 @@ polys.plot <- function(pc,z=NULL,scheme="heat",lab="",...) {
ylim <- zlim
plot(0,0,ylim=ylim,xlim=xlim,type="n",xaxt="n",bty="n",xlab="",ylab=lab,...)
for (i in 1:length(pc)) {
coli <- round((z[i] - zlim[1])/(zlim[2]-zlim[1])*100)
coli <- round((z[i] - zlim[1])/(zlim[2]-zlim[1])*n.col)+1
poly2(pc[[i]],col=scheme[coli])
}
......@@ -621,7 +621,8 @@ plot.mrf.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (!x$plot.me||is.null(x$xt$polys)) return(NULL) ## shouldn't or can't plot
## get basic plot data
raw <- data[x$term][[1]]
dat <- data.frame(x=factor(names(x$xt$polys),levels=levels(x$knots)));names(dat) <- x$term
dat <- data.frame(x=factor(names(x$xt$polys),levels=levels(x$knots)))
names(dat) <- x$term
X <- PredictMat(x,dat) # prediction matrix for this term
if (is.null(xlab)) xlabel<- "" else xlabel <- xlab
if (is.null(ylab)) ylabel <- "" else ylabel <- ylab
......@@ -1161,7 +1162,7 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
{ class(x) <- c("gam","glm","lm") # needed to get termplot to call model.frame.glm
if (is.null(select)) {
attr(x,"para.only") <- TRUE
termplot(x,se=se,rug=rug,col.se=1,col.term=1)
termplot(x,se=se,rug=rug,col.se=1,col.term=1,...)
} else { # figure out which plot is required
if (select > m) {
select <- select - m # i.e. which parametric term
......@@ -1169,7 +1170,7 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
term.labels <- term.labels[order==1]
if (select <= length(term.labels)) {
# if (interactive() && m &&i%%ppp==0)
termplot(x,terms=term.labels[select],se=se,rug=rug,col.se=1,col.term=1)
termplot(x,terms=term.labels[select],se=se,rug=rug,col.se=1,col.term=1,...)
}
}
}
......
......@@ -1909,14 +1909,6 @@ smooth.construct.re.smooth.spec <- function(object,data,knots)
object$te.ok <- 2 ## these terms are suitable as te marginals, but
## can not be plotted
## The Nychka inversion argument for p-values does not apply to
## random effects, because the eigen spectrum of the penalty is too
## uniform for the truncation argument to work in all cases,
## and because the interpretation of the Bayesian cov matrix
## as including accross the function bias is also not right,
## since these are proper random effects (i.e. what is there to
## be biased here)! The theoretical problems can result in very
## low power in practice.
object$random <- TRUE ## treat as a random effect for p-value comp.
......@@ -1928,7 +1920,7 @@ smooth.construct.re.smooth.spec <- function(object,data,knots)
Predict.matrix.random.effect<-function(object,data)
# prediction method function for the p.spline smooth class
# prediction method function for the random effect class
{ X <- model.matrix(object$form,data)
X
}
......@@ -2028,7 +2020,8 @@ smooth.construct.mrf.smooth.spec <- function(object, data, knots) {
if (sum(!levels(x)%in%levels(k)))
stop("data contain regions that are not contained in the knot specification")
levels(x) <- levels(k) ## to allow for regions with no data
##levels(x) <- levels(k) ## to allow for regions with no data
x <- factor(x,levels=levels(k)) ## to allow for regions with no data
object$X <- model.matrix(~x-1,) ## model matrix
......@@ -2119,8 +2112,8 @@ smooth.construct.mrf.smooth.spec <- function(object, data, knots) {
Predict.matrix.mrf.smooth<-function(object, data) {
x <- as.factor(data[[object$term]])
levels(x) <- levels(object$knots)
x <- factor(data[[object$term]],levels=levels(object$knots))
##levels(x) <- levels(object$knots)
X <- model.matrix(~x-1)
if (!is.null(object$P)) X <- X%*%object$P
X
......
......@@ -8,6 +8,26 @@ ISSUES:
* still not really happy with nesting constraints, especially under rank
deficiency.
1.7-20
* '...' now passd to termplot by plot.gam (thanks Andreas Eckner).
* fix to null deviance computation for binomial when n>1, matrix response
used and an offset is present. (Thanks to Tim Miller)
* Some pruning of unused code from recov and reTest.
* recov modified to stop it returning a numerically non-symmetric Ve, and
causing occasional failures of summary.gam with "re" terms.
* MRF smooth bug. Region ordering could become confused under some
circumstances due to incorrect setting of factor levels. Corrected
thanks to detailed bug report from Andreas Bender.
* polys.plot colour/grey scale bug. Could ask for colour 0 from colour
scheme, and therefore fail. Fixed.
1.7-19
** summary.gam and anova.gam now use an improved p-value computation for
......
......@@ -50,7 +50,7 @@ modified from Venables and Ripley's \code{negative.binomial} family.
}
\references{
Dunn, P.K. and G.K. Smith (2005) Series evaluation of Tweedie exponential dispersion model densities.
Dunn, P.K. and G.K. Smyth (2005) Series evaluation of Tweedie exponential dispersion model densities.
Statistics and Computing 15:267-280
Tweedie, M. C. K. (1984). An index which distinguishes between
......
......@@ -11,8 +11,8 @@ p-spline(s) or cubic regression spline(s). Discrete P-spline type penalties are
of the basis, but the penalties themselves have a basis representation, allowing the strength of the
penalty to vary with the covariates. The coefficients of the penalty basis are the smoothing parameters.
When invoking an adaptive smoother the \code{k} argument specifies the dimension of the smoothing basis,
while the \code{m} argument specifies the dimension of the penalty basis. For an adaptive smooth of two
When invoking an adaptive smoother the \code{k} argument specifies the dimension of the smoothing basis (default 40 in 1D, 15 in 2D),
while the \code{m} argument specifies the dimension of the penalty basis (default 5 in 1D, 3 in 2D). For an adaptive smooth of two
variables \code{k} is taken as the dimension of both marginal bases: different marginal basis dimensions can be
specified by making \code{k} a two element vector. Similarly, in the two dimensional case \code{m} is the
dimension of both marginal bases for the penalties, unless it is a two element vector, which specifies different
......
......@@ -15,6 +15,8 @@ straight line.) \code{s(x,bs="cc")} specifies a cyclic penalized cubic regressio
`Cardinal' spline bases are used: Wood (2006) sections 4.1.2 and 4.1.3 gives full details. These bases have
very low setup costs. For a given basis dimension, \code{k}, they typically perform a little less well
then thin plate regression splines, but a little better than p-splines. See \code{\link{te}} to use these bases in tensor product smooths of several variables.
Default \code{k} is 10.
}
\usage{
......
......@@ -13,6 +13,7 @@ the basis coefficients. Cyclic P-splines are specified by model terms like \code
These bases can be used in tensor product smooths (see \code{\link{te}}).
The advantage of P-splines is the flexible way that penalty and basis order can be mixed. This often provides a useful way of `taming' an otherwise poorly behave smooth. However, in regular use, splines with derivative based penalties (e.g. \code{"tp"} or \code{"cr"} bases) tend to result in slightly better MSE performance, presumably because the good approximation theoretic properties of splines are rather closely connected to the use of derivative penalties.
}
\usage{
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment