Commit 93453c8f authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Debian changes 1.8-13-1

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

  * New upstream release
parents 59ffb32e fe6ab539
...@@ -2,10 +2,98 @@ ...@@ -2,10 +2,98 @@
*** denotes really big changes *** denotes really big changes
Currently deprecated and liable to be removed: Currently deprecated and liable to be removed:
- bam(...,sparse=TRUE) [1.8-5]
- single penalty tensor product smooths. - single penalty tensor product smooths.
- p.type!=0 in summary.gam. - p.type!=0 in summary.gam.
1.8-13
* Added help file ?one.se.rule on the `one standard error rule' for obtaining
smoother models.
* bam(...,discrete=TRUE) no longer complains about more coefficients than data.
* 's', 'te', 'ti' and 't2' modified to allow user to specify that the smooth
should pass through zero at a specified point. See ?identifiability.
* anova.gam modified to use more appropriate reference degrees of freedom
for multiple model call, where possible. Also fixed to allow multiple
formulae models and to use -2*logLik in place of `deviance' for
general.family models.
* offsets allowed with multinomial, ziplss and gaulss families.
* gevlss family implementing generalized extreme value location, sclae and
shape models.
* Faster code used in 'uniquecombs'. Speeds up discretization step in
'bam(...,discrete=TRUE)'. Could still be improved for multi-column case.
* modification to 'smoothCon' to allow resetting of smooth supplied
constraints - enables fix of bug in bam handling of 't2' terms, where
parameterization of penalty and model matrix did not previously match
properly.
* clarification of `exclude' argument to predict.gam in help files.
* modification to 'plot.gam' etc, so that 'ylim' is no longer shifted by
'shift'.
* ylim and ... handling improved for 'fs' plot method (thanks Dave Miller)
* gam.check now recognises RStudio and plots appropriately.
* bam(...,sparse=TRUE) removed - not efficient, because of unavoidability
of dense off diagonal terms in X'X or equivalent. Deprecated since 1.8-5.
* tweak to initial.sp/g to avoid infinite loop in s.p. initialization, in
rather unusual circumstances. Thanks to Mark Bravington.
* bam and gam have `drop.intercept' argument to force the parametric terms not
to include a constant in their span, even when there are factor variables.
* Fix in Vb.corr (2nd order edf correction) for fixed smoothing parameter case.
* added 'all.vars1' to enable terms like s(x$y) in model formulae.
* modification to gam.fit4 to ignore 'start' if it is immediately worse than
'null.coef'.
* cSplineDes can now accept a 'derivs' argument.
* added drop.intercept handling for multiple formulae (mod by Matteo Fasiolo).
* 'gam.side' fix to avoid duplication of smooth model matrices to be tested
against, when checking for numerical degeneracy. Problem could occasionally
cause a failure (especially with bam), when the total matrix to be tested
against ended upo with more columns than rows.
* 4 occurances of as.name("model.frame") changed to quote(stats::model.frame)
* fix in predict.bamd discrete prediction code to be a bit more relaxed about
use of as.factor, etc in formulae.
* fix in predict.gam handling of 'na.action' to avoid trying to get type of
na.action from name of na.action function - this was fragile to nesting
and could cause predict.bam to fail in the presence of NA's.
* fix of gam.outer so that general families (e.g. cox.ph) can have all
their smoothing parameters supplied without then ignoring the penalties!
* fix in multiple formulae handling of fixed smoothing parameters.
* Fix of bug in zlim handling in vis.gam perspective plot with standard
errors. Thanks Petra Kuhnert.
* probit link added to 'jagam' (suggested by Kenneth Knoblauch).
* 'Sl.' routines revised to allow operation with non-linearly parameterized
smoothers.
* bug fix in Hessian computation in gam.fit5 - leading diagonal of Hessian of
log|Hp| could be wrong where Hp is penalized Hessian.
* better use of crossprod in gamlss.gH
1.8-12 1.8-12
......
Package: mgcv Package: mgcv
Version: 1.8-12 Version: 1.8-13
Author: Simon Wood <simon.wood@r-project.org> Author: Simon Wood <simon.wood@r-project.org>
Maintainer: 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 Title: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness
...@@ -16,6 +16,6 @@ LazyLoad: yes ...@@ -16,6 +16,6 @@ LazyLoad: yes
ByteCompile: yes ByteCompile: yes
License: GPL (>= 2) License: GPL (>= 2)
NeedsCompilation: yes NeedsCompilation: yes
Packaged: 2016-03-02 18:21:25 UTC; sw283 Packaged: 2016-07-20 08:12:21 UTC; sw283
Repository: CRAN Repository: CRAN
Date/Publication: 2016-03-03 07:57:37 Date/Publication: 2016-07-21 19:18:06
This diff is collapsed.
...@@ -9,7 +9,7 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity, ...@@ -9,7 +9,7 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
gam2objective, gam2objective,
gamm, gam.check, gam.control,gam.fit3, gamm, gam.check, gam.control,gam.fit3,
gam.fit, gam.outer,gam.vcomp, gamSim , gam.fit, gam.outer,gam.vcomp, gamSim ,
gaulss,gam.side,get.var, gaulss,gam.side,get.var,gevlss,
influence.gam, influence.gam,
in.out,inSide,interpret.gam,initial.sp, in.out,inSide,interpret.gam,initial.sp,
jagam, jagam,
...@@ -44,7 +44,7 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity, ...@@ -44,7 +44,7 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
residuals.gam,rig,rTweedie,rmvn, residuals.gam,rig,rTweedie,rmvn,
Rrank,s,scat,sdiag,"sdiag<-", Rrank,s,scat,sdiag,"sdiag<-",
sim2jam, sim2jam,
slanczos, slanczos,smooth2random,
smoothCon,smooth.construct,smooth.construct2, smoothCon,smooth.construct,smooth.construct2,
smooth.construct.bs.smooth.spec, smooth.construct.bs.smooth.spec,
smooth.construct.cc.smooth.spec, smooth.construct.cc.smooth.spec,
......
This diff is collapsed.
...@@ -33,7 +33,9 @@ cox.ph <- function (link = "identity") { ...@@ -33,7 +33,9 @@ cox.ph <- function (link = "identity") {
y.order <- order(G$y,decreasing=TRUE) y.order <- order(G$y,decreasing=TRUE)
G$family$data$y.order <- y.order G$family$data$y.order <- y.order
G$y <- G$y[y.order] G$y <- G$y[y.order]
attrX <- attributes(G$X)
G$X <- G$X[y.order,,drop=FALSE] G$X <- G$X[y.order,,drop=FALSE]
attributes(G$X) <- attrX
G$w <- G$w[y.order] G$w <- G$w[y.order]
}) })
...@@ -108,7 +110,7 @@ cox.ph <- function (link = "identity") { ...@@ -108,7 +110,7 @@ cox.ph <- function (link = "identity") {
rd <- qf <- NULL ## these functions currently undefined for Cox PH rd <- qf <- NULL ## these functions currently undefined for Cox PH
ll <- function(y,X,coef,wt,family,deriv=0,d1b=0,d2b=0,Hp=NULL,rank=0,fh=NULL,D=NULL) { ll <- function(y,X,coef,wt,family,offset=NULL,deriv=0,d1b=0,d2b=0,Hp=NULL,rank=0,fh=NULL,D=NULL) {
## function defining the cox model log lik. ## function defining the cox model log lik.
## Calls C code "coxlpl" ## Calls C code "coxlpl"
## deriv codes: 0 - evaluate the log likelihood ## deriv codes: 0 - evaluate the log likelihood
...@@ -122,7 +124,7 @@ cox.ph <- function (link = "identity") { ...@@ -122,7 +124,7 @@ cox.ph <- function (link = "identity") {
## or its Choleski factor ## or its Choleski factor
## D is the diagonal pre-conditioning matrix used to obtain Hp ## D is the diagonal pre-conditioning matrix used to obtain Hp
## if Hr is the raw Hp then Hp = D*t(D*Hr) ## if Hr is the raw Hp then Hp = D*t(D*Hr)
if (!is.null(offset)&&sum(offset!=0)) stop("cox.ph does not yet handle offsets")
##tr <- sort(unique(y),decreasing=TRUE) ##tr <- sort(unique(y),decreasing=TRUE)
tr <- unique(y) tr <- unique(y)
r <- match(y,tr) r <- match(y,tr)
......
...@@ -795,7 +795,7 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) { ...@@ -795,7 +795,7 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) {
y1 <- y + (y == 0) y1 <- y + (y == 0)
theta <- if (p == 1) log(y1/mu) else (y1^(1 - p) - mu^(1 - p))/(1 - p) theta <- if (p == 1) log(y1/mu) else (y1^(1 - p) - mu^(1 - p))/(1 - p)
kappa <- if (p == 2) log(y1/mu) else (y^(2 - p) - mu^(2 - p))/(2 - p) kappa <- if (p == 2) log(y1/mu) else (y^(2 - p) - mu^(2 - p))/(2 - p)
2 * (y * theta - kappa) * wt pmax(2 * (y * theta - kappa) * wt,0)
} }
Dd <- function(y, mu, theta, wt, level=0) { Dd <- function(y, mu, theta, wt, level=0) {
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
...@@ -427,7 +427,7 @@ smooth2random.t2.smooth <- function(object,vnames,type=1) { ...@@ -427,7 +427,7 @@ smooth2random.t2.smooth <- function(object,vnames,type=1) {
## 2. rind: and index vector such that if br is the vector of ## 2. rind: and index vector such that if br is the vector of
## random coefficients for the term, br[rind] is the coefs in ## random coefficients for the term, br[rind] is the coefs in
## order for this term. rinc - dummy here. ## order for this term. rinc - dummy here.
## 3. A matrix, U, that transforms coefs, in order [rand1, rand2,... fix] ## 3. A matrix, trans.D, that transforms coefs, in order [rand1, rand2,... fix]
## back to original parameterization. If null, then not needed. ## back to original parameterization. If null, then not needed.
## 4. A matrix Xf for the fixed effects, if any. ## 4. A matrix Xf for the fixed effects, if any.
## 5. fixed TRUE/FALSE if its fixed or not. If fixed the other stuff is ## 5. fixed TRUE/FALSE if its fixed or not. If fixed the other stuff is
...@@ -483,7 +483,7 @@ smooth2random.mgcv.smooth <- function(object,vnames,type=1) { ...@@ -483,7 +483,7 @@ smooth2random.mgcv.smooth <- function(object,vnames,type=1) {
## Returns 1. a list of random effects, including grouping factors, and ## Returns 1. a list of random effects, including grouping factors, and
## a fixed effects matrix. Grouping factors, model matrix and ## a fixed effects matrix. Grouping factors, model matrix and
## model matrix name attached as attributes, to each element. ## model matrix name attached as attributes, to each element.
## 2. rind: and index vector such that if br is the vector of ## 2. rind: an index vector such that if br is the vector of
## random coefficients for the term, br[rind] is the coefs in ## random coefficients for the term, br[rind] is the coefs in
## order for this term. rinc - dummy here. ## order for this term. rinc - dummy here.
## 3. A matrix, U, + vec D that transforms coefs, in order [rand1, rand2,... fix] ## 3. A matrix, U, + vec D that transforms coefs, in order [rand1, rand2,... fix]
...@@ -546,7 +546,7 @@ smooth2random.tensor.smooth <- function(object,vnames,type=1) { ...@@ -546,7 +546,7 @@ smooth2random.tensor.smooth <- function(object,vnames,type=1) {
## Returns 1. a list of random effects, including grouping factors, and ## Returns 1. a list of random effects, including grouping factors, and
## a fixed effects matrix. Grouping factors, model matrix and ## a fixed effects matrix. Grouping factors, model matrix and
## model matrix name attached as attributes, to each element. ## model matrix name attached as attributes, to each element.
## 2. rind: and index vector such that if br is the vector of ## 2. rind: an index vector such that if br is the vector of
## random coefficients for the term, br[rind] is the coefs in ## random coefficients for the term, br[rind] is the coefs in
## order for this term. rinc - dummy here. ## order for this term. rinc - dummy here.
## 3. A matrix, U, that transforms coefs, in order [rand1, rand2,... fix] ## 3. A matrix, U, that transforms coefs, in order [rand1, rand2,... fix]
...@@ -1201,7 +1201,10 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis ...@@ -1201,7 +1201,10 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
# random terms. correlation describes the correlation structure. This routine is basically an interface # random terms. correlation describes the correlation structure. This routine is basically an interface
# between the basis constructors provided in mgcv and the gammPQL routine used to estimate the model. # between the basis constructors provided in mgcv and the gammPQL routine used to estimate the model.
{ if (inherits(family,"extended.family")) warning("family are not designed for use with gamm!") { if (inherits(family,"extended.family")) warning("family are not designed for use with gamm!")
## lmeControl turns sigma=NULL into sigma=0, but if you supply sigma=0 rejects it,
## which will break the standard handling of the control list. Following line fixes.
## but actually Martin M has now fixed lmeControl...
##if (!is.null(control$sigma)&&control$sigma==0) control$sigma <- NULL
control <- do.call("lmeControl",control) control <- do.call("lmeControl",control)
# check that random is a named list # check that random is a named list
if (!is.null(random)) if (!is.null(random))
...@@ -1248,7 +1251,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis ...@@ -1248,7 +1251,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
mf$min.sp <- mf$H <- mf$gamma <- mf$fit <- mf$niterPQL <- mf$verbosePQL <- mf$G <- mf$method <- mf$... <- NULL mf$min.sp <- mf$H <- mf$gamma <- mf$fit <- mf$niterPQL <- mf$verbosePQL <- mf$G <- mf$method <- mf$... <- NULL
} }
mf$drop.unused.levels <- drop.unused.levels mf$drop.unused.levels <- drop.unused.levels
mf[[1]] <- as.name("model.frame") mf[[1]] <- quote(stats::model.frame) ## as.name("model.frame")
pmf <- mf pmf <- mf
gmf <- eval(mf, parent.frame()) # the model frame now contains all the data, for the gam part only gmf <- eval(mf, parent.frame()) # the model frame now contains all the data, for the gam part only
gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly
...@@ -1268,7 +1271,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis ...@@ -1268,7 +1271,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
## summarize the *raw* input variables ## summarize the *raw* input variables
## note can't use get_all_vars here -- buggy with matrices ## note can't use get_all_vars here -- buggy with matrices
vars <- all.vars(gp$fake.formula[-2]) ## drop response here vars <- all.vars1(gp$fake.formula[-2]) ## drop response here
inp <- parse(text = paste("list(", paste(vars, collapse = ","),")")) inp <- parse(text = paste("list(", paste(vars, collapse = ","),")"))
dl <- eval(inp, data, parent.frame()) dl <- eval(inp, data, parent.frame())
names(dl) <- vars ## list of all variables needed names(dl) <- vars ## list of all variables needed
......
...@@ -12,8 +12,8 @@ write.jagslp <- function(resp,family,file,use.weights,offset=FALSE) { ...@@ -12,8 +12,8 @@ write.jagslp <- function(resp,family,file,use.weights,offset=FALSE) {
## write the JAGS code for the linear predictor ## write the JAGS code for the linear predictor
## and response distribution. ## and response distribution.
iltab <- ## table of inverse link functions iltab <- ## table of inverse link functions
c("eta[i]","exp(eta[i])","ilogit(eta[i])","1/eta[i]","eta[i]^2") c("eta[i]","exp(eta[i])","ilogit(eta[i])","phi(eta[i])","1/eta[i]","eta[i]^2")
names(iltab) <- c("identity","log","logit","inverse","sqrt") names(iltab) <- c("identity","log","logit","probit","inverse","sqrt")
if (!family$link%in%names(iltab)) stop("sorry link not yet handled") if (!family$link%in%names(iltab)) stop("sorry link not yet handled")
## code linear predictor and expected response... ## code linear predictor and expected response...
...@@ -114,7 +114,7 @@ sp.prior = "gamma",diagonalize=FALSE) { ...@@ -114,7 +114,7 @@ sp.prior = "gamma",diagonalize=FALSE) {
mf$family <- mf$knots <- mf$sp <- mf$file <- mf$control <- mf$family <- mf$knots <- mf$sp <- mf$file <- mf$control <-
mf$centred <- mf$sp.prior <- mf$diagonalize <- NULL mf$centred <- mf$sp.prior <- mf$diagonalize <- NULL
mf$drop.unused.levels <- drop.unused.levels mf$drop.unused.levels <- drop.unused.levels
mf[[1]]<-as.name("model.frame") mf[[1]] <- quote(stats::model.frame) ##as.name("model.frame")
pmf <- mf pmf <- mf
pmf$formula <- gp$pf pmf$formula <- gp$pf
...@@ -165,7 +165,7 @@ sp.prior = "gamma",diagonalize=FALSE) { ...@@ -165,7 +165,7 @@ sp.prior = "gamma",diagonalize=FALSE) {
## get initial values, for use by JAGS, and to guess suitable values for ## get initial values, for use by JAGS, and to guess suitable values for
## uninformative priors... ## uninformative priors...
lambda <- initial.spg(G$X,G$y,G$w,family,G$S,G$off,G$L) ## initial sp values lambda <- initial.spg(G$X,G$y,G$w,family,G$S,G$off,offset=G$offset,L=G$L) ## initial sp values
jags.ini <- list() jags.ini <- list()
lam <- if (is.null(G$L)) lambda else G$L%*%lambda lam <- if (is.null(G$L)) lambda else G$L%*%lambda
jin <- jini(G,lam) jin <- jini(G,lam)
......
This diff is collapsed.
...@@ -96,11 +96,13 @@ mvn <- function(d=2) { ...@@ -96,11 +96,13 @@ mvn <- function(d=2) {
nbeta <- ncol(G$X) nbeta <- ncol(G$X)
ntheta <- ydim*(ydim+1)/2 ## number of cov matrix factor params ntheta <- ydim*(ydim+1)/2 ## number of cov matrix factor params
lpi <- attr(G$X,"lpi") lpi <- attr(G$X,"lpi")
#offs <- attr(G$X,"offset")
XX <- crossprod(G$X) XX <- crossprod(G$X)
G$X <- cbind(G$X,matrix(0,nrow(G$X),ntheta)) ## add dummy columns to G$X G$X <- cbind(G$X,matrix(0,nrow(G$X),ntheta)) ## add dummy columns to G$X
#G$cmX <- c(G$cmX,rep(0,ntheta)) ## and corresponding column means #G$cmX <- c(G$cmX,rep(0,ntheta)) ## and corresponding column means
G$term.names <- c(G$term.names,paste("R",1:ntheta,sep=".")) G$term.names <- c(G$term.names,paste("R",1:ntheta,sep="."))
attr(G$X,"lpi") <- lpi attr(G$X,"lpi") <- lpi
#offs -> attr(G$X,"offset")
attr(G$X,"XX") <- XX attr(G$X,"XX") <- XX
## pad out sqrt of balanced penalty matrix to account for extra params ## pad out sqrt of balanced penalty matrix to account for extra params
attr(G$Sl,"E") <- cbind(attr(G$Sl,"E"),matrix(0,nbeta,ntheta)) attr(G$Sl,"E") <- cbind(attr(G$Sl,"E"),matrix(0,nbeta,ntheta))
...@@ -157,7 +159,7 @@ mvn <- function(d=2) { ...@@ -157,7 +159,7 @@ mvn <- function(d=2) {
##rd <- qf <- NULL ## these functions currently undefined for ##rd <- qf <- NULL ## these functions currently undefined for
ll <- function(y,X,coef,wt,family,deriv=0,d1b=NULL,d2b=NULL,Hp=NULL,rank=0,fh=NULL,D=NULL) { ll <- function(y,X,coef,wt,family,offset=NULL,deriv=0,d1b=NULL,d2b=NULL,Hp=NULL,rank=0,fh=NULL,D=NULL) {
## function defining the Multivariate Normal model log lik. ## function defining the Multivariate Normal model log lik.
## Calls C code "mvn_ll" ## Calls C code "mvn_ll"
## deriv codes: 0 - evaluate the log likelihood ## deriv codes: 0 - evaluate the log likelihood
...@@ -171,6 +173,9 @@ mvn <- function(d=2) { ...@@ -171,6 +173,9 @@ mvn <- function(d=2) {
## or its Choleski factor ## or its Choleski factor
## D is the diagonal pre-conditioning matrix used to obtain Hp ## D is the diagonal pre-conditioning matrix used to obtain Hp
## if Hr is the raw Hp then Hp = D*t(D*Hr) ## if Hr is the raw Hp then Hp = D*t(D*Hr)
if (!is.null(offset)) {
for (i in 1:length(offset)) if (sum(offset[[i]]!=0)) stop("mvn does not yet handle offsets")
}
lpi <- attr(X,"lpi") ## lpi[[k]] is index of model matrix columns for dim k lpi <- attr(X,"lpi") ## lpi[[k]] is index of model matrix columns for dim k
overlap <- attr(lpi,"overlap") ## do dimensions share terms? overlap <- attr(lpi,"overlap") ## do dimensions share terms?
drop <- attr(X,"drop") drop <- attr(X,"drop")
......
...@@ -267,7 +267,7 @@ gam.check <- function(b, old.style=FALSE, ...@@ -267,7 +267,7 @@ gam.check <- function(b, old.style=FALSE,
napredict(b$na.action, b$linear.predictors[,1]) else napredict(b$na.action, b$linear.predictors[,1]) else
napredict(b$na.action, b$linear.predictors) napredict(b$na.action, b$linear.predictors)
## if (b$method%in%c("GCV","GACV","UBRE","REML","ML","P-ML","P-REML","mle.REML","mle.ML","PQL")) { ## if (b$method%in%c("GCV","GACV","UBRE","REML","ML","P-ML","P-REML","mle.REML","mle.ML","PQL")) {
old.par<-par(mfrow=c(2,2)) if (is.null(.Platform$GUI) || .Platform$GUI != "RStudio") old.par <- par(mfrow=c(2,2))
if (old.style) if (old.style)
qqnorm(resid,...) qqnorm(resid,...)
else else
...@@ -326,7 +326,7 @@ gam.check <- function(b, old.style=FALSE, ...@@ -326,7 +326,7 @@ gam.check <- function(b, old.style=FALSE,
cat("indicate that k is too low, especially if edf is close to k\'.\n\n") cat("indicate that k is too low, especially if edf is close to k\'.\n\n")
printCoefmat(kchck,digits=3); printCoefmat(kchck,digits=3);
} }
par(old.par) if (is.null(.Platform$GUI) ||.Platform$GUI != "RStudio") par(old.par)
## } else plot(linpred,resid,xlab="linear predictor",ylab="residuals",...) ## } else plot(linpred,resid,xlab="linear predictor",ylab="residuals",...)
} ## end of gam.check } ## end of gam.check
...@@ -675,7 +675,6 @@ plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult= ...@@ -675,7 +675,6 @@ plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=
fac <- rep(x$flev,rep(n,nf)) fac <- rep(x$flev,rep(n,nf))
dat <- data.frame(fac,xx) dat <- data.frame(fac,xx)
names(dat) <- c(x$fterm,x$base$term) names(dat) <- c(x$fterm,x$base$term)
# X <- Predict.matrix.fs.interaction(x,dat)
X <- PredictMat(x,dat) X <- PredictMat(x,dat)
if (is.null(xlab)) xlabel <- x$base$term else xlabel <- xlab if (is.null(xlab)) xlabel <- x$base$term else xlabel <- xlab
if (is.null(ylab)) ylabel <- label else ylabel <- ylab if (is.null(ylab)) ylabel <- label else ylabel <- ylab
...@@ -683,7 +682,8 @@ plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult= ...@@ -683,7 +682,8 @@ plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=
main="",x=xx,n=n,nf=nf)) main="",x=xx,n=n,nf=nf))
} else { ## produce the plot } else { ## produce the plot
ind <- 1:P$n ind <- 1:P$n
plot(P$x[ind],P$fit[ind],ylim=range(P$fit),xlab=P$xlab,ylab=P$ylab,type="l") if(is.null(ylim)) ylim <- range(P$fit)
plot(P$x[ind],P$fit[ind],ylim=ylim,xlab=P$xlab,ylab=P$ylab,type="l",...)
if (P$nf>1) for (i in 2:P$nf) { if (P$nf>1) for (i in 2:P$nf) {
ind <- ind + P$n ind <- ind + P$n
if (scheme==0) lines(P$x,P$fit[ind],lty=i,col=i) else if (scheme==0) lines(P$x,P$fit[ind],lty=i,col=i) else
...@@ -867,17 +867,17 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2, ...@@ -867,17 +867,17 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (min.r < ylimit[1]) ylimit[1] <- min.r if (min.r < ylimit[1]) ylimit[1] <- min.r
} }
} }
if (!is.null(ylim)) ylimit <- ylim ylimit <- if (is.null(ylim)) ylimit <- ylimit + shift else ylim
## plot the smooth... ## plot the smooth...
if (shade) { if (shade) {
plot(P$x,trans(P$fit+shift),type="n",xlab=P$xlab,ylim=trans(ylimit+shift), plot(P$x,trans(P$fit+shift),type="n",xlab=P$xlab,ylim=trans(ylimit),
xlim=P$xlim,ylab=P$ylab,main=P$main,...) xlim=P$xlim,ylab=P$ylab,main=P$main,...)
polygon(c(P$x,P$x[n:1],P$x[1]), polygon(c(P$x,P$x[n:1],P$x[1]),
trans(c(ul,ll[n:1],ul[1])+shift),col = shade.col,border = NA) trans(c(ul,ll[n:1],ul[1])+shift),col = shade.col,border = NA)
lines(P$x,trans(P$fit+shift),...) lines(P$x,trans(P$fit+shift),...)
} else { ## ordinary plot } else { ## ordinary plot
plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,ylim=trans(ylimit+shift),xlim=P$xlim, plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,ylim=trans(ylimit),xlim=P$xlim,
ylab=P$ylab,main=P$main,...) ylab=P$ylab,main=P$main,...)
if (is.null(list(...)[["lty"]])) { if (is.null(list(...)[["lty"]])) {
lines(P$x,trans(ul+shift),lty=2,...) lines(P$x,trans(ul+shift),lty=2,...)
...@@ -936,9 +936,10 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2, ...@@ -936,9 +936,10 @@ plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
if (scale==0&&is.null(ylim)) { if (scale==0&&is.null(ylim)) {
if (partial.resids) ylimit <- range(P$p.resid,na.rm=TRUE) else ylimit <-range(P$fit) if (partial.resids) ylimit <- range(P$p.resid,na.rm=TRUE) else ylimit <-range(P$fit)
} }
if (!is.null(ylim)) ylimit <- ylim ylimit <- if (is.null(ylim)) ylimit <- ylimit + shift else ylim
plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab, plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,
ylab=P$ylab,ylim=trans(ylimit+shift),xlim=P$xlim,main=P$main,...) ylab=P$ylab,ylim=trans(ylimit),xlim=P$xlim,main=P$main,...)
if (rug) { if (rug) {
if (jit) rug(jitter(as.numeric(P$raw)),...) if (jit) rug(jitter(as.numeric(P$raw)),...)
else rug(as.numeric(P$raw),...) else rug(as.numeric(P$raw),...)
...@@ -1442,15 +1443,15 @@ vis.gam <- function(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,color="he ...@@ -1442,15 +1443,15 @@ vis.gam <- function(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,color="he
lo.col <- "green" lo.col <- "green"
hi.col <- "red" hi.col <- "red"
} }
if (!is.null(zlim)) if (!is.null(zlim)) {
{ if (length(zlim)!=2||zlim[1]>=zlim[2]) stop("Something wrong with zlim") if (length(zlim)!=2||zlim[1]>=zlim[2]) stop("Something wrong with zlim")
min.z<-zlim[1] min.z<-zlim[1]
max.z<-zlim[2] max.z<-zlim[2]
} else } else {
{ z.max<-max(fv$fit+fv$se.fit*se,na.rm=TRUE) max.z<-max(fv$fit+fv$se.fit*se,na.rm=TRUE)
z.min<-min(fv$fit-fv$se.fit*se,na.rm=TRUE) min.z<-min(fv$fit-fv$se.fit*se,na.rm=TRUE)
zlim<-c(min.z,max.z)
} }
zlim<-c(z.min,z.max)
z<-fv$fit-fv$se.fit*se;z<-matrix(z,n.grid,n.grid) z<-fv$fit-fv$se.fit*se;z<-matrix(z,n.grid,n.grid)
if (plot.type=="contour") warning("sorry no option for contouring with errors: try plot.gam") if (plot.type=="contour") warning("sorry no option for contouring with errors: try plot.gam")
......
This diff is collapsed.
mgcv (1.8-13-1) unstable; urgency=medium
* New upstream release
-- Dirk Eddelbuettel <edd@debian.org> Thu, 21 Jul 2016 20:09:11 -0500
mgcv (1.8-12-2) unstable; urgency=medium mgcv (1.8-12-2) unstable; urgency=medium
* debian/compat: Created (Closes: #828882) * debian/compat: Created (Closes: #828882)
......
...@@ -74,7 +74,10 @@ Tweedie, M. C. K. (1984). An index which distinguishes between ...@@ -74,7 +74,10 @@ Tweedie, M. C. K. (1984). An index which distinguishes between
Golden Jubilee International Conference (Eds. J. K. Ghosh and J. Golden Jubilee International Conference (Eds. J. K. Ghosh and J.
Roy), pp. 579-604. Calcutta: Indian Statistical Institute. Roy), pp. 579-604. Calcutta: Indian Statistical Institute.
Wood, S.N., N. Pya and B. Saefken (2015), Smoothing parameter and model selection for general smooth models. \url{http://arxiv.org/abs/1511.03864} Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and
model selection for general smooth models.
Journal of the American Statistical Association.
\url{http://arxiv.org/abs/1511.03864}
} }
\seealso{\code{\link{ldTweedie}}, \code{\link{rTweedie}}} \seealso{\code{\link{ldTweedie}}, \code{\link{rTweedie}}}
......
...@@ -7,10 +7,8 @@ ...@@ -7,10 +7,8 @@
\code{gam} objects. For a single fitted \code{gam} object, Wald tests of \code{gam} objects. For a single fitted \code{gam} object, Wald tests of
the significance of each parametric and smooth term are performed, so interpretation the significance of each parametric and smooth term are performed, so interpretation
is analogous to \code{\link{drop1}} rather than \code{anova.lm} (i.e. it's like type III ANOVA, is analogous to \code{\link{drop1}} rather than \code{anova.lm} (i.e. it's like type III ANOVA,
rather than a sequential type I ANOVA). Otherwise rather than a sequential type I ANOVA). Otherwise the fitted models are compared using an analysis of deviance table: this latter approach should not be use to test the significance of terms which can be penalized
the fitted models are compared using an analysis of deviance table: this latter approach to zero. Models to be compared should be fitted to the same data using the same smoothing parameter selection method.
should not be use to test the significance of terms which can be penalized
to zero. See details.
} }
\usage{ \usage{
\method{anova}{gam}(object, ..., dispersion = NULL, test = NULL, \method{anova}{gam}(object, ..., dispersion = NULL, test = NULL,
...@@ -32,7 +30,8 @@ p-values. See \code{\link{summary.gam}} for details.} ...@@ -32,7 +30,8 @@ p-values. See \code{\link{summary.gam}} for details.}
} }
\details{ If more than one fitted model is provided than \code{anova.glm} is \details{ If more than one fitted model is provided than \code{anova.glm} is
used, with the difference in model degrees of freedom being taken as the difference used, with the difference in model degrees of freedom being taken as the difference
in effective degress of freedom. The p-values resulting from this are only approximate, in effective degress of freedom (when possible this is a smoothing parameter uncertainty corrected version).
The p-values resulting from this are only approximate,
and must be used with care. The approximation is most accurate when the comparison and must be used with care. The approximation is most accurate when the comparison
relates to unpenalized terms, or smoothers with a null space of dimension greater than zero. relates to unpenalized terms, or smoothers with a null space of dimension greater than zero.
(Basically we require that the difference terms could be well approximated by unpenalized (Basically we require that the difference terms could be well approximated by unpenalized
...@@ -82,7 +81,7 @@ Wood, S.N. (2013b) A simple test for random effects in regression models. Biomet ...@@ -82,7 +81,7 @@ Wood, S.N. (2013b) A simple test for random effects in regression models. Biomet
\author{ Simon N. Wood \email{simon.wood@r-project.org} with substantial \author{ Simon N. Wood \email{simon.wood@r-project.org} with substantial
improvements by Henric Nilsson.} improvements by Henric Nilsson.}
\section{WARNING}{ If models 'a' and 'b' differ only in terms with no un-penalized components then \section{WARNING}{ If models 'a' and 'b' differ only in terms with no un-penalized components (such as random effects) then
p values from anova(a,b) are unreliable, and usually much too low. p values from anova(a,b) are unreliable, and usually much too low.
Default P-values will usually be wrong for parametric terms penalized using `paraPen': use freq=TRUE Default P-values will usually be wrong for parametric terms penalized using `paraPen': use freq=TRUE
......
...@@ -21,8 +21,8 @@ bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL, ...@@ -21,8 +21,8 @@ bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
na.action=na.omit, offset=NULL,method="fREML",control=list(), na.action=na.omit, offset=NULL,method="fREML",control=list(),
select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL, select=FALSE,scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL,
paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE, paraPen=NULL,chunk.size=10000,rho=0,AR.start=NULL,discrete=FALSE,
sparse=FALSE,cluster=NULL,nthreads=NA,gc.level=1,use.chol=FALSE, cluster=NULL,nthreads=NA,gc.level=1,use.chol=FALSE,
samfrac=1,drop.unused.levels=TRUE,G=NULL,fit=TRUE,...) samfrac=1,drop.unused.levels=TRUE,G=NULL,fit=TRUE,drop.intercept=NULL,...)
} }
%- maybe also `usage' for other objects documented here. %- maybe also `usage' for other objects documented here.
...@@ -129,10 +129,6 @@ there are no breaks in AR1 correlaion.} ...@@ -129,10 +129,6 @@ there are no breaks in AR1 correlaion.}
\item{discrete}{with \code{method="fREML"} it is possible to discretize covariates for storage and efficiency reasons. \item{discrete}{with \code{method="fREML"} it is possible to discretize covariates for storage and efficiency reasons.
If \code{discrete} is \code{TRUE}, a number or a vector of numbers for each smoother term, then discretization happens. If numbers are supplied they give the number of discretization bins. Experimental at present. } If \code{discrete} is \code{TRUE}, a number or a vector of numbers for each smoother term, then discretization happens. If numbers are supplied they give the number of discretization bins. Experimental at present. }
\item{sparse}{Deprecated. If all smooths are P-splines and all tensor products are of the form \code{te(...,bs="ps",np=FALSE)}
then in principle computation could be made faster using sparse matrix methods, and you could set this to \code{TRUE}.
In practice the speed up is disappointing, and the computation is less well conditioned than the default. See details.}
\item{cluster}{\code{bam} can compute the computationally dominant QR decomposition in parallel using \link[parallel]{parLapply} \item{cluster}{\code{bam} can compute the computationally dominant QR decomposition in parallel using \link[parallel]{parLapply}
from the \code{parallel} package, if it is supplied with a cluster on which to do this (a cluster here can be some cores of a from the \code{parallel} package, if it is supplied with a cluster on which to do this (a cluster here can be some cores of a
single machine). See details and example code. single machine). See details and example code.
...@@ -160,6 +156,9 @@ involving factor variables you might want to turn this off. Only do so if you kn ...@@ -160,6 +156,9 @@ involving factor variables you might want to turn this off. Only do so if you kn
\item{fit}{if \code{FALSE} then the model is set up for fitting but not estimated, and an object is returned, suitable for passing as the \code{G} argument to \code{bam}.} \item{fit}{if \code{FALSE} then the model is set up for fitting but not estimated, and an object is returned, suitable for passing as the \code{G} argument to \code{bam}.}
\item{drop.intercept}{Set to \code{TRUE} to force the model to really not have the a constant in the parametric model part,
even with factor variables present.}
\item{...}{further arguments for \item{...}{further arguments for
passing on e.g. to \code{gam.fit} (such as \code{mustart}). } passing on e.g. to \code{gam.fit} (such as \code{mustart}). }
...@@ -170,8 +169,7 @@ An object of class \code{"gam"} as described in \code{\link{gamObject}}. ...@@ -170,8 +169,7 @@ An object of class \code{"gam"} as described in \code{\link{gamObject}}.
} }
\details{ When \code{discrete=FALSE}, \code{bam} operates by first setting up the basis characteristics for the smooths, using a representative subsample \details{ When \code{discrete=FALSE}, \code{bam} operates by first setting up the basis characteristics for the smooths, using a representative subsample of the data. Then the model matrix is constructed in blocks using \code{\link{predict.gam}}. For each block the factor R,
of the data. Then the model matrix is constructed in blocks using \code{\link{predict.gam}}. For each block the factor R,
from the QR decomposition of the whole model matrix is updated, along with Q'y. and the sum of squares of y. At the end of from the QR decomposition of the whole model matrix is updated, along with Q'y. and the sum of squares of y. At the end of
block processing, fitting takes place, without the need to ever form the whole model matrix. block processing, fitting takes place, without the need to ever form the whole model matrix.
...@@ -192,21 +190,14 @@ If \code{cluster} is provided as a cluster set up using \code{\link[parallel]{ma ...@@ -192,21 +190,14 @@ If \code{cluster} is provided as a cluster set up using \code{\link[parallel]{ma
When \code{discrete=TRUE} the covariate data are first discretized. Discretization takes place on a smooth by smooth basis, or in the case of tensor product smooths (or any smooth that can be represented as such, such as random effects), separately for each marginal smooth. The required spline bases are then evaluated at the discrete values, and stored, along with index vectors indicating which original observation they relate to. Fitting is by a version of performance oriented iteration/PQL using REML smoothing parameter selection on each iterative working model (as for the default method). The iteration is based on the derivatives of the REML score, without computing the score itself, allowing the expensive computations to be reduced to one parallel block Cholesky decomposition per iteration (plus two basic operations of equal cost, but easily parallelized). Unlike standard POI/PQL, only one step of the smoothing parameter update for the working model is taken at each step (rather than iterating to the optimal set of smoothing parameters for each working model). At each step a weighted model matrix crossproduct of the model matrix is required - this is efficiently computed from the pre-computed basis functions evaluated at the discretized covariate values. Efficient computation with tensor product terms means that some terms within a tensor product may be re-ordered for maximum efficiency. Parallel computation is controlled using the \code{nthreads} argument. For this method no cluster computation is used, and the \code{parallel} package is not required. When \code{discrete=TRUE} the covariate data are first discretized. Discretization takes place on a smooth by smooth basis, or in the case of tensor product smooths (or any smooth that can be represented as such, such as random effects), separately for each marginal smooth. The required spline bases are then evaluated at the discrete values, and stored, along with index vectors indicating which original observation they relate to. Fitting is by a version of performance oriented iteration/PQL using REML smoothing parameter selection on each iterative working model (as for the default method). The iteration is based on the derivatives of the REML score, without computing the score itself, allowing the expensive computations to be reduced to one parallel block Cholesky decomposition per iteration (plus two basic operations of equal cost, but easily parallelized). Unlike standard POI/PQL, only one step of the smoothing parameter update for the working model is taken at each step (rather than iterating to the optimal set of smoothing parameters for each working model). At each step a weighted model matrix crossproduct of the model matrix is required - this is efficiently computed from the pre-computed basis functions evaluated at the discretized covariate values. Efficient computation with tensor product terms means that some terms within a tensor product may be re-ordered for maximum efficiency. Parallel computation is controlled using the \code{nthreads} argument. For this method no cluster computation is used, and the \code{parallel} package is not required.
If the deprecated argument \code{sparse=TRUE} then QR updating is replaced by an alternative scheme, in which the model matrix is stored whole
as a sparse matrix. This only makes sense if all smooths are P-splines and all tensor products are of the
form \code{te(...,bs="ps",np=FALSE)}, but no check is made. The computations are then based on the Choleski decomposition of
the crossproduct of the sparse model matrix. Although this crossproduct is nearly dense, sparsity should make its
formation efficient, which is useful as it is the leading order term in the operations count. However there is no benefit
in using sparse methods to form the Choleski decomposition, given that the crossproduct is dense.
In practice the sparse matrix handling overheads mean that modest or no speed ups are produced
by this approach, while the computation is less stable than the default, and the memory footprint often higher
(but please let the author know if you find an example where the speedup is really worthwhile).
} }