Commit 8e9a69da authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.7-26

parent bd9ae805
** denotes quite substantial/important changes
*** denotes really big changes
1.7-26
* Namespace fixes.
1.7-25
* code added to allow openMP based multi-threading in gam fits (see
?gam.control and ?"mgcv-parallel").
* bam now allows AR1 error model to be split blockwise. See argument
'AR.start'.
* magic.post.proc made more efficient (one of two O(np^2) steps removed).
* var.summary now coerces character to factor.
* bugs fixed whereby etastart etc were not passed to initial.spg and
get.null.coefs. Thanks to Gavin Simpson.
* reformulate removed from predict.gam to avoid (slow) repeated parser
calls.
* gaussian(link="log") initialization fixed so that negative data
does not make it fail, via fix.family patching function.
* bug fix in plot method for "fs" basis - ignored any side conditions.
Thanks to Martijn Weiling and Jacolien van Rij.
* gamm now checks whether smooths nested in factors have illegal side
conditions, and halts if so (re-ordering formula can help).
* anova.glmlist no longer called.
* Compiled code now uses R_chck_calloc and R_chk_free for memory management
to avoid the possibility of unfriendly exit on running out of memory.
* fix in gam.side which would fail with unpenalized interactions in the
presence of main effects.
1.7-24
* Examples pruned in negbin, smooth.construct.ad.smooth.spec and bam help
......
Package: mgcv
Version: 1.7-24
Version: 1.7-26
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
estimation
Description: Routines for GAMs and other generalized ridge regression
with multiple smoothing parameter selection by GCV, REML or
UBRE/AIC. Also GAMMs. Includes a gam() function.
Description: Routines for GAMs and other generalized ridge regression
with multiple smoothing parameter selection by GCV, REML
or UBRE/AIC. Also GAMMs. Includes a gam() function.
Priority: recommended
Depends: R (>= 2.14.0), stats, graphics
Imports: nlme, methods, Matrix
Suggests: nlme (>= 3.1-64), splines, Matrix, parallel
Depends: R (>= 2.14.0), nlme (>= 3.1-64)
Imports: methods, stats, graphics, Matrix
Suggests: splines, parallel
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2013-06-05 07:28:12 UTC; sw283
Packaged: 2013-09-06 09:30:03 UTC; ripley
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2013-06-05 10:34:04
Date/Publication: 2013-09-06 11:31:13
This diff is collapsed.
c5fca2e3b064e4b18ad5c835193c7e19 *DESCRIPTION
5f405da54d89b1e236e272dad685cc2c *NAMESPACE
26423bcf418f00a9a89c06ddaa0646b5 *R/bam.r
2d679b4646af9013c3e6b11a61e171a6 *R/fast-REML.r
7873e574394156cd1ebb920c693bf3f0 *R/gam.fit3.r
d364a4f120922c8c4bd91fd70611687b *ChangeLog
82baebca68c3f9b9b0615e9b7daeca0b *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
4466d00a9e94717f775d55c8fead7f91 *NAMESPACE
5b5d847f492e7c5a687bafe83f241a8c *R/bam.r
6f5c9798cb6e42b3431febdd3149a26e *R/fast-REML.r
b65bae33a48ef86e7459192fb74cbee4 *R/gam.fit3.r
902657a0ee2dedc3fdfa501bf3b37c5b *R/gam.sim.r
6045d40194b509250344c55436845f8e *R/gamm.r
0421a1ed94ac15480a0a5a55e3886638 *R/mgcv.r
2370f3d746b680472cd51ccb9dc10284 *R/plots.r
7c562d025d05df4a5002cd0bd2c0030c *R/smooth.r
bc0821492c2db83bc43ecef857b455f9 *R/soap.r
fb66d6c18398411a99ffcb788b854f13 *R/sparse.r
01a1ab28361fc9a574887a74a9309edd *changeLog
d74688ecfa7912cd40f5373154da27d2 *R/gamm.r
dfda3f652d3d245a26a9109f480bbe4e *R/mgcv.r
2c8567dee7cc29736e260a9518361c6a *R/misc.r
f8b84b98c0524c2664885294650af3b8 *R/plots.r
b0b6e71c15742608f2dec02c1d8ae71b *R/smooth.r
90f72ef91dfc2a25c262fc002f781bf7 *R/soap.r
0c5d2acdc985d2af9de6bb736d4df2ab *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
88d77139cc983317b6acd8c5f1252ab9 *gnugpl2.txt
f693920e12f8a6f2b6cab93648628150 *index
9388a0979dab4fe1988b777ae4dcfb0a *inst/CITATION
5bf12ddc0dab9daae72271b96a15c539 *man/Predict.matrix.Rd
c51c9b8c9c73f81895176ded39b91394 *man/Predict.matrix.cr.smooth.Rd
e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
3e5e30b44947d3ddf00fcd55be219038 *man/Tweedie.Rd
27b7ae82bab920543222c4428d33cee1 *man/anova.gam.Rd
7c6fcc602c8278f96ef8cacd28a9c8ea *man/bam.Rd
b385d6d5419d0d6aefe03af1a79d5c4e *man/bam.update.Rd
7408512e0b41ba735bc79fce357282ad *man/bam.Rd
71bde8b8caa24a36529ce7e0ac3165d8 *man/bam.update.Rd
bb5ec26743d46d7ce1dbf128bceab35a *man/cSplineDes.Rd
9b4d616d1b6c4a46ca77d16cded3f806 *man/choose.k.Rd
c03748964ef606621418e428ae49b103 *man/columb.Rd
......@@ -37,7 +37,7 @@ d87ff6287d7e343ea24d2053f4724b35 *man/formula.gam.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
d7392fabfe137e327d89cc9686946740 *man/gam.Rd
fe61dd0efab9e920c17335faf3d5764c *man/gam.check.Rd
96c9417e4ac5d79ec9ed3f363adfc4e9 *man/gam.control.Rd
9ba1c43c08de6de712138a0741a94e68 *man/gam.control.Rd
44db24b66ce63bc16d2c8bc3f5b42ac5 *man/gam.convergence.Rd
58ab3b3d6f4fd0d008d73c3c4e6d3305 *man/gam.fit.Rd
21339a5d1eb8c83679dd9022ab682b5e *man/gam.fit3.Rd
......@@ -47,7 +47,7 @@ e969287d1a5c281faa7eb6cfce31a7c5 *man/gam.outer.Rd
956bf1a6ac1361dd0403c28153b03a9b *man/gam.side.Rd
b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
a66a814cc4c6f806e824751fda519ae0 *man/gam2objective.Rd
807a8e175d84ba3eac87f8ef2931f6e4 *man/gamObject.Rd
eacda367e277e56bc4cc58153dfa5318 *man/gamObject.Rd
0ac5fb78c9db628ce554a8f68588058c *man/gamSim.Rd
9c301060b3c35bc6ab83e6f38ce8d3b6 *man/gamm.Rd
39ae109127110152e97dc9f79e08bb14 *man/get.var.Rd
......@@ -60,10 +60,11 @@ aba56a0341ba9526a302e39d33aa9042 *man/interpret.gam.Rd
58e73ac26b93dc9d28bb27c8699e12cf *man/linear.functional.terms.Rd
5de18c3ad064a5bda4f9027d9455170a *man/logLik.gam.Rd
611f5f6acac9c5f40869c01cf7f75dd3 *man/ls.size.Rd
396548e4a40652c54e80751acbfb2e2c *man/magic.Rd
496388445d8cde9b8e0c3917cbe7461d *man/magic.post.proc.Rd
19a3f7da03e76d7d2fc0a4e544564c72 *man/magic.Rd
5169af4be5fccf9fa79b6de08e9ea035 *man/magic.post.proc.Rd
9427464e63d70da15fc468961ef0c29b *man/mgcv-FAQ.Rd
2db8cd353119cff94b57ff8ae06a6f5d *man/mgcv-package.Rd
8365171cd0280a0aa5ad90c30349c949 *man/mgcv-parallel.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
bc9b89db7e7ff246749551c16f5f1f07 *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
......@@ -123,20 +124,20 @@ d5a0f198090ecbfedaa6549f2918b997 *po/R-po.po
1a4a267ddcb87bb83f09c291d3e97523 *po/fr.po
813514ea4e046ecb4563eb3ae8aa202a *po/mgcv.pot
749ac663240fd6eaa2b725602d47ef2a *po/po.po
cd54024d76a9b53dc17ef26323fc053f *src/Makevars
94a2bcbb75cc60e8460e72ed154678c9 *src/gdi.c
c49d921d97e78dcd26088ec6326e6940 *src/Makevars
29b51c84b16b1692061d363065beac3b *src/gdi.c
49af97195accb65adc75620183d39a4c *src/general.h
622474c4de77d29db86e9a3f73dc7862 *src/init.c
1635205c9d5bb9a5713663775f47fd2a *src/magic.c
ad5838098ac3fe9d20443fcad3e9a62d *src/mat.c
de0ae24ea5cb533640a3ab57e0383595 *src/matrix.c
2a9d03aeee4de9228502620734b6c84c *src/init.c
e12a89a18a1aa847fe958b02f9dc9f80 *src/magic.c
0ed50431782a747557935ad347e4e64d *src/mat.c
6c9101f596440a09ecc3966a2374efd2 *src/matrix.c
0f8448f67d16668f9027084a2d9a1b52 *src/matrix.h
cb199160696be2d8550b11502542b951 *src/mgcv.c
2b38840c84b6f8d6c2dce7603df8b1e9 *src/mgcv.h
b9328a6d6aa86380fe52a971deccd3d4 *src/misc.c
7e0ba698a21a01150fda519661ef9857 *src/qp.c
96e43aa805215a0463fd13fd4b3d41f8 *src/mgcv.c
94660efe79de3659f210b023819afed0 *src/mgcv.h
c6b8c3093e2970654237c81ccfab1473 *src/misc.c
97bc53816f97fc68a37b2e0716a6dc0a *src/qp.c
cd563899be5b09897d1bf36a7889caa0 *src/qp.h
dcd47b55e6a2e5f21eac7ddbfb357fba *src/soap.c
e9cab4a461eb8e086a0e4834cbf16f30 *src/sparse-smooth.c
a0de1174b314d662859a6ba6f1d7f03b *src/tprs.c
15ed450da804a9dbd0833459af1f2474 *src/soap.c
449046a076c579cd0fafc137d250e65e *src/sparse-smooth.c
42f015e636b5fbd2a3a3bb46cbe24eaa *src/tprs.c
5352d5d2298acd9b03ee1895933d4fb4 *src/tprs.h
......@@ -69,11 +69,12 @@ importFrom(graphics,axis,box,contour,hist,lines,mtext, par, persp,plot,points,
# runif,termplot,terms,terms.formula, uniroot,
# var,vcov)
importFrom(stats,anova,influence,cooks.distance,logLik,vcov,residuals,predict,
model.matrix)
importFrom(nlme,Dim,corMatrix,logDet,pdConstruct,pdFactor,pdMatrix)
#importFrom(Matrix,t,mean,colMeans,colSums,chol,solve,diag)
import(Matrix)
importFrom(stats,anova,influence,cooks.distance,logLik,vcov,residuals,predict,model.matrix)
#importFrom(nlme,Dim,corMatrix,logDet,pdConstruct,pdFactor,pdMatrix,getGroupsFormula,lme,varFixed,lmeControl)
import(nlme)
importMethodsFrom(Matrix,t,colMeans,colSums,chol,solve,lu,expand)
importFrom(Matrix,Diagonal,sparseMatrix,Matrix)
#import(Matrix)
importFrom(methods,cbind2)
S3method(anova, gam)
......
This diff is collapsed.
......@@ -263,7 +263,7 @@ ldetS <- function(Sl,rho,fixed,np,root=FALSE) {
ind <- k.sp:(k.sp+m-1) ## index for smoothing parameters
## call gam.reparam to deal with this block
## in a stable way...
grp <- mgcv:::gam.reparam(Sl[[b]]$rS,lsp=rho[ind],deriv=2)
grp <- gam.reparam(Sl[[b]]$rS,lsp=rho[ind],deriv=2)
Sl[[b]]$lambda <- exp(rho[ind])
ldS <- ldS + grp$det
## next deal with the derivatives...
......@@ -661,7 +661,7 @@ ident.test <- function(X,E) {
## coeff vector, with dropped coefs re-nstated as zeroes.
Xnorm <- norm(X,type="F")
qrx <- qr(rbind(X/Xnorm,E),LAPACK=TRUE) ## pivoted QR
rank <- mgcv:::Rrank(qr.R(qrx),tol=.Machine$double.eps^.75)
rank <- Rrank(qr.R(qrx),tol=.Machine$double.eps^.75)
drop <- qrx$pivot[-(1:rank)] ## index of un-identifiable coefs
undrop <- 1:ncol(X)
if (length(drop)>0) undrop <- undrop[-drop]
......@@ -781,7 +781,7 @@ Sl.postproc <- function(Sl,fit,undrop,X0,cov=FALSE,scale = -1) {
## 3. At this stage reweight the model matrix in G if needed
## (e.g. in IRLS) to get X
## 4. um <- Sl.Xprep(Sl,X) to deal with identifiability and re-para.
## 5. initial smoothing parameters from mgcv:::initial.sp(X,G$S,G$off),
## 5. initial smoothing parameters from initial.sp(X,G$S,G$off),
## initial phi from, say variance of y over 10??
## 6. fit <- fast.REML.fit(um$Sl,um$X,G$y,rho,L=G$L,rho.0=G$lsp0,
## log.phi=log.phi,phi.fixed=FALSE/TRUE,Mp=um$Mp)
......
## R routines for gam fitting with calculation of derivatives w.r.t. sp.s
## (c) Simon Wood 2004-2009
## (c) Simon Wood 2004-2013
## These routines are for type 3 gam fitting. The basic idea is that a P-IRLS
## is run to convergence, and only then is a scheme for evaluating the
......@@ -81,32 +81,40 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
mustart = NULL, offset = rep(0, nobs),U1=diag(ncol(x)), Mp=-1, family = gaussian(),
control = gam.control(), intercept = TRUE,deriv=2,
gamma=1,scale=1,printWarn=TRUE,scoreType="REML",null.coef=rep(0,ncol(x)),
pearson.extra=0,dev.extra=0,n.true=-1,...)
## Version with new reparameterization and truncation strategy.
## ISSUES:
## 1. More efficient use of re-parameterization strategy could be employed,
## as repara is O(nq^2). Could pass in a repara object, containing sp's
## at last repara, X, rS, E, Eb in transformed form + T. If object is NULL
## that would signal need to create it. It would only be re-created if the log sps
## have moved far enough that some component is more than 3? from sps at last re-para.
## repara object would be returned at end of routine.
## --- note that this would require det|S| calculation to be performed separately
## in cases where re-parameterization is not done.
## --- not clear if this is really worth the effort: point is to avoid forming X%*%T
## every time routine is called, but this is only saving one O(np^2) operation
## from several to many.
## --- also this idea is difficult to make work with optimizers other than mine.
## This version is designed to allow iterative weights to be negative. This means that
## it deals with weights, rather than sqrt weights.
## deriv, sp, S, rS, H added to arg list.
## need to modify family before call.
## pearson.extra is an extra component to add to the pearson statistic in the P-REML/P-ML
## case, only.
## dev.extra is an extra component to add to the deviance in the REML and ML cases only.
## Similarly, n.true is to be used in place of the length(y) in ML/REML calculations,
## and the scale.est only.
{ if (family$link==family$canonical) fisher <- TRUE else fisher=FALSE
pearson.extra=0,dev.extra=0,n.true=-1,...) {
## Inputs:
## * x model matrix
## * y response
## * sp log smoothing parameters
## * Eb square root of nicely balanced total penalty matrix used for rank detection
## * UrS list of penalty square roots in range space of overall penalty. UrS[[i]]%*%t(UrS[[i]])
## is penalty. See 'estimate.gam' for more.
## * weights prior weights (reciprocal variance scale)
## * start initial values for parameters
## * etastart initial values for eta
## * mustart initial values for mu... only one of last## * control - control list.
## * intercept - indicates whether model has one.
## * deriv - order 0,1 or 2 derivatives are to be returned (lower is cheaper!)
## * gamma - multiplier for effective degrees of freedom in GCV/UBRE.
## * scale parameter. Negative signals to estimate.
## * printWarn print or supress?
## * scoreType - type of smoothness selection to use.
## * null.coef - coefficients for a null model, in order to be able to check for immediate
## divergence.
## * pearson.extra is an extra component to add to the pearson statistic in the P-REML/P-ML
## case, only.
## * dev.extra is an extra component to add to the deviance in the REML and ML cases only.
## * n.true is to be used in place of the length(y) in ML/REML calculations,
## and the scale.est only.
##
## Version with new reparameterization and truncation strategy. Allows iterative weights
## to be negative. Basically the workhorse routine for Wood (2011) JRSSB.
## A much modified version of glm.fit. Purpose is to estimate regression coefficients
## and compute a smoothness selection score along with its derivatives.
##
if (control$trace) { t0 <- proc.time();tc <- 0}
if (family$link==family$canonical) fisher <- TRUE else fisher=FALSE
## ... if canonical Newton = Fisher, but Fisher cheaper!
if (scale>0) scale.known <- TRUE else scale.known <- FALSE
if (!scale.known&&scoreType%in%c("REML","ML")) { ## the final element of sp is actually log(scale)
......@@ -194,7 +202,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
valideta <- function(eta) TRUE
validmu <- family$validmu
if (is.null(validmu))
validmu <- function(mu) TRUE
validmu <- function(mu) TRUE
if (is.null(mustart)) {
eval(family$initialize)
}
......@@ -302,19 +310,22 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
## Here a Fortran call has been replaced by pls_fit1 call
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),
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))
penalty=as.double(1),rank.tol=as.double(rank.tol),nt=as.integer(control$nthreads))
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),
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))
penalty=as.double(1),rank.tol=as.double(rank.tol),nt=as.integer(control$nthreads))
if (control$trace) tc <- tc + sum((proc.time()-t1)[c(1,4)])
}
start <- oo$y[1:ncol(x)];
......@@ -517,7 +528,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
if (scoreType%in%c("ML","P-ML")) {REML <- -1;remlInd <- 0}
if (REML==0) rSncol <- unlist(lapply(rS,ncol)) else rSncol <- unlist(lapply(UrS,ncol))
if (control$trace) t1 <- proc.time()
oo <- .C(C_gdi1,X=as.double(x[good,]),E=as.double(Sr),Eb = as.double(Eb),
rS = as.double(unlist(rS)),U1=as.double(U1),sp=as.double(exp(sp)),
z=as.double(z),w=as.double(w),wf=as.double(wf),alpha=as.double(alpha),
......@@ -532,9 +543,11 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
p=as.integer(ncol(x)),M=as.integer(nSp),Mp=as.integer(Mp),Enrow = as.integer(rows.E),
rSncol=as.integer(rSncol),deriv=as.integer(deriv.sp),
REML = as.integer(REML),fisher=as.integer(fisher),
fixed.penalty = as.integer(rp$fixed.penalty))
if (control$trace) cat("done!\n")
fixed.penalty = as.integer(rp$fixed.penalty),nthreads=as.integer(control$nthreads))
if (control$trace) {
tg <- sum((proc.time()-t1)[c(1,4)])
cat("done!\n")
}
rV <- matrix(oo$rV,ncol(x),ncol(x)) ## rV%*%t(rV)*scale gives covariance matrix
......@@ -691,7 +704,11 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
nulldf <- n.ok - as.integer(intercept)
aic.model <- aic(y, n, mu, weights, dev) # note: incomplete 2*edf needs to be added
if (control$trace) {
t1 <- proc.time()
at <- sum((t1-t0)[c(1,4)])
cat("Proportion time in C: ",(tc+tg)/at," ls:",tc/at," gdi:",tg/at,"\n")
}
list(coefficients = coef, residuals = residuals, fitted.values = mu,
family = family, linear.predictors = eta, deviance = dev,
......@@ -1232,7 +1249,7 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
else if (i==200) ct <- "iteration limit reached"
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
bfgs0 <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
control,gamma,scale,conv.tol=1e-6,maxNstep=5,maxSstep=2,
......@@ -2039,6 +2056,19 @@ fix.family.ls<-function(fam)
stop("family not recognised")
}
fix.family <- function(fam) {
## allows families to be patched...
if (fam$family[1]=="gaussian") { ## sensible starting values given link...
fam$initialize <- expression({
n <- rep.int(1, nobs)
if (family$link == "inverse") mustart <- y + (y==0)*sd(y)*.01 else
if (family$link == "log") mustart <- pmax(y,.01*sd(y)) else
mustart <- y
})
}
fam
}
negbin <- function (theta = stop("'theta' must be specified"), link = "log") {
## modified from Venables and Ripley's MASS library to work with gam.fit3,
......@@ -2117,7 +2147,7 @@ negbin <- function (theta = stop("'theta' must be specified"), link = "log") {
linkinv = stats$linkinv, variance = variance,dvar=dvar,d2var=d2var,d3var=d3var, dev.resids = dev.resids,
aic = aic, mu.eta = stats$mu.eta, initialize = initialize,ls=ls,
validmu = validmu, valideta = stats$valideta,getTheta = getTheta,qf=qf,rd=rd,canonical="log"), class = "family")
}
} ## negbin
......@@ -2268,7 +2298,7 @@ Tweedie <- function(p=1,link=power(0)) {
if (!is.null(stats$name))
linktemp <- stats$name
} else {
stop(gettextf("link \"%s\" not available for poisson family.",
stop(gettextf("link \"%s\" not available for Tweedie family.",
linktemp, collapse = ""),domain = NA)
}
}
......
......@@ -518,7 +518,7 @@ smooth2random.fs.interaction <- function(object,vnames,type=1) {
## penalties are not interwoven, but blocked (i.e. this ordering is
## as for gamm case).
if (object$fixed) return(list(fixed=TRUE,Xf=object$X))
if (type == 2) require(Matrix)
##if (type == 2) require(Matrix)
colx <- ncol(object$X)
diagU <- rep(1,colx)
ind <- 1:colx
......@@ -783,7 +783,7 @@ smooth2random.tensor.smooth <- function(object,vnames,type=1) {
gamm.setup<-function(formula,pterms,data=stop("No data supplied to gamm.setup"),knots=NULL,
gamm.setup <- function(formula,pterms,data=stop("No data supplied to gamm.setup"),knots=NULL,
parametric.only=FALSE,absorb.cons=FALSE)
## set up the model matrix, penalty matrices and auxilliary information about the smoothing bases
## needed for a gamm fit.
......@@ -835,6 +835,7 @@ gamm.setup<-function(formula,pterms,data=stop("No data supplied to gamm.setup"),
} else {
if (is.null(f.name)) f.name <- G$smooth[[i]]$fterm
else if (f.name!=G$smooth[[i]]$fterm) stop("only one level of smooth nesting is supported by gamm")
if (!is.null(attr(G$smooth[[i]],"del.index"))) stop("side conditions not allowed for nested smooths")
}
if (k < G$m) pord[(k+1):G$m] <- (1:G$m)[!done]
## .... ordered so that nested smooths are last
......@@ -989,7 +990,7 @@ extract.lme.cov2<-function(b,data,start.level=1)
## The above has been tested at nlme_3.1-109
## grps <- nlme::getGroups(b$modelStruct$corStruct) # replace grps (but not n.levels)
getGroupsFormula(b$modelStruct$corStruct)
vnames <- all.vars(nlme:::getGroupsFormula(b$modelStruct$corStruct))
vnames <- all.vars(nlme::getGroupsFormula(b$modelStruct$corStruct))
## attr(b$modelStruct$corStruct,"formula") # would include covariates
lab <- paste(eval(parse(text=vnames[1]),envir=b$data))
if (length(vnames)>1) for (i in 2:length(vnames)) {
......@@ -1362,7 +1363,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
# between the basis constructors provided in mgcv and the gammPQL routine used to estimate the model.
# NOTE: need to fill out the gam object properly
{
if (!require("nlme")) stop("gamm() requires package nlme to be installed")
##if (!require("nlme")) stop("gamm() requires package nlme to be installed")
control <- do.call("lmeControl",control)
# check that random is a named list
if (!is.null(random))
......@@ -1526,6 +1527,9 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
df.null=nrow(G$X),y=G$y,terms=gam.terms,pterms=pTerms,xlevels=G$xlevels,
contrasts=G$contrasts,assign=G$assign,na.action=attr(mf,"na.action"),
cmX=G$cmX,var.summary=G$var.summary,scale.estimated=TRUE)
pvars <- all.vars(delete.response(object$terms))
object$pred.formula <- if (length(pvars)>0) reformulate(pvars) else NULL
#######################################################
## Transform parameters back to the original space....
......
This diff is collapsed.
## (c) Simon N. Wood 2011-2013
## Many of the following are simple wrappers for C functions, used largely
## for testing purposes
block.reorder <- function(x,n.blocks=1,reverse=FALSE) {
## takes a matrix x divides it into n.blocks row-wise blocks, and re-orders
## so that the blocks are stored one after the other.
## e.g. library(mgcv); x <- matrix(1:18,6,3);xb <- mgcv:::block.reorder(x,2)
## x;xb;mgcv:::block.reorder(xb,2,TRUE)
r = nrow(x);cols = ncol(x);
if (n.blocks <= 1) return(x);
if (r%%n.blocks) {
nb = ceiling(r/n.blocks)
} else nb = r/n.blocks;
oo <- .C(C_row_block_reorder,x=as.double(x),as.integer(r),as.integer(cols),
as.integer(nb),as.integer(reverse));
matrix(oo$x,r,cols)
}
pqr <- function(x,nt=1) {
## parallel QR decomposition, using openMP in C, and up to nt threads (only if worthwhile)
## library(mgcv);n <- 20;p<-4;X <- matrix(runif(n*p),n,p);er <- mgcv:::pqr(X,nt=2)
x.c <- ncol(x);r <- nrow(x)
oo <- .C(C_mgcv_pqr,x=as.double(c(x,rep(0,nt*x.c^2))),as.integer(r),as.integer(x.c),
pivot=as.integer(rep(0,x.c)), tau=as.double(rep(0,(nt+1)*x.c)),as.integer(nt))
list(x=oo$x,r=r,c=x.c,tau=oo$tau,pivot=oo$pivot,nt=nt)
}
pqr.R <- function(x) {
## x is an object returned by pqr. This extracts the R factor...
## e.g. as pqr then...
## R <- mgcv:::pqr.R(er); R0 <- qr.R(qr(X,tol=0))
## svd(R)$d;svd(R0)$d
oo <- .C(C_getRpqr,R=as.double(rep(0,x$c^2)),as.double(x$x),as.integer(x$r),as.integer(x$c),
as.integer(x$c),as.integer(x$nt))
matrix(oo$R,x$c,x$c)
}
pqr.qy <- function(x,a,tr=FALSE) {
## x contains a parallel QR decomp as computed by pqr. a is a matrix. computes
## Qa or Q'a depending on tr.
## e.g. as above, then...
## a <- diag(p);Q <- mgcv:::pqr.qy(er,a);crossprod(Q)
## X[,er$pivot+1];Q%*%R
## Qt <- mgcv:::pqr.qy(er,diag(n),TRUE);Qt%*%t(Qt);range(Q-t(Qt))
## Q <- qr.Q(qr(X,tol=0));z <- runif(n);y0<-t(Q)%*%z
## mgcv:::pqr.qy(er,z,TRUE)->y
## z <- runif(p);y0<-Q%*%z;mgcv:::pqr.qy(er,z)->y
if (is.matrix(a)) a.c <- ncol(a) else a.c <- 1
if (tr) {
if (is.matrix(a)) { if (nrow(a) != x$r) stop("a has wrong number of rows") }
else if (length(a) != x$r) stop("a has wrong number of rows")
} else {
if (is.matrix(a)) { if (nrow(a) != x$c) stop("a has wrong number of rows") }
else if (length(a) != x$c) stop("a has wrong number of rows")
a <- c(a,rep(0,a.c*(x$r-x$c)))
}
oo <- .C(C_mgcv_pqrqy,a=as.double(a),as.double(x$x),as.double(x$tau),as.integer(x$r),
as.integer(x$c),as.integer(a.c),as.integer(tr),as.integer(x$nt))
if (tr) return(matrix(oo$a[1:(a.c*x$c)],x$c,a.c)) else
return(matrix(oo$a,x$r,a.c))
}
pmmult <- function(A,B,tA=FALSE,tB=FALSE,nt=1) {
## parallel matrix multiplication (not for use on vectors or thin matrices)
## library(mgcv);r <- 10;c <- 5;n <- 8
## A <- matrix(runif(r*n),r,n);B <- matrix(runif(n*c),n,c);range(A%*%B-mgcv:::pmmult(A,B,nt=1))
## A <- matrix(runif(r*n),n,r);B <- matrix(runif(n*c),n,c);range(t(A)%*%B-mgcv:::pmmult(A,B,TRUE,FALSE,nt=1))
## A <- matrix(runif(r*n),n,r);B <- matrix(runif(n*c),c,n);range(t(A)%*%t(B)-mgcv:::pmmult(A,B,TRUE,TRUE,nt=1))
## A <- matrix(runif(r*n),r,n);B <- matrix(runif(n*c),c,n);range(A%*%t(B)-mgcv:::pmmult(A,B,FALSE,TRUE,nt=1))
if (tA) { n = nrow(A);r = ncol(A)} else {n = ncol(A);r = nrow(A)}
if (tB) { c = nrow(B)} else {c = ncol(B)}
C <- rep(0,r * c)
oo <- .C(C_mgcv_pmmult,C=as.double(C),as.double(A),as.double(B),as.integer(tA),as.integer(tB),as.integer(r),
as.integer(c),as.integer(n),as.integer(nt));
matrix(oo$C,r,c)
}
\ No newline at end of file
......@@ -191,7 +191,7 @@ k.check <- function(b,subsample=5000,n.rep=400) {
} else modf <- b$model
nr <- length(rsd)
for (k in 1:m) { ## work through smooths
dat <- as.data.frame(mgcv:::ExtractData(b$smooth[[k]],modf,NULL)$data)
dat <- as.data.frame(ExtractData(b$smooth[[k]],modf,NULL)$data)
snames[k] <- b$smooth[[k]]$label
ind <- b$smooth[[k]]$first.para:b$smooth[[k]]$last.para
kc[k] <- length(ind)
......@@ -230,7 +230,7 @@ k.check <- function(b,subsample=5000,n.rep=400) {
}
}
nn <- 3
ni <- mgcv:::nearest(nn,as.matrix(dat))$ni
ni <- nearest(nn,as.matrix(dat))$ni
e <- rsd - rsd[ni[,1]]
for (j in 2:nn) e <- c(e,rsd-rsd[ni[,j]])
v.obs[k] <- mean(e^2)/2
......@@ -660,7 +660,8 @@ plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=
fac <- rep(x$flev,rep(n,nf))
dat <- data.frame(fac,xx)
names(dat) <- c(x$fterm,x$base$term)
X <- Predict.matrix.fs.interaction(x,dat)
# X <- Predict.matrix.fs.interaction(x,dat)
X <- PredictMat(x,dat)
if (is.null(xlab)) xlabel <- x$base$term else xlabel <- xlab
if (is.null(ylab)) ylabel <- label else ylabel <- ylab
return(list(X=X,scale=TRUE,se=FALSE,raw=raw,xlab=xlabel,ylab=ylabel,
......
......@@ -31,7 +31,7 @@ nat.param <- function(X,S,rank=NULL,type=0,tol=.Machine$double.eps^.8,unit.fnorm
## test code:
## x <- runif(100)
## sm <- smoothCon(s(x,bs="cr"),data=data.frame(x=x),knots=NULL,absorb.cons=FALSE)[[1]]
## np <- mgcv:::nat.param(sm$X,sm$S[[1]],type=3)
## np <- nat.param(sm$X,sm$S[[1]],type=3)
## range(np$X-sm$X%*%np$P)
if (type==2||type==3) { ## no need for QR step
er <- eigen(S,symmetric=TRUE)
......@@ -2133,7 +2133,7 @@ smooth.construct.mrf.smooth.spec <- function(object, data, knots) {
## natural parameterization given in Wood (2006) 4.1.14
if (object$bs.dim<length(levels(k))) { ## use low rank approx
rp <- mgcv:::nat.param(object$X,object$S[[1]],type=0)
rp <- nat.param(object$X,object$S[[1]],type=0)
np <- ncol(object$X)
## now retain only bs.dim least penalized elements
## of basis, which are the final bs.dim cols of rp$X
......
......@@ -171,8 +171,7 @@ crunch.knots <- function(G,knots,x0,y0,dx,dy)
setup.soap <- function(bnd,knots,nmax=100,k=10,bndSpec=NULL) {
## setup soap film smooth - nmax is number of grid cells for longest side
## it's important that grid cells are square!
require(Matrix)
## check boundary...
if (!inherits(bnd,"list")) stop("bnd must be a list.")
......@@ -272,7 +271,7 @@ soap.basis <- function(sd,x=NA,y=NA,film=TRUE,wiggly=TRUE,penalty=TRUE,plot=FALS
## If plot==TRUE then then data suitable for plotting are returned at the resolution
## of the solution grid. Then beta contains either the coefficients, or a single number
## representing the single basis function to return (0 for the offset).
require(Matrix)
if (!plot) {
indout <- inSide(sd$bnd,x,y); n <- length(x)
} else {
......
## (c) Simon N. Wood 2011
## (c) Simon N. Wood 2011-2013
## functions for sparse smoothing.
tri2nei <- function(T) {
## Each row of matrix T gives the indices of the points making up
## one triangle in d dimensions. T has d+1 columns. This routine
......@@ -19,7 +20,6 @@ tri.pen <- function(X,T) {
## finds a sparse approximate TPS penalty, based on the points in X,
## with triangulation T. Rows of X are points. Rows of T index vertices
## of triangles in X.
require(Matrix)
nn <- tri2nei(T) ## get neighbour list from T
## now obtain generalized FD penalty...
n <- nrow(X);d <- ncol(X);
......@@ -95,31 +95,34 @@ apply.spline <- function(spl,y) {
## Routines for sparse thin plate splines...
## kd tree/k nearest neighbout related routines....
kd.vis <- function(X,cex=.5) {
## code obtains and visualizes a kd tree for points in rows of X
kd.vis <- function(kd,X,cex=.5) {
## code visualizes a kd tree for points in rows of X
## kd <- kd.tree(X) produces correct tree.
if (ncol(X)!=2) stop("only deals with 2D case")
n <- nrow(X)
d <- ncol(X)
ind <- rind <- rep(0,n)
m <- 2;
while (m<n) m <- m*2
nb <- min(m-1,2*n-m/2-1)
lo <- hi <- rep(0,nb*d)
oo <- .C(C_Rkdtree,as.double(X),as.integer(n),as.integer(d),lo = as.double(lo),hi = as.double(hi),
ind = as.integer(ind), rind = as.integer(rind));
lo <- matrix(oo$lo,nb,d)
hi <- matrix(oo$hi,nb,d)
nb <- kd$idat[1]
dd <- matrix(kd$ddat[-1],nb,2*d,byrow=TRUE)
lo <- dd[,1:d];hi <- dd[,1:d+d]
rm(dd)
ll <- min(X[,1]); ul<- max(X[,1])
w <- ul-ll
ind <- lo[,1] < ll-w/10;lo[ind,1] <- ll-w/10
ind <- hi[,1] > ul+w/10;hi[ind,1] <- ul+w/10
ll <- min(X[,2]);ul <- max(X[,2])
w <- ul-ll
ind <- lo[,2] < ll-w/10;lo[ind,2] <- ll-w/10
ind <- hi[,2] > ul+w/10;hi[ind,2] <- ul+w/10
plot(X[,1],X[,2],pch=19,cex=cex,col=2)
for (i in 1:nb) {
rect(lo[i,1],lo[i,2],hi[i,1],hi[i,2])
}
#points(X[,1],X[,2],pch=19,cex=cex,col=2)
}
nearest <- function(k,X,gt.zero = FALSE,get.a=FALSE,balanced=FALSE,cut.off=5) {
nearest <- function(k,X,gt.zero = FALSE,get.a=FALSE) {
## The rows of X contain coordinates of points.
## For each point, this routine finds its k nearest
## neighbours, returning a list of 2, n by k matrices:
......@@ -129,14 +132,9 @@ nearest <- function(k,X,gt.zero = FALSE,get.a=FALSE,balanced=FALSE,cut.off=5) {
## neighbours.
## a - area associated with each point, if get.a is TRUE
## ties are broken arbitrarily.
## if balanced = TRUE, then some nearest neighbours are sacrificed
## for neighbours chosen to be on either side of the box in each
## direction in this case k>2*ncol(X). These neighbours are only used
## if closer than cut.off*max(k nearest distances).
## gt.zero indicates that neighbours must have distances greater
## than zero...
require(mgcv)
if