Commit c7dd87ef authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.7-1

parent 5447230a
Package: mgcv
Version: 1.7-0
Version: 1.7-1
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: GAMs with GCV/AIC/REML smoothness estimation and GAMMs by PQL
......@@ -12,6 +12,6 @@ Imports: graphics, stats, nlme, Matrix
Suggests: nlme (>= 3.1-64), splines, Matrix
LazyLoad: yes
License: GPL (>= 2)
Packaged: 2010-11-02 10:24:50 UTC; simon
Packaged: 2010-11-05 16:43:24 UTC; simon
Repository: CRAN
Date/Publication: 2010-11-03 07:41:48
Date/Publication: 2010-11-08 07:31:31
......@@ -3,6 +3,20 @@
## With contributions from Henric Nilsson
Rrank <- function(R,tol=.Machine$double.eps^.9) {
## Finds rank of upper triangular matrix R, by estimating condition
## number of upper rank by rank block, and reducing rank until this is
## acceptably low... assumes R pivoted
rank <- m <- ncol(R)
ok <- TRUE
while (ok) {
Rcond <- .C(C_R_cond,R=as.double(R),r=as.integer(m),c=as.integer(rank),
work=as.double(rep(0,4*m)),Rcond=as.double(1))$Rcond
if (Rcond*tol<1) ok <- FALSE else rank <- rank - 1
}
rank
}
rig <- function(n,mean,scale) {
## inverse guassian deviates generated by algorithm 5.7 of
## Gentle, 2003. scale = 1/lambda.
......@@ -362,9 +376,11 @@ gam.side <- function(sm,Xp,tol=.Machine$double.eps^.5)
## Need to test for intercept or equivalent in Xp
intercept <- FALSE
if (ncol(Xp)) {
f <- rep(1,nrow(Xp))
ff <- qr.fitted(qr(Xp),f)
if (max(abs(ff-f))<.Machine$double.eps*100) intercept <- TRUE
if (sum(apply(Xp,2,sd)<.Machine$double.eps^.75)>0) intercept <- TRUE else {
f <- rep(1,nrow(Xp))
ff <- qr.fitted(qr(Xp),f)
if (max(abs(ff-f))<.Machine$double.eps^.75) intercept <- TRUE
}
}
sm.id <- as.list(v.names)
......@@ -739,7 +755,26 @@ gam.setup <- function(formula,pterms,data=stop("No data supplied to gam.setup"),
G$cmX <- colMeans(X) ## useful for componentwise CI construction
} else {
G$cmX <- colMeans(Xp)
G$P <- qr.coef(qr(Xp),X) ## transform from fit params to prediction params
## transform from fit params to prediction params...
G$P <- qr.coef(qr(Xp),X) ## old code assumes always full rank!!
qrx <- qr(Xp,LAPACK=TRUE)
R <- qr.R(qrx)
p <- ncol(R)
rank <- Rrank(R) ## rank of Xp/R
QtX <- qr.qty(qrx,X)[1:rank,]
if (rank<p) { ## rank deficient
R <- R[1:rank,]
qrr <- qr(t(R),tol=0)
R <- qr.R(qrr)
G$P <- forwardsolve(t(R),QtX)
} else {
G$P <- backsolve(R,QtX)
}
if (rank<p) {
G$P <- qr.qy(qrr,rbind(G$P,matrix(0,p-rank,p)))
}
G$P[qrx$pivot,] <- G$P
}
G$cmX[-(1:G$nsdf)] <- 0 ## zero the smooth parts here
G$X <- X;rm(X)
......@@ -1378,16 +1413,19 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,...) {
k<-k+1
}
}
names(object$coefficients) <- term.names # note - won't work on matrices!!
##names(object$coefficients) <- term.names # note - won't work on matrices!!
names(object$edf) <- term.names
names(object$edf1) <- term.names
if (!is.null(G$P)) {
object$coefficients <- G$P %*% object$coefficients
object$coefficients <- as.numeric(G$P %*% object$coefficients)
object$Vp <- G$P %*% object$Vp %*% t(G$P)
object$Ve <- G$P %*% object$Ve %*% t(G$P)
rownames(object$Vp) <- colnames(object$Vp) <- term.names
rownames(object$Ve) <- colnames(object$Ve) <- term.names
}
names(object$coefficients) <- term.names
object
}
......@@ -2282,6 +2320,7 @@ concurvity <- function(b,full=TRUE) {
m <- length(b$smooth)
if (m<1) stop("nothing to do for this model")
X <- model.matrix(b)
X <- X[rowSums(is.na(X))==0,]
## this step speeds up remaining computation...
X <- qr.R(qr(X,tol=0,LAPACK=FALSE))
stop <- start <- rep(1,m)
......
......@@ -334,7 +334,7 @@ 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=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
......
......@@ -2,6 +2,20 @@
*** denotes really big changes
ISSUES:
1.7-1
* post fitting constraint modification would fail if model matrix was
rank deficient until penalized. This was an issue when mixing new t2
terms with "re" type random effects. Fixed.
* plot.mrf.smooth bug fix. There was an implicit assumption that the
`polys' list was ordered in the same way as the levels of the covariate
of the smooth. fixed.
* gam.side intercept detection could occasionally fail. Improved.
* concurvity would fail if model matrix contained NA's. Fixed.
1.7-0
** `t2' alternative tensor product smooths added. These can be used with
......
......@@ -162,7 +162,7 @@ generalized additive mixed models. Biometrics 62(4):1025-1036
\examples{
# following shows how tensor pruduct deals nicely with
# following shows how tensor product deals nicely with
# badly scaled covariates (range of x 5\% of range of z )
test1<-function(x,z,sx=0.3,sz=0.4)
{ x<-x*20
......
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