Commit 91641428 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.6-2

parent fe0211b0
Package: mgcv
Version: 1.6-1
Version: 1.6-2
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: 2009-12-01 16:17:22 UTC; simon
Packaged: 2010-04-24 15:31:36 UTC; simon
Repository: CRAN
Date/Publication: 2009-12-01 17:19:02
Date/Publication: 2010-04-27 19:09:58
useDynLib(mgcv, .registration = TRUE, .fixes = "C_")
export(anova.gam, bam, cSplineDes,
export(anova.gam, bam, bam.update, cSplineDes,
exclude.too.far,extract.lme.cov, extract.lme.cov2,
formXtViX, full.score, formula.gam,fixDependence,fix.family.link,
fix.family.var, fix.family.ls, gam, gam2derivative,
gam2objective,
gamm, gam.check, gam.control,gam.fit3,
gam.fit, gam.outer, gamSim , influence.gam,
gam.fit, gam.outer,gam.vcomp, gamSim , influence.gam,
interpret.gam,
gam.side,
get.var,ldTweedie,
......@@ -25,6 +25,7 @@ export(anova.gam, bam, cSplineDes,
Predict.matrix.tprs.smooth,
Predict.matrix.ts.smooth,
Predict.matrix.pspline.smooth,
Predict.matrix.random.effect,
residuals.gam,rig, s,
smoothCon,smooth.construct,smooth.construct2,
smooth.construct.cc.smooth.spec,
......@@ -35,6 +36,7 @@ export(anova.gam, bam, cSplineDes,
smooth.construct.tp.smooth.spec,
smooth.construct.ts.smooth.spec,
smooth.construct.ps.smooth.spec,
smooth.construct.re.smooth.spec,
smooth.construct.ad.smooth.spec,
summary.gam,sp.vcov,
te,tensor.prod.model.matrix,tensor.prod.penalties,
......
......@@ -205,6 +205,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, etastart = NUL
object$Vp <- post$Vb
object$sig2 <- object$scale <- fit$scale
object$sp <- fit$sp
names(object$sp) <- names(G$sp)
class(object)<-c("gam")
}
......@@ -277,7 +278,7 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method)
G$n <- n
G$y <- mf[[gp$response]]
} else {
} else { ## n <= chunk.size
qrx <- qr.update(sqrt(G$w)*G$X,sqrt(G$w)*G$y)
}
......@@ -321,14 +322,16 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method)
} else {
}
G$smooth <- G$X <- NULL
object$gamma <- gamma;object$G <- G;object$qrx <- qrx ## to allow updating of the model
object
}
bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,na.action,
bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,na.action=na.omit,
offset=NULL,method="REML",control=gam.control(),scale=0,gamma=1,knots=NULL,
sp=NULL,min.sp=NULL,paraPen=NULL,chunk.size=10000,...)
......@@ -473,6 +476,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
object$weights <- object$prior.weights
object$xlevels <- G$xlevels
object$y <- object$model[[gp$response]]
object$NA.action <- na.action ## version to use in bam.update
gc()
## note that predict.gam assumes that it must be ok not to split the
......@@ -490,6 +494,126 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
}
bam.update <- function(b,data,chunk.size=10000) {
## update the strictly additive model `b' in the light of new data in `data'
## Need to update modelframe (b$model)
if (is.null(b$qrx)) {
stop("Model can not be updated")
}
gp<-interpret.gam(b$formula) # interpret the formula
X <- predict(b,newdata=data,type="lpmatrix",na.action=b$NA.action) ## extra part of model matrix
cnames <- names(b$coefficients)
## now get the new data in model frame form...
if ("(weights)"%in%names(b$model)) {
mf <- model.frame(gp$fake.formula,data,weights=weights,xlev=b$xlev,na.action=b$NA.action)
w <- mf[["(weights)"]]
} else {
mf <- model.frame(gp$fake.formula,data,xlev=b$xlev,na.action=b$NA.action)
w <- rep(1,nrow(mf))
}
b$model <- rbind(b$model,mf)
## get response and offset...
off.col <- attr(attr(b$model,"terms"),"offset")
if (is.null(off.col)) offset <- rep(0,nrow(mf)) else offset <- mf[,off.col]
y <- mf[,attr(attr(b$model,"terms"),"response")] - offset
## update G
b$G$y <- c(b$G$y,y)
b$G$offset <- c(b$G$offset,offset)
b$G$w <- c(b$G$w,w)
b$G$n <- nrow(b$model)
n <- b$G$n;
## update the qr decomposition...
w <- sqrt(w)
b$qrx <- mgcv:::qr.update(w*X,w*y,b$qrx$R,b$qrx$f,b$qrx$y.norm2)
## now do the refit...
rss.extra <- b$qrx$y.norm2 - sum(b$qrx$f^2)
if (b$method=="GCV"||b$method=="UBRE") method <- "GCV.Cp" else method <- b$method
if (method=="GCV.Cp") {
if (b$method=="GCV") scale <- -1 else scale = b$sig2
fit <- magic(b$qrx$f,b$qrx$R,b$sp,b$G$S,b$G$off,L=b$G$L,lsp0=b$G$lsp0,rank=b$G$rank,
H=b$G$H,C=b$G$C,gamma=b$gamma,scale=scale,gcv=(scale<=0),
extra.rss=rss.extra,n.score=n)
post <- magic.post.proc(b$qrx$R,fit,b$qrx$f*0+1)
b$y <- b$G$y;b$offset <- b$G$offset; b$G$w -> b$weights -> b$prior.weights;
} else { ## method is "REML" or "ML"
y <- b$G$y; w <- b$G$w;offset <- b$G$offset
b$G$y <- b$qrx$f
b$G$w <- b$G$y*0+1
b$G$X <- b$qrx$R
b$G$n <- length(b$G$y)
b$G$offset <- b$G$y*0
b$G$dev.extra <- rss.extra
b$G$pearson.extra <- rss.extra
b$G$n.true <- n
if (b$scale.estimated) scale <- -1 else scale = b$sig2
in.out <- list(sp=b$sp,scale=b$gcv.ubre)
object <- gam(G=b$G,method=method,gamma=b$gamma,scale=scale,in.out=in.out)
offset -> b$G$offset -> b$offset
w -> b$G$w -> b$weights -> b$prior.weights; n -> b$G$n
y -> b$G$y -> b$y;
}
if (method=="GCV.Cp") {
b$coefficients <- fit$b
b$edf <- post$edf
b$full.sp <- fit$sp.full
b$gcv.ubre <- fit$score
b$hat <- post$hat
b$mgcv.conv <- fit$gcv.info
b$optimizer="magic"
b$rank <- fit$gcv.info$rank
b$Ve <- post$Ve
b$Vp <- post$Vb
b$sig2 <- b$scale <- fit$scale
b$sp <- fit$sp
} else { ## REML or ML
b$coefficients <- object$coefficients
b$edf <- object$edf
b$full.sp <- object$sp.full
b$gcv.ubre <- object$gcv.ubre
b$hat <- object$hat
b$outer.info <- object$outer.info
b$rank <- object$rank
b$Ve <- object$Ve
b$Vp <- object$Vp
b$sig2 <- b$scale <- object$sig2
b$sp <- object$sp
}
b$G$X <- NULL
b$linear.predictors <- as.numeric(predict.gam(b,newdata=b$model,block.size=chunk.size))
b$fitted.values <- b$linear.predictor ## strictly additive only!
b$residuals <- sqrt(b$family$dev.resids(b$y,b$fitted.values,b$weights)) *
sign(b$y-b$fitted.values)
b$deviance <- sum(b$residuals^2)
b$aic <- b$family$aic(b$y,1,b$fitted.values,b$weights,b$deviance)
b$null.deviance <- sum(b$family$dev.resids(b$y,mean(b$y),b$weights))
names(b$coefficients) <- names(b$edf) <- cnames
b
} ## end of bam.update
#### ISSUES:
## ? negative binomial support --- docs say it's there...
......
......@@ -360,8 +360,10 @@ gamm.setup<-function(formula,pterms,data=stop("No data supplied to gamm.setup"),
G$Xf <- G$X # full GAM model matrix, treating smooths as fixed effects
random<-list()
random.i<-0
if (G$nsdf>0) ind <- 1:G$nsdf else ind <- rep(0,0)
X <- G$X[,1:G$nsdf,drop=FALSE] # accumulate fixed effects into here
X <- G$X[,ind,drop=FALSE] # accumulate fixed effects into here
xlab <- rep("",0)
if (G$m)
......@@ -1095,48 +1097,15 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
first <- last + 1
}
S<-S/ret$lme$sigma^2 # X'V^{-1}X divided by \sigma^2, so should S be
Vb <- chol2inv(chol(XVX+S)) # covariance matrix - in constraint space
# need to project out of constraint space
if (FALSE) { ## old code before constraint absorption
Vp <- matrix(Vb[1:G$nsdf,],G$nsdf,ncol(Vb))
X <- matrix(XVX[1:G$nsdf,],G$nsdf,ncol(XVX))
first <- G$nsdf+1
if (G$m >0) for (i in 1:G$m)
{ if (is.null(G$smooth[[i]]$C)) nc <- 0 else nc <- nrow(G$smooth[[i]]$C)
last <- first+object$smooth[[i]]$df-1 # old code: nrow(object$smooth[[i]]$ZSZ)-1
{ V <- rbind(matrix(0,nc,ncol(Vp)),Vb[first:last,])
if (nc) V <- qr.qy(G$smooth[[i]]$qrc,V)
Vp <- rbind(Vp,V) # cov matrix
V <- rbind(matrix(0,nc,ncol(Vp)),XVX[first:last,])
if (nc) V <- qr.qy(G$smooth[[i]]$qrc,V)
X <- rbind(X,V) # X'V^{-1}X
}
first <- last+1
}
Vb<-matrix(Vp[,1:G$nsdf],nrow(Vp),G$nsdf)
Z <- matrix(X[,1:G$nsdf],nrow(X),G$nsdf)
first<-G$nsdf+1
if (G$m>0) for (i in 1:G$m)
{ last <- first + object$smooth[[i]]$df-1
if (is.null(G$smooth[[i]]$C)) nc <- 0 else nc <- nrow(G$smooth[[i]]$C)
{ V <- cbind(matrix(0,nrow(Vb),nc),Vp[,first:last])
if (nc) V <- t(qr.qy(G$smooth[[i]]$qrc,t(V)))
Vb<-cbind(Vb,V)
V <- cbind(matrix(0,nrow(Vb),nc),X[,first:last])
if (nc) V <- t(qr.qy(G$smooth[[i]]$qrc,t(V)))
Z<-cbind(Z,V)
}
first <- last+1
}
# then get edf diag(Vb%*%Z)
object$edf<-rowSums(Vb*t(Z)) } else {
object$edf<-rowSums(Vb*t(XVX))
}
## Vb <- chol2inv(chol(XVX+S)) # covariance matrix - in constraint space
## replacement, more stable version...
ev <- eigen(XVX+S,symmetric=TRUE)
ind <- ev$values != 0
iv <- ev$values;iv[ind] <- 1/ev$values[ind]
Vb <- ev$vectors%*%(iv*t(ev$vectors))
object$edf<-rowSums(Vb*t(XVX))
object$sig2 <- ret$lme$sigma^2
if (lme.used) { object$method <- paste("lme.",method,sep="")}
......
......@@ -567,7 +567,8 @@ gam.setup <- function(formula,pterms,data=stop("No data supplied to gam.setup"),
temp.term <- split$smooth.spec[[i]]$term
for (j in 1:length(temp.term)) id.list[[id]]$data[[j]] <- cbind(id.list[[id]]$data[[j]],
get.var(temp.term[j],data,vecMat=FALSE))
} else { ## new id
} else { ## new id
id.list[[id]] <- list(sm.i=i) ## start the array of indices of smooths with this id
id.list[[id]]$data <- list()
## need to collect together all data for which this basis will be used,
......@@ -669,7 +670,7 @@ gam.setup <- function(formula,pterms,data=stop("No data supplied to gam.setup"),
## min.sp must be length nrow(L) to make sense
## sp must be length ncol(L) --- need to partition
## L into columns relating to free log smoothing paramters,
## L into columns relating to free log smoothing parameters,
## and columns, L0, corresponding to values supplied in sp.
## lsp0 = L0%*%log(sp[sp>=0]) [need to fudge sp==0 case by
## setting log(0) to, e.g. 10*log(.Machine$double.xmin)]
......@@ -985,7 +986,6 @@ gam.negbin <- function(lsp,fscale,family,control,method,optimizer,gamma,G,scale,
gam.outer <- function(lsp,fscale,family,control,method,optimizer,criterion,scale,gamma,G,...)
# function for smoothing parameter estimation by outer optimization. i.e.
# P-IRLS scheme iterated to convergence for each trial set of smoothing
......@@ -1157,6 +1157,13 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,...) {
else G$sig2 <- -1 #gcv
} else {G$sig2 <- scale}
if (fam.name=="quasi"||fam.name=="quasipoisson"||fam.name=="quasibinomial") {
## REML/ML invalid with quasi families
if (method=="REML") method <- "P-REML"
if (method=="ML") method <- "P-ML"
if (method=="P-ML"||method=="P-REML") warning("RE/ML is not recommended for quasi families")
}
if (reml) { ## then RE(ML) selection, but which variant?
criterion <- method
if (fam.name == "binomial"||fam.name == "poisson") scale <- 1
......@@ -1171,11 +1178,6 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,...) {
}
}
if (fam.name=="quasi"||fam.name=="quasipoisson"||fam.name=="quasibinomial") {
## REML/ML invalid with quasi families
if (method=="REML") method <- "P-REML"
if (method=="ML") method <- "P-ML"
}
if (substr(fam.name,1,17)=="Negative Binomial") {
scale <- 1; ## no choise
......@@ -1459,16 +1461,16 @@ gam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
}
gam.check <- function(b)
gam.check <- function(b,...)
# takes a fitted gam object and produces some standard diagnostic plots
{ if (b$method%in%c("GCV","GACV","UBRE","REML","ML","P-ML","P-REML"))
{ old.par<-par(mfrow=c(2,2))
sc.name<-b$method
qqnorm(residuals(b))
qqnorm(residuals(b),...)
plot(b$linear.predictors,residuals(b),main="Resids vs. linear pred.",
xlab="linear predictor",ylab="residuals");
hist(residuals(b),xlab="Residuals",main="Histogram of residuals");
plot(fitted(b),b$y,xlab="Fitted Values",ylab="Response",main="Response vs. Fitted Values")
xlab="linear predictor",ylab="residuals",...);
hist(residuals(b),xlab="Residuals",main="Histogram of residuals",...);
plot(fitted(b),b$y,xlab="Fitted Values",ylab="Response",main="Response vs. Fitted Values",...)
## now summarize convergence information
cat("\nMethod:",b$method," Optimizer:",b$optimizer)
......@@ -1504,7 +1506,7 @@ gam.check <- function(b)
cat("\n")
par(old.par)
} else ## probably a `gamm' `gam' object
plot(b$linear.predictor,residuals(b),xlab="linear predictor",ylab="residuals")
plot(b$linear.predictor,residuals(b),xlab="linear predictor",ylab="residuals",...)
}
print.gam<-function (x,...)
......@@ -2325,7 +2327,9 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
}
} ## end of sp.contour
# start of main function
#########################
## start of main function
#########################
w.resid<-NULL
if (length(residuals)>1) # residuals supplied
{ if (length(residuals)==length(x$residuals))
......@@ -2351,7 +2355,9 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
# plot should ignore all "by" variables
# sort out number of pages and plots per page
n.plots <- m + n.para
n.plots <- n.para
if (m>0) for (i in 1:m) n.plots <- n.plots + as.numeric(x$smooth[[i]]$plot.me)
if (pages>n.plots) pages<-n.plots
if (pages<0) pages<-0
if (pages!=0) # figure out how to display things
......@@ -2395,7 +2401,10 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
}
pd<-list();
i<-1 # needs a value if no smooths, but parametric terms ...
## First the loop to get the data for the plots...
if (m>0) for (i in 1:m) # work through smooth terms
if (x$smooth[[i]]$plot.me)
{ if (x$smooth[[i]]$dim==1)
{ raw<-x$model[x$smooth[[i]]$term]
xx<-seq(min(raw),max(raw),length=n) # generate x sequence for prediction
......@@ -2482,14 +2491,15 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
pd[[i]]<-pd.item;rm(pd.item)
} else
{ pd[[i]]<-list(dim=x$smooth[[i]]$dim)}
}
} ## end of loop creating plot data
# now plot .....
## now plot .....
if (se) # pd$fit and pd$se
{ k<-0
if (scale==-1&&is.null(ylim)) # getting common scale for 1-d terms
if (m>0) for (i in 1:m)
if (x$smooth[[i]]$plot.me)
{ if (pd[[i]]$dim==1)
{ ul<-pd[[i]]$fit+pd[[i]]$se
ll<-pd[[i]]$fit-pd[[i]]$se
......@@ -2509,6 +2519,7 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
}
j<-1
if (m>0) for (i in 1:m)
if (x$smooth[[i]]$plot.me)
{ if (is.null(select)||i==select)
{ ##if (interactive()&& is.null(select) && pd[[i]]$dim<3 && i>1&&(i-1)%%ppp==0)
##readline("Press return for next page....")
......@@ -2636,6 +2647,8 @@ plot.gam <- function(x,residuals=FALSE,rug=TRUE,se=TRUE,pages=0,select=NULL,scal
j<-j+pd[[i]]$dim
}
}
if (n.para>0) # plot parameteric terms
{ class(x) <- c("gam","glm","lm") # needed to get termplot to call model.frame.glm
if (is.null(select)) {
......@@ -3000,6 +3013,75 @@ sp.vcov <- function(x) {
} else return(NULL)
}
gam.vcomp <- function(x,rescale=TRUE,conf.lev=.95) {
## Routine to convert smoothing parameters to variance components
## 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 (rescale) { ## undo any rescaling of S[[i]] that may have been done
k <- 1;m <- length(x$smooth)
idx <- rep("",0)
if (m>0) for (i in 1:m) { ## loop through all smooths
if (!is.null(x$smooth[[i]]$id)) { ## smooth has an id
if (x$smooth[[i]]$id%in%idx) {
ok <- FALSE ## id already dealt with --- ignore smooth
} else {
idx <- c(idx,x$smooth[[i]]$id) ## add id to id list
ok <- TRUE
}
} else { ok <- TRUE} ## no id so proceed
if (ok) for (j in 1:length(x$smooth[[i]]$S.scale)) {
x$sp[k] <- x$sp[k] / x$smooth[[i]]$S.scale[j]
k <- k + 1
}
} ## finished rescaling
}
## variance components (original scale)
vc <- c(scale/x$sp)
names(vc) <- names(x$sp)
## If a Hessian exists, get CI's for variance components...
if (x$method%in%c("ML","P-ML","REML","P-REML")&&!is.null(x$outer.info$hess)) {
H <- x$outer.info$hess ## the hessian w.r.t. log sps and log scale
if (ncol(H)>length(x$sp)) scale.est <- TRUE else scale.est <- FALSE
## get derivs of log sqrt var comps wrt log sp and log scale....
J <- matrix(0,nrow(H),ncol(H))
if (scale.est) {
diag(J) <- -2
J[,ncol(J)] <- 2
vc <- c(vc,scale);names(vc) <- c(names(x$sp),"scale")
} else {
diag(J) <- -0.5
}
H <- t(J)%*%H%*%J ## hessian of log sqrt variances
eh <- eigen(H,symmetric=TRUE)
ind <- eh$values>max(eh$values)*.Machine$double.eps^75 ## index of non zero eigenvalues
rank <- sum(ind) ## rank of hessian
iv <- eh$values*0;iv[ind] <- 1/eh$values[ind]
V <- eh$vectors%*%(iv*t(eh$vectors)) ## cov matrix for log sqrt variances
lsd <- log(sqrt(vc)) ## log sqrt variances
sd.lsd <- sqrt(diag(V))
if (conf.lev<=0||conf.lev>=1) conf.lev <- 0.95
crit <- qnorm(1-(1-conf.lev)/2)
ll <- lsd - crit * sd.lsd
ul <- lsd + crit * sd.lsd
res <- cbind(exp(lsd),exp(ll),exp(ul))
rownames(res) <- names(vc)
colnames(res) <- c("std.dev","lower","upper")
cat("\n")
cat(paste("Standard deviations and",conf.lev,"confidence intervals:\n\n"))
print(res)
cat("\nRank: ");cat(rank);cat("/");cat(ncol(H));cat("\n")
invisible(res)
} else {
return(vc)
}
} ## end of gam.vcomp
vcov.gam <- function(object, freq = FALSE, dispersion = NULL, ...)
## supplied by Henric Nilsson <henric.nilsson@statisticon.se>
{ if (freq)
......@@ -3413,7 +3495,7 @@ initial.spg <- function(X,y,weights,family,S,off,L=NULL,lsp0=NULL,type=1) {
start <- etastart <- mustart <- NULL
nobs <- nrow(X)
eval(family$initialize)
w <- weights*family$mu.eta(family$linkfun(mustart))^2/family$variance(mustart)
w <- as.numeric(weights*family$mu.eta(family$linkfun(mustart))^2/family$variance(mustart))
w <- sqrt(w)
if (type==1) { ## what PI would have used
lambda <- initial.sp(w*X,S,off)
......@@ -3617,7 +3699,8 @@ print.mgcv.version <- function()
version <- version[pmatch("Version",version)]
um <- strsplit(version," ")[[1]]
version <- um[nchar(um)>0][2]
cat(paste("This is mgcv ",version,". For overview type `help(\"mgcv-package\")'.\n",sep=""))
hello <- paste("This is mgcv ",version,". For overview type 'help(\"mgcv-package\")'.",sep="")
packageStartupMessage(hello)
}
set.mgcv.options <- function()
......
......@@ -316,6 +316,7 @@ smooth.construct.tensor.smooth.spec<-function(object,data,knots)
}
object$margin[[i]]<-smooth.construct(object$margin[[i]],dat,knt)
Xm[[i]]<-object$margin[[i]]$X
if (!is.null(object$margin[[i]]$te.ok) && !object$margin[[i]]$te.ok) stop("attempt to use unsuitable marginal smooth class")
if (length(object$margin[[i]]$S)>1)
stop("Sorry, tensor products of smooths with multiple penalties are not supported.")
Sm[[i]]<-object$margin[[i]]$S[[1]]
......@@ -1122,11 +1123,59 @@ smooth.construct.ad.smooth.spec<-function(object,data,knots)
} ## adaptive penalty finished
} ## penalized case finished
}
pspl$te.ok <- FALSE ## not suitable as a tensor product marginal
pspl
}
#################################
# Random effects terms start here
#################################
smooth.construct.re.smooth.spec<-function(object,data,knots)
## a simple random effects constructor method function
## basic idea is that s(x,f,z,...,bs="re") generates model matrix
## corresponding to ~ x:f:z: ... - 1. Corresponding coefficients
## have an identity penalty.
{
## id's with factor variables are problematic - should terms have
## same levels, or just same number of levels, for example?
## => ruled out
if (!is.null(object$id)) stop("random effects don't work with ids.")
form <- as.formula(paste("~",paste(object$term,collapse=":"),"-1"))
object$X <- model.matrix(form,data)
object$bs.dim <- ncol(object$X)
## now construct penalty
object$S <- list(diag(object$bs.dim)) # get penalty
object$rank <- object$bs.dim # penalty rank
object$null.space.dim <- 0 # dimension of unpenalized space
object$C <- matrix(0,0,ncol(object$X)) # null constraint matrix
## need to store formula (levels taken care of by calling function)
object$form <- form
object$plot.me <- FALSE ## "re" terms should not be plotted by plot.gam
object$te.ok <- FALSE ## these terms are not suitable as te marginals
class(object)<-"random.effect" # Give object a class
object
}
Predict.matrix.random.effect<-function(object,data)
# prediction method function for the p.spline smooth class
{ require(splines)
X <- model.matrix(object$form,data)
X
}
......@@ -1251,6 +1300,10 @@ smoothCon <- function(object,data,knots,absorb.cons=FALSE,scale.penalty=TRUE,n=n
{ sm <- smooth.construct3(object,data,knots)
if (!is.null(attr(sm,"qrc"))) warning("smooth objects should not have a qrc attribute.")
## add plotting indicator if not present.
## plot.me tells `plot.gam' whether or not to plot the term
if (is.null(sm$plot.me)) sm$plot.me <- TRUE
## automatically produce centering constraint...
## must be done here on original model matrix to ensure same
## basis for all `id' linked terms
......@@ -1293,12 +1346,16 @@ smoothCon <- function(object,data,knots,absorb.cons=FALSE,scale.penalty=TRUE,n=n
## but to do otherwise would mess up the meaning of smoothing parameters
## sufficiently that linking terms via `id's would not work properly (they
## would have the same basis, but different penalties)
sm$S.scale <- rep(1,length(sm$S))
if (scale.penalty && length(sm$S)>0 && is.null(sm$no.rescale)) # then the penalty coefficient matrix is rescaled
{ maXX <- mean(abs(t(sm$X)%*%sm$X)) # `size' of X'X
for (i in 1:length(sm$S)) {
maS <- mean(abs(sm$S[[i]]))
sm$S[[i]] <- sm$S[[i]] * maXX / maS
}
maS <- mean(abs(sm$S[[i]])) / maXX
sm$S[[i]] <- sm$S[[i]] / maS
sm$S.scale[i] <- maS ## multiply S[[i]] by this to get original S[[i]]
}
}
## check whether different data to be used for basis setup
......
** denotes quite substantial/important changes
*** denotes really big changes
1.6-2
** Random effect support for `gam' improved with the addition of
a random effect "smooth" class, and a function for extracting
variance components. See ?random.effects and links for details.
* smooths now contain extra elements: S.scale records the scale factor
used in any linear rescaling of a penalty matrix; plot.me indicates
whether `plot.gam' should attempt to plot the term; te.ok indicates
whether the smooth is a suitable marginal for a tensor product.
* Fix in `gamm.setup' -- models with no fixed effects (except smooths)
could fail to fit properly, because of an indexing error (caused odd
results with indicator by variables)
* help files have had various misuses on `itemize' fixed.
* initialization could fail if response was actually a 1D array. fixed.
* New function `bam.update' allows efficient updating of very large
strictly additive models fitted by `bam' when additional data become
available.
* gam now warns if RE/ML used with quasi families.
* gam.check now accepts graphics parameters.
* fixed problem in welcome message that messed up ESS.
1.6-1
* Bug in cSplineDes caused bases to be non-cyclic unless first
......
......@@ -2,7 +2,7 @@
\alias{anova.gam}
\alias{print.anova.gam}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Hypothesis tests related to GAM fits}
\title{Approximate hypothesis tests related to GAM fits}
\description{ Performs hypothesis tests relating to one or more fitted
\code{gam} objects. For a single fitted \code{gam} object, Wald tests of
the significance of each parametric and smooth term are performed. Otherwise
......
......@@ -15,7 +15,7 @@ for large datasets.
}
\usage{
bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
na.action, offset=NULL,method="REML",control=gam.control(),
na.action=na.omit, offset=NULL,method="REML",control=gam.control(),
scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL,paraPen=NULL,
chunk.size=10000,...)
}
......
\name{bam.update}
\alias{bam.update}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Update a strictly additive bam model for new data.}
\description{ Strictly additive models fitted by \code{\link{bam}} can be efficiently updated as new data becomes available,
by simply updateing the QR decomposition on which estimation is based, and re-optimizing the smoothing parameters, starting
from the previous estimates. This routine implements this.
}
\usage{
bam.update(b,data,chunk.size=10000)
}
%- maybe also `usage' for other objects documented here.