...
 
Commits (2)
** denotes quite substantial/important changes
*** denotes really big changes
Liable to change in future release:
- method="GCV.Cp" as default in gam. Default will
change to "REML". (1.8-24, June 2018)
Currently deprecated and liable to be removed:
- gam performance iteration (1.8-19, Sept 2017)
- argument 'pers' in plot.gam (1.8-23, Nov 2018)
- argument 'pers' in plot.gam (1.8-23, Nov 2017)
Issues:
* t2 in bam(...,discrete=TRUE) - not treated as tensor products at
present, and reparameterization needs checking (also for bam).
1.8-24
* Extended Fellner Schall optimizer now avaialable for all families with 'gam'
using gam(...,optimizer="efs").
* Change to default behaviour of plot.gam when 'seWithMean=TRUE', and of
predict.gam when 'type="iterms"'. The extra uncertainty added to CIs or
standard errors now reflects the uncertainty in the mean in all other model
terms, not just the uncertanity in the mean of the fixed effects as before.
See ?plot.gam and ?predict.gam (including for how to get the old behaviour).
* 're' smooths can now accept matrix arguments: see ?linear.functional.terms.
* cox.ph now allows an offset to be provided.
* Fix in smoothCon for bug in case in which only a single coefficient is
involved in a sum-to-zero constraint. Could cause failure in e.g. t2 with
cc marginal.
* Model terms s, te etc are now always evaluated in mgcv workspace explicitly
to avoid masking problems in obscure circumstances.
* 'mrf' smooth documentation modified to make it clearer how to specify 'nb',
and code modified so that it is now possible to speficy the neighbour
structure using names rather than indices.
* 'bfgs' fix to handle extended families.
* plot.gam modified to only prompt (via devAskNewPage) for a new page after
the first page is used up.
* export 'k.check'.
* Fix to 'Rrank'. Previously a matrix R with more columns than rows could
cause a segfault.
* Fix to non-finite likelihood handling in gam.fit5.
* Fix in bgam.fitd to step reduce under indefinite deviance and to ensure
penalty evaluation is round off negative proof.
* newton slighty modified to avoid (small) chance of all sp's being dropped
for indef likelihood.
1.8-23
......
Package: mgcv
Version: 1.8-23
Version: 1.8-24
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with Automatic Smoothness
......@@ -18,6 +18,6 @@ LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
NeedsCompilation: yes
Packaged: 2018-01-08 15:15:43 UTC; sw283
Packaged: 2018-06-18 13:08:20 UTC; sw283
Repository: CRAN
Date/Publication: 2018-01-15 00:23:36 UTC
Date/Publication: 2018-06-18 19:07:59 UTC
This diff is collapsed.
......@@ -14,7 +14,7 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
gaulss,gam.side,get.var,gevlss,
influence.gam,
in.out,inSide,interpret.gam,initial.sp,
jagam,ldetS,
jagam,k.check,ldetS,
ldTweedie,
logLik.gam,ls.size,
magic, magic.post.proc, model.matrix.gam,mini.roots,
......
......@@ -211,7 +211,7 @@ discrete.mf <- function(gp,mf,names.pmf,m=NULL,full=TRUE) {
## * nr records the number of unique discretized covariate values
## i.e. the number of rows before the padding starts
## * k.start contains the starting column in index vector k, for
## each variable.
## each variable. The final element is the column beyond the last one.
## * k is the index matrix. The ith record of the 1st column of the
## jth variable is in row k[i,k.start[j]] of the corresponding
## column of mf.
......@@ -527,7 +527,7 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
for (iter in 1L:control$maxit) { ## main fitting loop
devold <- dev
dev <- 0
## accumulate the QR decomposition of the weighted model matrix
if (iter==1||!additive) {
qrx <- list()
......@@ -536,10 +536,13 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
eta <- Xbd(G$Xd,coef,G$kd,G$ks,G$ts,G$dt,G$v,G$qc,G$drop) + offset
lsp.full <- G$lsp0
if (n.sp>0) lsp.full <- lsp.full + if (is.null(G$L)) lsp[1:n.sp] else G$L %*% lsp[1:n.sp]
Sb <- Sl.Sb(Sl,lsp.full,prop$beta) ## store S beta to allow rapid step halving
#Sb <- Sl.Sb(Sl,lsp.full,prop$beta) ## store S beta to allow rapid step halving
rSb <- Sl.rSb(Sl,lsp.full,prop$beta) ## store S beta to allow rapid step halving
if (iter>2) {
Sb0 <- Sl.Sb(Sl,lsp.full,b0)
bSb0 <- sum(b0*Sb0) ## penalty at start of beta step
#Sb0 <- Sl.Sb(Sl,lsp.full,b0)
#bSb0 <- sum(b0*Sb0) ## penalty at start of beta step
rSb0 <- Sl.rSb(Sl,lsp.full,b0)
bSb0 <- sum(rSb0^2)
## get deviance at step start, with current theta if efam
dev0 <- if (efam) sum(family$dev.resids(G$y,mu0,G$w,theta)) else
sum(family$dev.resids(G$y,mu0,G$w))
......@@ -551,10 +554,12 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
dev <- if (efam) sum(family$dev.resids(G$y,mu,G$w,theta)) else
sum(family$dev.resids(G$y,mu,G$w))
if (iter>2) { ## coef step length control
bSb <- sum(prop$beta*Sb) ## penalty at end of beta step
if (dev0 + bSb0 < dev + bSb && kk < 30) { ## beta step not improving current pen dev
#bSb <- sum(prop$beta*Sb) ## penalty at end of beta step
bSb <- sum(rSb^2) ## penalty at end of beta step
if ((!is.finite(dev) || dev0 + bSb0 < dev + bSb) && kk < 30) { ## beta step not improving current pen dev
coef <- (coef0 + coef)/2 ## halve the step
Sb <- (Sb0 + Sb)/2
#Sb <- (Sb0 + Sb)/2
rSb <- (rSb0 + rSb)/2
eta <- (eta0 + eta)/2
prop$beta <- (b0 + prop$beta)/2
kk <- kk + 1
......@@ -567,7 +572,8 @@ bgam.fitd <- function (G, mf, gp ,scale , coef=NULL,etastart = NULL,
eta0 <- eta
mu0 <- mu
b0 <- prop$beta ## beta repara
dev <- dev + sum(prop$beta*Sb) ## add penalty to deviance
#dev <- dev + sum(prop$beta*Sb) ## add penalty to deviance
dev <- dev + sum(rSb^2)
} else reml <- dev ## for convergence checking
if (efam) { ## extended family
......@@ -2113,7 +2119,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
} else { ## not a tensor smooth
v[[kb]] <- rep(0,0)
dt[kb] <- dt[kb] + 1
G$Xd[[k]] <- G$X[1:dk$nr[k],G$smooth[[i]]$first.para:G$smooth[[i]]$last.para]
G$Xd[[k]] <- G$X[1:dk$nr[k],G$smooth[[i]]$first.para:G$smooth[[i]]$last.para,drop=FALSE]
k <- k + 1
}
kb <- kb + 1
......
......@@ -39,6 +39,7 @@ cox.ph <- function (link = "identity") {
G$X <- G$X[y.order,,drop=FALSE]
attributes(G$X) <- attrX
G$w <- G$w[y.order]
G$offset <- G$offset[y.order]
})
postproc <- expression({
......@@ -46,7 +47,7 @@ cox.ph <- function (link = "identity") {
## baseline hazard estimation...
## first get the estimated hazard and prediction information...
G$X <- Sl.initial.repara(G$Sl,G$X,inverse=TRUE,cov=FALSE,both.sides=FALSE)
object$family$data <- G$family$hazard(G$y,G$X,object$coefficients,G$w)
object$family$data <- G$family$hazard(G$y,G$X,object$coefficients,G$w,G$offset)
rumblefish <- G$family$hazard(G$y,matrix(0,nrow(G$X),0),object$coefficients,G$w)
s0.base <- exp(-rumblefish$h[rumblefish$r]) ## no model baseline survival
s0.base[s0.base >= 1] <- 1 - 2*.Machine$double.eps ## avoid NA later
......@@ -69,7 +70,7 @@ cox.ph <- function (link = "identity") {
if (is.null(start)) start <- rep(0,ncol(x))
})
hazard <- function(y, X,beta,wt) {
hazard <- function(y, X,beta,wt,offset=0) {
## get the baseline hazard function information, given times in descending order in y
## model matrix (same ordering) in X, coefs in beta and censoring in wt (1 = death, 0
## = censoring)
......@@ -96,7 +97,7 @@ cox.ph <- function (link = "identity") {
}
#tr <- unique(y);nt <- length(tr)
p <- ncol(X)
eta <- as.double(X%*%beta)
eta <- as.double(X%*%beta) + offset
if (ns==1) {
r <- match(y,tr)
oo <- .C("coxpp",eta,A=as.double(X),as.integer(r),d=as.integer(wt),
......@@ -156,9 +157,10 @@ cox.ph <- function (link = "identity") {
}
if (sum(is.na(y))>0) stop("NA times supplied for cox.ph prediction")
X <- X[ii,,drop=FALSE];y <- y[ii];
n <- nrow(X)
if (is.null(off)) off <- rep(0,n)
if (is.null(strat)) {
n <- nrow(X)
oo <- .C("coxpred",as.double(X),t=as.double(y),as.double(beta),as.double(Vb),
oo <- .C("coxpred",as.double(X),t=as.double(y),as.double(beta),as.double(off),as.double(Vb),
a=as.double(family$data$a),h=as.double(family$data$h),q=as.double(family$data$q),
tr = as.double(family$data$tr),
n=as.integer(n),p=as.integer(ncol(X)),nt = as.integer(family$data$nt),
......@@ -175,7 +177,7 @@ cox.ph <- function (link = "identity") {
ind <- which(strat==pstrata[i]) ## prediction data index
trind <- which(family$data$tr.strat == pstrata[i])
n <- length(ind)
oo <- .C("coxpred",as.double(X[ind,]),t=as.double(y[ind]),as.double(beta),as.double(Vb),
oo <- .C("coxpred",as.double(X[ind,]),t=as.double(y[ind]),as.double(beta),as.double(off),as.double(Vb),
a=as.double(family$data$a[,trind]),h=as.double(family$data$h[trind]),q=as.double(family$data$q[trind]),
tr = as.double(family$data$tr[trind]),
n=as.integer(n),p=as.integer(p),nt = as.integer(length(trind)),
......@@ -205,7 +207,8 @@ cox.ph <- function (link = "identity") {
## or its Choleski factor
## D is the diagonal pre-conditioning matrix used to obtain Hp
## 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")
#!! if (!is.null(offset)&&sum(offset!=0)) stop("cox.ph does not yet handle offsets")
if (is.matrix(y)) {
## first column is time, second is *numeric* code indicating strata
......@@ -217,7 +220,7 @@ cox.ph <- function (link = "identity") {
p <- ncol(X)
deriv <- deriv - 1
mu <- X%*%coef
mu <- X%*%coef + offset
g <- rep(0,p);H <- rep(0,p*p)
if (deriv > 0) {
M <- ncol(d1b)
......
## (c) Simon N. Wood (ocat, tw, nb, ziP) & Natalya Pya (scat, beta),
## (c) Simon N. Wood & Natalya Pya (scat, beta),
## 2013-2017. Released under GPL2.
## See gam.fit4.r for testing functions fmud.test and fetad.test.
estimate.theta <- function(theta,family,y,mu,scale=1,wt=1,tol=1e-7,attachH=FALSE) {
## Simple Newton iteration to estimate theta for an extended family,
## given y and mu. To be iterated with estimation of mu given theta.
## If used within a PIRLS loop then divergence testing of coef update
## will have to re-compute pdev after theta update.
## Not clear best way to handle scale - could oprimize here as well
## Not clear best way to handle scale - could optimize here as well
if (!inherits(family,"extended.family")) stop("not an extended family")
......@@ -115,15 +116,15 @@ find.null.dev <- function(family,y,eta,offset,weights) {
## function fix.family.link.extended.family with functions
## gkg where k is 2,3 or 4 giving the kth derivative of the
## link over the first derivative of the link to the power k.
## for non standard links these functions muct be supplied.
## for non standard links these functions must be supplied.
## dev.resids - function computing deviance residuals.
## Dd - function returning derivatives of deviance residuals w.r.t. mu and theta.
## aic - function computing twice - log likelihood for 2df to be added to.
## aic - function computing twice -ve log likelihood for 2df to be added to.
## initialize - expression to be evaluated in gam.fit4 and initial.spg
## to initialize mu or eta.
## preinitialize - optional function of y and family, returning a list with optional elements
## Theta - intitial Theta and y - modified y for use in fitting (see e.g. ocat and betar)
## postproc - function with arguments family,y,prior.weights,fitted,linear.predictors,offset,intercept
## postproc - function with arguments family, y, prior.weights, fitted, linear.predictors, offset, intercept
## to call after fit to compute (optionally) the label for the family, deviance and null deviance.
## See ocat for simple example and betar or ziP for complicated. Called in estimate.gam.
## ls - function to evaluated log saturated likelihood and derivatives w.r.t.
......@@ -140,7 +141,7 @@ find.null.dev <- function(family,y,eta,offset,weights) {
## predict - optional function for predicting from model, called by predict.gam.
## family$data - optional list storing any family specific data for use, e.g. in predict
## function. - deprecated (commented out below - appears to be used nowhere)
## scale - < 0 to estimate. ignored if NULL
## extended family object for ordered categorical
......@@ -254,7 +255,7 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
} ## end of dev.resids
Dd <- function(y, mu, theta, wt=NULL, level=0) {
## derivatives of the deviance...
## derivatives of the ocat deviance...
# F <- function(x) { ## e^(x)/(1+e^x) without overflow
# h <- ind <- x > 0; h[ind] <- 1/(exp(-x[ind]) + 1)
# x <- exp(x[!ind]); h[!ind] <- (x/(1+x))
......@@ -387,7 +388,8 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
oo$Dth[ind,k] <- Da0[ind]*etk
oo$Dmuth[ind,k] <- Dmua0[ind]*etk
oo$Dmu2th[ind,k] <- Dmu2a0[ind]*etk
}
}
oo$EDmu2th <- oo$Dmu2th
}
if (level >1) { ## and the second derivative components
oo$Dmu4 <- 2*((3*b^2 + 4*a*c)/f + a2*(6*a2/f - 12*b)/f2 - d)/f
......@@ -440,7 +442,7 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
}
}
oo
} ## end of Dd
} ## end of Dd (ocat)
aic <- function(y, mu, theta=NULL, wt, dev) {
......@@ -473,77 +475,7 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
ls <- function(y,w,theta,scale) {
## the log saturated likelihood function.
#! actually only first line used since re-def as 0
#vec <- !is.null(attr(theta,"vec.grad"))
#lsth1 <- if (vec) matrix(0,length(y),R-2) else rep(0,R-2)
return(list(ls=0,lsth1=rep(0,R-2),lsth2=matrix(0,R-2,R-2)))
F <- function(x) {
h <- ind <- x > 0; h[ind] <- 1/(exp(-x[ind]) + 1)
x <- exp(x[!ind]); h[!ind] <- (x/(1+x))
h
}
Fdiff <- function(a,b) {
## cancellation resistent F(b)-F(a), b>a
h <- rep(1,length(b)); h[b>0] <- -1; eb <- exp(b*h)
h <- h*0+1; h[a>0] <- -1; ea <- exp(a*h)
ind <- b<0;bi <- eb[ind];ai <- ea[ind]
h[ind] <- bi/(1+bi) - ai/(1+ai)
ind1 <- a>0;bi <- eb[ind1];ai <- ea[ind1]
h[ind1] <- (ai-bi)/((ai+1)*(bi+1))
ind <- !ind & !ind1;bi <- eb[ind];ai <- ea[ind]
h[ind] <- (1-ai*bi)/((bi+1)*(ai+1))
h
}
R = length(theta)+2
alpha <- rep(0,R+1) ## the thresholds
alpha[1] <- -Inf;alpha[R+1] <- Inf
alpha[2] <- -1
if (R > 2) {
ind <- 3:R
alpha[ind] <- alpha[2] + cumsum(exp(theta))
}
al1 <- alpha[y+1];al0 = alpha[y]
g1 <- F((al1-al0)/2);g0 <- F((al0-al1)/2)
##A <- pmax(g1 - g0,.Machine$double.eps)
A <- Fdiff((al0-al1)/2,(al1-al0)/2)
ls <- sum(log(A))
B <- g1^2 - g1 + g0^2 - g0
C <- 2 * g1^3 - 3 * g1^2 + g1 - 2 * g0^3 + 3 * g0^2 - g0
Da0 <- .5 * B/A ; Da1 <- -0.5 *B/A
Da0a0 <- .25 * C/A - .25 * B^2/A^2
Da1a1 <- .25 * C/A - .25 * B^2/A^2
Da0a1 <- - .25 * C/A + .25 * B^2/A^2
i <- 0
n2d <- (R-2)*(R-1)/2
n <- length(y)
Dth <- matrix(0,n,R-2)
Dth2 <- matrix(0,n,n2d)
for (j in 1:(R-2)) for (k in j:(R-2)) {
i <- i + 1 ## the second deriv col
ind <- y >= j ## rest are zero
ar1.k <- ar.k <- rep(exp(theta[k]),n)
ar.k[y==R | y <= k] <- 0; ar1.k[y<k+2] <- 0
ar.j <- ar1.j <- rep(exp(theta[j]),n)
ar.j[y==R | y <= j] <- 0; ar1.j[y<j+2] <- 0
ar.kj <- ar1.kj <- rep(0,n)
if (k==j) {
ar.kj[y>k&y<R] <- exp(theta[k])
ar1.kj[y>k+1] <- exp(theta[k])
Dth[ind,k] <- Da1[ind]*ar.k[ind] + Da0[ind]*ar1.k[ind]
}
Dth2[,i] <- Da1a1*ar.k*ar.j + Da0a1*ar.k*ar1.j + Da1 * ar.kj +
Da0a0*ar1.k*ar1.j + Da0a1*ar1.k*ar.j + Da0 * ar1.kj
}
lsth2=colSums(Dth2)
if (R>2) {
ls2 <- matrix(0,R-2,R-2);ii <- 0
for (i in 1:(R-2)) for (j in i:(R-2)) {
ii <- ii + 1
ls2[i,j] <- ls2[j,i] <- lsth2[ii]
}
}
list(ls=ls,lsth1=colSums(Dth),lsth2=ls2)
} ## end of ls
## initialization is interesting -- needs to be with reference to initial cut-points
......@@ -748,7 +680,7 @@ nb <- function (theta = NULL, link = "log") {
}
Dd <- function(y, mu, theta, wt, level=0) {
## derivatives of the deviance...
## derivatives of the nb deviance...
##ltheta <- theta
theta <- exp(theta)
yth <- y + theta
......@@ -765,6 +697,7 @@ nb <- function (theta = NULL, link = "log") {
r$Dmuth <- 2 * wt * theta * (1 - yth/muth)/muth
r$Dmu3 <- 4 * wt * (yth/muth^3 - y/mu^3)
r$Dmu2th <- 2 * wt * theta * (2*yth/muth - 1)/muth^2
r$EDmu2th <- 2 * wt / muth^2
}
if (level>1) { ## whole damn lot
r$Dmu4 <- 2 * wt * (6*y/mu^4 - 6*yth/muth^4)
......@@ -918,7 +851,7 @@ tw <- function (theta = NULL, link = "log",a=1.01,b=1.99) {
}
Dd <- function(y, mu, theta, wt, level=0) {
## derivatives of the deviance...
## derivatives of the tw deviance...
a <- get(".a");b <- get(".b")
th <- theta
p <- if (th>0) (b+a*exp(-th))/(1+exp(-th)) else (b*exp(th)+a)/(exp(th)+1)
......@@ -1116,7 +1049,7 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
r$Dth <- 2 * wt *theta*(-mu*log.yoney - log1p(-y)+ mu*psi0.muth+onemu*psi0.onemuth -psi0.th)
r$Dmuth <- r$Dmu + 2 * wt * theta^2*(mu*psi1.muth -onemu*psi1.onemuth)
r$Dmu3 <- 2 * wt *theta^3 * (psi2.muth - psi2.onemuth)
r$Dmu2th <- 2* r$Dmu2 + 2 * wt * theta^3* (mu*psi2.muth + onemu*psi2.onemuth)
r$EDmu2th <- r$Dmu2th <- 2* r$Dmu2 + 2 * wt * theta^3* (mu*psi2.muth + onemu*psi2.onemuth)
}
if (level>1) { ## whole lot
r$Dmu4 <- 2 * wt *theta^4 * (psi3.muth+psi3.onemuth)
......@@ -1387,7 +1320,7 @@ scat <- function (theta = NULL, link = "identity",min.df = 3) {
}
Dd <- function(y, mu, theta, wt, level=0) {
## derivatives of the deviance...
## derivatives of the scat deviance...
## ltheta <- theta
min.df <- get(".min.df")
nu <- exp(theta[1])+min.df; sig <- exp(theta[2])
......@@ -1847,6 +1780,7 @@ ziP <- function (theta = NULL, link = "identity",b=0) {
mu <- object$linear.predictors
wts <- object$prior.weights
res <- object$family$dev.resids(y,mu,wts)
## next line is correct as function returns -2*saturated.log.lik
res <- res - object$family$saturated.ll(y,object$family,wts)
fv <- predict.gam(object,type="response")
s <- attr(res,"sign")
......@@ -1854,7 +1788,7 @@ ziP <- function (theta = NULL, link = "identity",b=0) {
res <- as.numeric(sqrt(pmax(res,0)) * s)
}
res
} ## residuals
} ## residuals (ziP)
predict <- function(family,se=FALSE,eta=NULL,y=NULL,X=NULL,
beta=NULL,off=NULL,Vb=NULL) {
......
......@@ -241,7 +241,7 @@ Sl.Sb <- function(Sl,rho,beta) {
if (length(Sl[[b]]$S)==1) { ## singleton - multiple of identity
a[ind] <- a[ind] + beta[ind] * exp(rho[k])
k <- k + 1
} else { ## mulit-S block
} else { ## multi-S block
for (j in 1:length(Sl[[b]]$S)) {
a[ind] <- a[ind] + exp(rho[k]) * (Sl[[b]]$S[[j]] %*% beta[ind])
k <- k + 1
......@@ -251,6 +251,37 @@ Sl.Sb <- function(Sl,rho,beta) {
a
} ## Sl.Sb
Sl.rSb <- function(Sl,rho,beta) {
## Computes vector 'a' containing all terms rS %*% beta stacked end to end.
## sum of squares of 'a' this is bSb, but 'a' is linear in beta
## Assumes initial re-parameterization has taken
## place, so single penalties are multiples of identity and uses S for
## multi-S blocks. Logic is identical to Sl.addS.
k <- 1 ## sp counter
kk <- 0 ## total length of returned vector
if (length(Sl)>0) for (b in 1:length(Sl)) {
ind <- (Sl[[b]]$start:Sl[[b]]$stop)[Sl[[b]]$ind]
kk <- kk + length(Sl[[b]]$S)*length(ind)
}
a <- rep(0,kk)
kk <- 0
if (length(Sl)>0) for (b in 1:length(Sl)) {
ind <- (Sl[[b]]$start:Sl[[b]]$stop)[Sl[[b]]$ind]
if (length(Sl[[b]]$S)==1) { ## singleton - multiple of identity
a[kk + 1:length(ind)] <- beta[ind] * exp(rho[k]/2)
k <- k + 1
kk <- kk + length(ind)
} else { ## multi-S block
for (j in 1:length(Sl[[b]]$S)) {
a[kk + 1:length(ind)] <- exp(rho[k]/2) * (beta[ind] %*% Sl[[b]]$rS[[j]])
k <- k + 1
kk <- kk + length(ind)
}
}
}
a
} ## Sl.rSb
Sl.initial.repara <- function(Sl,X,inverse=FALSE,both.sides=TRUE,cov=TRUE,nt=1) {
## Routine to apply initial Sl re-parameterization to model matrix X,
## or, if inverse==TRUE, to apply inverse re-para to parameter vector
......@@ -879,7 +910,7 @@ Sl.fitChol <- function(Sl,XX,f,rho,yy=0,L=NULL,rho0=0,log.phi=0,phi.fixed=TRUE,
uconv.ind <- (abs(reml1) > tol)|(abs(diag(reml2))>tol)
hess <- reml2
grad <- reml1
if (length(grad)>0) {
if (length(grad)>0&&sum(uconv.ind)>0) {
if (sum(uconv.ind)!=ncol(reml2)) {
reml1 <- reml1[uconv.ind]
reml2 <- reml2[uconv.ind,uconv.ind]
......
......@@ -118,7 +118,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
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)
if (!scale.known&&scoreType%in%c("REML","ML","EFS")) { ## the final element of sp is actually log(scale)
nsp <- length(sp)
scale <- exp(sp[nsp])
sp <- sp[-nsp]
......@@ -137,9 +137,10 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
q <- ncol(x)
if (length(UrS)) { ## find a stable reparameterization...
grderiv <- deriv*as.numeric(scoreType%in%c("REML","ML","P-REML","P-ML"))
grderiv <- if (scoreType=="EFS") 1 else deriv*as.numeric(scoreType%in%c("REML","ML","P-REML","P-ML"))
rp <- gam.reparam(UrS,sp,grderiv) ## note also detects fixed penalty if present
## Following is for debugging only...
# deriv.check <- FALSE
......@@ -177,7 +178,8 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
rows.E <- q-Mp
Sr <- cbind(rp$E,matrix(0,nrow(rp$E),Mp))
St <- rbind(cbind(rp$S,matrix(0,nrow(rp$S),Mp)),matrix(0,Mp,q))
} else {
} else {
grderiv <- 0
T <- diag(q);
St <- matrix(0,q,q)
rSncol <- sp <- rows.E <- Eb <- Sr <- 0
......@@ -185,7 +187,12 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
rp <- list(det=0,det1 = rep(0,0),det2 = rep(0,0),fixed.penalty=FALSE)
}
iter <- 0;coef <- rep(0,ncol(x))
if (scoreType=="EFS") {
scoreType <- "REML" ## basically optimizing REML
deriv <- 0 ## only derivatives of log|S|_+ required (see above)
}
conv <- FALSE
n <- nobs <- NROW(y) ## n is just to keep codetools happy
if (n.true <= 0) n.true <- nobs ## n.true is used in criteria in place of nobs
......@@ -773,7 +780,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
boundary = boundary,D1=D1,D2=D2,P=P,P1=P1,P2=P2,trA=trA,trA1=trA1,trA2=trA2,
GCV=GCV,GCV1=GCV1,GCV2=GCV2,GACV=GACV,GACV1=GACV1,GACV2=GACV2,UBRE=UBRE,
UBRE1=UBRE1,UBRE2=UBRE2,REML=REML,REML1=REML1,REML2=REML2,rV=rV,db.drho=db.drho,
dw.drho=dw.drho,dVkk = matrix(oo$dVkk,nSp,nSp),
dw.drho=dw.drho,dVkk = matrix(oo$dVkk,nSp,nSp),ldetS1 = if (grderiv) rp$det1 else 0,
scale.est=scale.est,reml.scale= reml.scale,aic=aic.model,rank=oo$rank.est,K=Kmat)
} ## end gam.fit3
......@@ -1338,8 +1345,10 @@ newton <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
## tiny relative to largest, as this space is likely to be poorly
## modelled on scale of Newton step...
uconv.ind <- uconv.ind & abs(grad)>max(abs(grad))*.001
uconv.ind1 <- uconv.ind & abs(grad)>max(abs(grad))*.001
if (sum(uconv.ind1)==0) uconv.ind1 <- uconv.ind ## nothing left reset
if (sum(uconv.ind)==0) uconv.ind[which(abs(grad)==max(abs(grad)))] <- TRUE ## need at least 1 to update
## exclude apparently converged gradients from computation
hess1 <- hess[uconv.ind,uconv.ind]
grad1 <- grad[uconv.ind]
......@@ -1776,8 +1785,22 @@ bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
score <- b$UBRE;grad <- t(L)%*%b$UBRE1
} else { ## default to deviance based GCV
score <- b$GCV;grad <- t(L)%*%b$GCV1;
}
## dVkk only refers to smoothing parameters, but sp may contain
## extra parameters at start and scale parameter at end. Have
## to reduce L accordingly...
if (!is.null(family$n.theta)&&family$n.theta>0) {
ind <- 1:family$n.theta
nind <- ncol(L) - family$n.theta - if (family$n.theta + nrow(b$dVkk)<nrow(L)) 1 else 0
spind <- if (nind>0) family$n.theta+1:nind else rep(0,0)
rspind <- family$n.theta + 1:nrow(b$dVkk)
} else {
nind <- ncol(L) - if (nrow(b$dVkk)<nrow(L)) 1 else 0
spind <- if (nind>0) 1:nind else rep(0,0) ## index of smooth parameters
rspind <- 1:nrow(b$dVkk)
}
L0 <- if (nrow(L)==nrow(b$dVkk)) L else L[-nrow(L),-ncol(L)]
L0 <- L[rspind,spind] ##if (nrow(L)!=nrow(b$dVkk)) L[spind,spind] else L
initial$dVkk <- diag(t(L0) %*% b$dVkk %*% L0)
initial$score <- score;initial$grad <- grad;
initial$scale.est <- b$scale.est
......@@ -2004,9 +2027,10 @@ bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
else score.scale <- abs(trial$scale.est) + abs(trial$score) ##trial$dev/nrow(X) + abs(trial$score)
uconv.ind <- abs(trial$grad) > score.scale*conv.tol
if (sum(uconv.ind)) converged <- FALSE
if (length(uconv.ind)>length(trial$dVkk)) trial$dVkk <- c(trial$dVkk,score.scale)
#if (length(uconv.ind)>length(trial$dVkk)) trial$dVkk <- c(trial$dVkk,score.scale)
## following must be tighter than convergence...
uconv.ind <- abs(trial$grad) > score.scale*conv.tol*.1 | abs(trial$dVkk) > score.scale * conv.tol*.1
uconv.ind <- abs(trial$grad) > score.scale*conv.tol*.1
uconv.ind[spind] <- uconv.ind[spind] | abs(trial$dVkk) > score.scale * conv.tol*.1
if (abs(initial$score-trial$score) > score.scale*conv.tol) {
if (!sum(uconv.ind)) uconv.ind <- uconv.ind | TRUE ## otherwise can't progress
converged <- FALSE
......@@ -2045,10 +2069,11 @@ bfgs <- function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
trial$dscore <- sum(trial$grad*step)
trial$scale.est <- b$scale.est
trial$dVkk <- diag(t(L0) %*% b$dVkk %*% L0) ## curvature testing matrix
if (length(uconv.ind)>length(trial$dVkk)) trial$dVkk <- c(trial$dVkk,score.scale)
#if (length(uconv.ind)>length(trial$dVkk)) trial$dVkk <- c(trial$dVkk,score.scale)
rm(b);counter <- counter + 1
## note that following rolls back until there is clear signal in derivs...
uconv.ind0 <- abs(trial$grad) > score.scale*conv.tol*20 | abs(trial$dVkk) > score.scale * conv.tol * 20
uconv.ind0 <- abs(trial$grad) > score.scale*conv.tol*20
uconv.ind0[spind] <- uconv.ind0[spind] | abs(trial$dVkk) > score.scale * conv.tol * 20
uconv.ind0 <- uconv.ind0 | uconv.ind ## make sure we don't start rolling back unproblematic sps
}
uconv.ind <- uconv.ind | TRUE
......
This diff is collapsed.
......@@ -987,7 +987,7 @@ extract.lme.cov2<-function(b,data,start.level=1)
}
}
list(V=V,ind=Cind)
}
} ## extract.lme.cov2
extract.lme.cov<-function(b,data,start.level=1)
# function to extract the response data covariance matrix from an lme fitted
......@@ -1071,7 +1071,7 @@ extract.lme.cov<-function(b,data,start.level=1)
V <- V+Z%*%Vr%*%t(Z)
}
V
}
} ## extract.lme.cov
formXtViX <- function(V,X)
## forms X'V^{-1}X as efficiently as possible given the structure of
......@@ -1147,7 +1147,7 @@ gammPQL <- function (fixed, random, family, data, correlation, weights,
zz <- eta + fit0$residuals - off
wz <- fit0$weights
fam <- family
## find non clashing name for pseudodata and insert in formula
zz.name <- new.name("zz",names(data))
eval(parse(text=paste("fixed[[2]] <- quote(",zz.name,")")))
......@@ -1188,6 +1188,13 @@ gammPQL <- function (fixed, random, family, data, correlation, weights,
if (!converged) warning("gamm not converged, try increasing niterPQL")
fit$y <- fit0$y
fit$w <- w ## prior weights
## would require full edf to be computable with re terms include!
#if (is.null(correlation)) { ## then a conditional AIC is possible
# y <- fit$y;weights <- w; nobs <- length(y)
# eval(fam$initialize)
# dev <- sum(fam$dev.resids(y,mu,w))
# fit$aic <- fam$aic(y,n,mu,w,dev)
#}
fit
}
......@@ -1350,15 +1357,18 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
length(offset.name)==0) lme.used <- TRUE else lme.used <- FALSE
if (lme.used&&!is.null(weights)&&!wisvf) lme.used <- FALSE
if (lme.used)
{ ## following construction is a work-around for problem in nlme 3-1.52
if (lme.used) { ## following construction is a work-around for problem in nlme 3-1.52
eval(parse(text=paste("ret$lme<-lme(",deparse(fixed.formula),
",random=rand,data=strip.offset(mf),correlation=correlation,",
"control=control,weights=weights,method=method)"
,sep="" )))
,sep="" )))
## need to be able to work out full edf for following to work...
# if (is.null(correlation)) { ## compute conditional aic precursor
# dev <- sum(family$dev.resids(G$y,fitted(ret$lme),weights))
# ret$lme$aic <- family$aic(G$y,1,fitted(ret$lme),weights,dev)
# }
##ret$lme<-lme(fixed.formula,random=rand,data=mf,correlation=correlation,control=control)
} else
{ ## Again, construction is a work around for nlme 3-1.52
} else { ## Again, construction is a work around for nlme 3-1.52
if (wisvf) stop("weights must be like glm weights for generalized case")
if (verbosePQL) cat("\n Maximum number of PQL iterations: ",niterPQL,"\n")
eval(parse(text=paste("ret$lme<-gammPQL(",deparse(fixed.formula),
......@@ -1366,8 +1376,7 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
"correlation=correlation,control=control,",
"weights=weights,niter=niterPQL,verbose=verbosePQL)",sep="")))
G$y <- ret$lme$y ## makes sure that binomial response is returned as a vector!
##ret$lme<-glmmPQL(fixed.formula,random=rand,data=mf,family=family,correlation=correlation,
## control=control,niter=niterPQL,verbose=verbosePQL)
}
### .... fitting finished
......@@ -1532,6 +1541,11 @@ gamm <- function(formula,random=NULL,correlation=NULL,family=gaussian(),data=lis
object$edf<-rowSums(Vb*t(XVX))
object$df.residual <- length(object$y) - sum(object$edf)
#if (!is.null(ret$lme$aic)) { ## finish the conditional AIC (only happens if no correlation)
# object$aic <- ret$lme$aic + sum(object$edf) ## requires edf for r.e. as well!
# ret$lme$aic<- NULL
#}
object$sig2 <- ret$lme$sigma^2
if (lme.used) { object$method <- paste("lme.",method,sep="")}
else { object$method <- "PQL"}
......
......@@ -5,7 +5,8 @@ Rrank <- function(R,tol=.Machine$double.eps^.9) {
## Finds rank of upper triangular matrix R, by estimating condition
## number of upper rank by rank block, and reducing rank until this is
## acceptably low... assumes R pivoted
rank <- m <- ncol(R)
m <- nrow(R)
rank <- min(m,ncol(R))
ok <- TRUE
while (ok) {
Rcond <- .C(C_R_cond,R=as.double(R),r=as.integer(m),c=as.integer(rank),
......@@ -250,12 +251,13 @@ interpret.gam0 <- function(gf,textra=NULL,extra.special=NULL)
## loadNamespace('mgcv'); k <- 10; mgcv::interpret.gam(y~s(x,k=k)) fails (can't find s)
## eval(parse(text=terms[i]),envir=p.env,enclos=loadNamespace('mgcv')) fails??
## following may supply namespace of mgcv explicitly if not on search path...
## If 's' etc are masked then we can fail even if mgcv on search path
#if (mgcvat) st <- eval(parse(text=terms[i]),envir=p.env) else {
st <- try(eval(parse(text=terms[i]),envir=p.env),silent=TRUE)
if (inherits(st,"try-error")) st <-
## If 's' etc are masked then we can fail even if mgcv on search path, hence paste
## of explicit mgcv reference into first attempt...
st <- try(eval(parse(text=paste("mgcv::",terms[i],sep="")),envir=p.env),silent=TRUE)
if (inherits(st,"try-error")) st <-
eval(parse(text=terms[i]),enclos=p.env,envir=mgcvns)
#}
if (!is.null(textra)) { ## modify the labels on smooths with textra
pos <- regexpr("(",st$lab,fixed=TRUE)[1]
st$label <- paste(substr(st$label,start=1,stop=pos-1),textra,
......@@ -1206,8 +1208,13 @@ gam.setup <- function(formula,pterms,
}
G$P[qrx$pivot,] <- G$P
}
if (G$nsdf>0) G$cmX[-(1:G$nsdf)] <- 0 ## zero the smooth parts here
else G$cmX <- G$cmX * 0
## cmX relates to computation of CIs incorportating uncertainty about the mean
## It may make more sense to incorporate all uncertainty about the mean,
## rather than just the uncertainty in the fixed effects mean. This means
## incorporating the mean of random effects and unconstrained smooths. Hence
## comment out the following.
#if (G$nsdf>0) G$cmX[-(1:G$nsdf)] <- 0 ## zero the smooth parts here
#else G$cmX <- G$cmX * 0
G$X <- X;rm(X)
n.p <- ncol(G$X)
# deal with penalties
......@@ -1436,10 +1443,10 @@ gam.outer <- function(lsp,fscale,family,control,method,optimizer,criterion,scale
if (optimizer[2]%in%c("nlm.fd")) .Deprecated(msg=paste("optimizer",optimizer[2],"is deprecated, please use newton or bfgs"))
if (optimizer[1]=="efs" && !inherits(family,"general.family")) {
warning("Extended Fellner Schall only implemented for general families")
optimizer <- c("outer","newton")
}
# if (optimizer[1]=="efs" && !inherits(family,"general.family")) {
# warning("Extended Fellner Schall only implemented for general families")
# optimizer <- c("outer","newton")
# }
if (length(lsp)==0) { ## no sp estimation to do -- run a fit instead
optimizer[2] <- "no.sps" ## will cause gam2objective to be called, below
......@@ -1469,8 +1476,15 @@ gam.outer <- function(lsp,fscale,family,control,method,optimizer,criterion,scale
if (optimizer[1]=="efs"&& optimizer[2] != "no.sps" ) { ## experimental extended efs
##warning("efs is still experimental!")
object <- efsud(x=G$X,y=G$y,lsp=lsp,Sl=G$Sl,weights=G$w,offset=G$offxset,family=family,
if (inherits(family,"general.family")) {
object <- efsud(x=G$X,y=G$y,lsp=lsp,Sl=G$Sl,weights=G$w,offset=G$offxset,family=family,
control=control,Mp=G$Mp,start=start)
} else {
family <- fix.family.ls(family)
object <- efsudr(x=G$X,y=G$y,lsp=lsp,Eb=G$Eb,UrS=G$UrS,weights=G$w,family=family,offset=G$offset,
start=start,
U1=G$U1, intercept = TRUE,scale=scale,Mp=G$Mp,control=control,n.true=G$n.true,...)
}
object$gcv.ubre <- object$REML
} else if (optimizer[2]=="newton"||optimizer[2]=="bfgs"){ ## the gam.fit3 method
if (optimizer[2]=="bfgs")
......@@ -1596,6 +1610,7 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,start=NU
}
if (!optimizer[1]%in%c("perf","outer","efs")) stop("unknown optimizer")
if (optimizer[1]=="efs") method <- "REML"
if (!method%in%c("GCV.Cp","GACV.Cp","REML","P-REML","ML","P-ML")) stop("unknown smoothness selection criterion")
G$family <- fix.family(G$family)
G$rS <- mini.roots(G$S,G$off,ncol(G$X),G$rank)
......@@ -2085,7 +2100,7 @@ gam.control <- function (nthreads=1,irls.reg=0.0,epsilon = 1e-7, maxit = 200,
mgcv.tol=1e-7,mgcv.half=15,trace =FALSE,
rank.tol=.Machine$double.eps^0.5,
nlm=list(),optim=list(),newton=list(),outerPIsteps=0,
idLinksBases=TRUE,scalePenalty=TRUE,
idLinksBases=TRUE,scalePenalty=TRUE,efs.lspmax=15,efs.tol=.1,
keepData=FALSE,scale.est="fletcher",edge.correct=FALSE)
# Control structure for a gam.
# irls.reg is the regularization parameter to use in the GAM fitting IRLS loop.
......@@ -2138,12 +2153,13 @@ gam.control <- function (nthreads=1,irls.reg=0.0,epsilon = 1e-7, maxit = 200,
# and optim defaults
if (is.null(optim$factr)) optim$factr <- 1e7
optim$factr <- abs(optim$factr)
if (efs.tol<=0) efs.tol <- .1
list(nthreads=round(nthreads),irls.reg=irls.reg,epsilon = epsilon, maxit = maxit,
trace = trace, mgcv.tol=mgcv.tol,mgcv.half=mgcv.half,
rank.tol=rank.tol,nlm=nlm,
optim=optim,newton=newton,outerPIsteps=outerPIsteps,
idLinksBases=idLinksBases,scalePenalty=scalePenalty,
idLinksBases=idLinksBases,scalePenalty=scalePenalty,efs.lspmax=efs.lspmax,efs.tol=efs.tol,
keepData=as.logical(keepData[1]),scale.est=scale.est,edge.correct=edge.correct)
}
......@@ -2509,7 +2525,7 @@ get.na.action <- function(na.action) {
predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclude=NULL,
block.size=NULL,newdata.guaranteed=FALSE,na.action=na.pass,
unconditional=FALSE,...) {
unconditional=FALSE,iterms.type=NULL,...) {
# This function is used for predicting from a GAM. 'object' is a gam object, newdata a dataframe to
# be used in prediction......
......@@ -2663,7 +2679,12 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclu
msg <- paste(paste(levn[!levn%in%levm],collapse=", "),"not in original fit",collapse="")
stop(msg)
}
newdata[[i]] <- factor(newdata[[i]],levels=levm) # set prediction levels to fit levels
## set prediction levels to fit levels...
if (is.matrix(newdata[[i]])) {
dum <- factor(newdata[[i]],levels=levm)
dim(dum) <- dim(newdata[[i]])
newdata[[i]] <- dum
} else newdata[[i]] <- factor(newdata[[i]],levels=levm)
}
if (type=="newdata") return(newdata)
......@@ -2851,6 +2872,7 @@ predict.gam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,exclu
if (type=="iterms"&& attr(object$smooth[[k]],"nCons")>0) { ## termwise se to "carry the intercept
## some general families, add parameters after cmX created, which are irrelevant to cmX...
if (length(object$cmX) < ncol(X)) object$cmX <- c(object$cmX,rep(0,ncol(X)-length(object$cmX)))
if (!is.null(iterms.type)&&iterms.type==2) object$cmX[-(1:object$nsdf)] <- 0 ## variability of fixed effects mean only
X1 <- matrix(object$cmX,nrow(X),ncol(X),byrow=TRUE)
meanL1 <- object$smooth[[k]]$meanL1
if (!is.null(meanL1)) X1 <- X1 / meanL1
......
## R plotting routines for the package mgcv (c) Simon Wood 2000-2017
# R plotting routines for the package mgcv (c) Simon Wood 2000-2017
## With contributions from Henric Nilsson
......@@ -1247,6 +1247,7 @@ plot.gam <- function(x,residuals=FALSE,rug=NULL,se=TRUE,pages=0,select=NULL,scal
## test whether mean variability to be added to variability (only for centred terms)
if (seWithMean && attr(x$smooth[[i]],"nCons")>0) {
if (length(x$cmX) < ncol(x$Vp)) x$cmX <- c(x$cmX,rep(0,ncol(x$Vp)-length(x$cmX)))
if (seWithMean==2) x$cmX[-(1:x$nsdf)] <- 0 ## variability of fixed effects mean only
X1 <- matrix(x$cmX,nrow(P$X),ncol(x$Vp),byrow=TRUE)
meanL1 <- x$smooth[[i]]$meanL1
if (!is.null(meanL1)) X1 <- X1 / meanL1
......@@ -1297,17 +1298,6 @@ plot.gam <- function(x,residuals=FALSE,rug=NULL,se=TRUE,pages=0,select=NULL,scal
} else
{ ppp<-1;oldpar<-par()}
if ((pages==0&&prod(par("mfcol"))<n.plots&&dev.interactive())||
pages>1&&dev.interactive()) ask <- TRUE else ask <- FALSE
if (!is.null(select)) {
ask <- FALSE
}
if (ask) {
oask <- devAskNewPage(TRUE)
on.exit(devAskNewPage(oask))
}
#####################################
## get a common scale, if required...
......@@ -1347,11 +1337,28 @@ plot.gam <- function(x,residuals=FALSE,rug=NULL,se=TRUE,pages=0,select=NULL,scal
## now plot smooths, by calling plot methods with plot data...
##############################################################
if ((pages==0&&prod(par("mfcol"))<n.plots&&dev.interactive())||
pages>1&&dev.interactive()) ask <- TRUE else ask <- FALSE
if (!is.null(select)) {
ask <- FALSE
}
# if (ask) { ## asks before plotting
# oask <- devAskNewPage(TRUE)
# on.exit(devAskNewPage(oask))
# }
if (m>0) for (i in 1:m) if (pd[[i]]$plot.me&&(is.null(select)||i==select)) {
plot(x$smooth[[i]],P=pd[[i]],partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,n3=n3,
pers=pers,theta=theta,phi=phi,jit=jit,xlab=xlab,ylab=ylab,main=main,
ylim=ylim,xlim=xlim,too.far=too.far,shade=shade,shade.col=shade.col,
shift=shift,trans=trans,by.resids=by.resids,scheme=scheme[i],...)
if (ask) { ## this is within loop so we don't get asked before it's necessary
oask <- devAskNewPage(TRUE)
on.exit(devAskNewPage(oask))
ask <- FALSE ## only need to do this once
}
} ## end of smooth plotting loop
......
......@@ -332,7 +332,7 @@ get.var <- function(txt,data,vecMat = TRUE)
ismat <- FALSE
} else ismat <- TRUE
} else ismat <- FALSE
if (vecMat&&is.matrix(x)) x <- as.numeric(x)
if (vecMat&&is.matrix(x)) x <- x[1:prod(dim(x))] ## modified from x <- as.numeric(x) to allow factors
if (ismat) attr(x,"matrix") <- TRUE
x
} ## get.var
......@@ -2545,7 +2545,14 @@ smooth.construct.mrf.smooth.spec <- function(object, data, knots) {
if (is.null(object$xt$nb)) { ## no neighbour list... construct one
if (is.null(object$xt$polys)) stop("no spatial information provided!")
object$xt$nb <- pol2nb(object$xt$polys)$nb
} ## now have a neighbour list
} else if (!is.numeric(object$xt$nb[[1]])) { ## user has (hopefully) supplied names not indices
nb.names <- names(object$xt$nb)
for (i in 1:length(nb.names)) {
object$xt$nb[[i]] <- which(nb.names %in% object$xt$nb[[i]])
}
}
## now have a neighbour list
a.name <- names(object$xt$nb)
if (all.equal(sort(a.name),sort(levels(k)))!=TRUE)
stop("mismatch between nb/polys supplied area names and data area names")
......@@ -3656,10 +3663,10 @@ smoothCon <- function(object,data,knots=NULL,absorb.cons=FALSE,scale.penalty=TRU
sml[[1]]$X <- as.numeric(by)*sm$X ## normal `by' handling
## Now do the summation stuff....
ind <- 1:n
X <- sml[[1]]$X[ind,]
X <- sml[[1]]$X[ind,,drop=FALSE]
for (i in 2:q) {
ind <- ind + n
X <- X + sml[[1]]$X[ind,]
X <- X + sml[[1]]$X[ind,,drop=FALSE]
}
sml[[1]]$X <- X
if (!is.null(offs)) { ## deal with any term specific offset (i.e. sum it too)
......@@ -3739,15 +3746,15 @@ smoothCon <- function(object,data,knots=NULL,absorb.cons=FALSE,scale.penalty=TRU
if (length(sm$S)>0)
for (l in 1:length(sm$S)) # some smooths have > 1 penalty
{ ZSZ <- sml[[i]]$S[[l]]
ZSZ[indi[1:nz],]<-qr.qty(qrc,sml[[i]]$S[[l]][indi,,drop=FALSE])[(nc+1):nx,]
if (nz>0) ZSZ[indi[1:nz],]<-qr.qty(qrc,sml[[i]]$S[[l]][indi,,drop=FALSE])[(nc+1):nx,]
ZSZ <- ZSZ[-indi[(nz+1):nx],]
ZSZ[,indi[1:nz]]<-t(qr.qty(qrc,t(ZSZ[,indi,drop=FALSE]))[(nc+1):nx,])
if (nz>0) ZSZ[,indi[1:nz]]<-t(qr.qty(qrc,t(ZSZ[,indi,drop=FALSE]))[(nc+1):nx,])
sml[[i]]$S[[l]] <- ZSZ[,-indi[(nz+1):nx],drop=FALSE] ## Z'SZ
## ZSZ<-qr.qty(qrc,sm$S[[l]])[(j+1):k,]
## sml[[i]]$S[[l]]<-t(qr.qty(qrc,t(ZSZ))[(j+1):k,]) ## Z'SZ
}
sml[[i]]$X[,indi[1:nz]]<-t(qr.qty(qrc,t(sml[[i]]$X[,indi,drop=FALSE]))[(nc+1):nx,])
if (nz>0) sml[[i]]$X[,indi[1:nz]]<-t(qr.qty(qrc,t(sml[[i]]$X[,indi,drop=FALSE]))[(nc+1):nx,])
sml[[i]]$X <- sml[[i]]$X[,-indi[(nz+1):nx]]
## sml[[i]]$X<-t(qr.qty(qrc,t(sml[[i]]$X))[(j+1):k,]) ## form XZ
attr(sml[[i]],"qrc") <- qrc
......@@ -4025,11 +4032,11 @@ PredictMat <- function(object,data,n=nrow(data))
nc <- j;nz <- nx - nc
if (sum(is.na(X))) {
ind <- !is.na(rowSums(X))
X[ind,indi[1:nz]]<-t(qr.qty(qrc,t(X[ind,indi,drop=FALSE]))[(nc+1):nx,])
if (nz>0) X[ind,indi[1:nz]]<-t(qr.qty(qrc,t(X[ind,indi,drop=FALSE]))[(nc+1):nx,])
X <- X[,-indi[(nz+1):nx]]
X[!ind,] <- NA
} else {
X[,indi[1:nz]]<-t(qr.qty(qrc,t(X[,indi,drop=FALSE]))[(nc+1):nx,,drop=FALSE])
if (nz>0) X[,indi[1:nz]]<-t(qr.qty(qrc,t(X[,indi,drop=FALSE]))[(nc+1):nx,,drop=FALSE])
X <- X[,-indi[(nz+1):nx]]
}
}
......
File mode changed from 100644 to 100755
......@@ -218,9 +218,9 @@ terms less the number of spline terms).
This routine is less stable than `gam' for the same dataset.
The negbin family is only supported for the *known theta* case.
With \code{discrete=TRUE}, \code{te} terms are efficiently computed, but \code{t2} are not.
The negbin family is only supported for the *known theta* case.
}
\seealso{\code{\link{mgcv.parallel}},
......@@ -236,9 +236,7 @@ The negbin family is only supported for the *known theta* case.
library(mgcv)
## See help("mgcv-parallel") for using bam in parallel
## Some examples are marked 'Not run' purely to keep
## checking load on CRAN down. Sample sizes are small for
## the same reason.
## Sample sizes are small for fast run times.
set.seed(3)
dat <- gamSim(1,n=25000,dist="normal",scale=20)
......@@ -263,6 +261,14 @@ system.time(b1 <- bam(y ~ s(x0,bs=bs)+s(x1,bs=bs)+s(x2,bs=bs,k=k),
data=dat,family=poisson()))
b1
## Similar using faster discrete method...
\donttest{
system.time(b2 <- bam(y ~ s(x0,bs=bs,k=k)+s(x1,bs=bs,k=k)+s(x2,bs=bs,k=k)+
s(x3,bs=bs,k=k),data=dat,family=poisson(),discrete=TRUE))
b2
}
}
......
......@@ -60,8 +60,8 @@ tdpois <- function(dat,event="z",et="futime",t="day",status="status1",
dap[1:(start-1),]
} ## tdpois
## The following takes too long for CRAN checking
## (but typically takes a minute or less).
## The following typically takes a minute or less...
\donttest{
## Convert pbcseq to equivalent Poisson form...
pbcseq$status1 <- as.numeric(pbcseq$status==2) ## death indicator
......
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
......@@ -87,7 +87,7 @@ parameter estimation criterion (given by \code{method}). \code{"perf"} (deprecat
for the more stable direct approach. \code{"outer"} can use several alternative optimizers, specified in the
second element of \code{optimizer}: \code{"newton"} (default), \code{"bfgs"}, \code{"optim"}, \code{"nlm"}
and \code{"nlm.fd"} (the latter is based entirely on finite differenced derivatives and is very slow). \code{"efs"}
for the extended Fellner Schall method - only currently available for general families.}
for the extended Fellner Schall method of Wood and Fasiolo (2017).}
\item{scale}{ If this is positive then it is taken as the known scale parameter. Negative signals that the
scale parameter is unknown. 0 signals that the scale parameter is 1 for Poisson and binomial and unknown otherwise.
......@@ -298,6 +298,10 @@ generalized additive mixed models. Biometrics 62(4):1025-1036
Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman
and Hall/CRC Press.
Wood, S.N. and M. Fasiolo (2017) A generalized Fellner-Schall method for smoothing
parameter optimization with application to Tweedie location, scale and shape models.
Biometrics 73 (4), 1071-1081
Wood S.N., F. Scheipl and J.J. Faraway (2012) Straightforward intermediate rank tensor product smoothing
in mixed models. Statistical Computing.
......
......@@ -12,10 +12,11 @@ choise of fitting method, see \code{\link{gam}} arguments \code{method} and \cod
\usage{
gam.control(nthreads=1,irls.reg=0.0,epsilon = 1e-07, maxit = 200,
mgcv.tol=1e-7,mgcv.half=15, trace = FALSE,
rank.tol=.Machine$double.eps^0.5,
nlm=list(),optim=list(),newton=list(),
outerPIsteps=0,idLinksBases=TRUE,scalePenalty=TRUE,
keepData=FALSE,scale.est="fletcher",edge.correct=FALSE)
rank.tol=.Machine$double.eps^0.5,nlm=list(),
optim=list(),newton=list(),outerPIsteps=0,
idLinksBases=TRUE,scalePenalty=TRUE,efs.lspmax=15,
efs.tol=.1,keepData=FALSE,scale.est="fletcher",
edge.correct=FALSE)
}
\arguments{
\item{nthreads}{Some parts of some smoothing parameter selection methods (e.g. REML) can use some
......@@ -72,6 +73,12 @@ of the penalty matrices of a smooth relative to its model matrix. This option re
the penalty matrices to accomodate this problem. Probably should be set to \code{FALSE}
if you are linking smoothing parameters but have set \code{idLinkBases} to \code{FALSE}.}
\item{efs.lspmax}{maximum log smoothing parameters to allow under extended Fellner Schall
smoothing parameter optimization.}
\item{efs.tol}{change in REML to count as negligible when testing for EFS convergence. If the
step is small and the last 3 steps led to a REML change smaller than this, then stop.}
\item{keepData}{Should a copy of the original \code{data} argument be kept in the \code{gam}
object? Strict compatibility with class \code{glm} would keep it, but it wastes space to
do so. }
......
......@@ -24,7 +24,9 @@ The fitted values for this family will be a three column matrix. The first colum
This family does not produce a null deviance. Note that the distribution for \eqn{\xi=0}{xi=0} is approximated by setting \eqn{\xi}{xi} to a small number.
The derivative system code for this family is mostly auto-generated, and the family is still somewhat experimental.
The derivative system code for this family is mostly auto-generated, and the family is still somewhat experimental.
The GEV distribution is rather challenging numerically, and for small datasets or poorly fitting models imporoved numerical robustness may be obtained by using the extended Fellner-Schall method of Wood and Fasiolo (2017) for smoothing parameter estimation. See examples.
}
\references{
......@@ -36,6 +38,8 @@ model selection for general smooth models.
Journal of the American Statistical Association 111, 1548-1575
\url{http://dx.doi.org/10.1080/01621459.2016.1180986}
Wood, S.N. and M. Fasiolo (2017) A generalized Fellner-Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models. Biometrics 73(4): 1071-1081.
\url{http://dx.doi.org/10.1111/biom.12666}
}
......@@ -66,6 +70,11 @@ dat <- data.frame(y,x0,x1,x2);pairs(dat)
## fit model....
b <- gam(list(y~s(x2),~s(x0),~s(x1)),family=gevlss,data=dat)
## same fit using the extended Fellner-Schall method which
## can provide improved numerical robustness...
b <- gam(list(y~s(x2),~s(x0),~s(x1)),family=gevlss,data=dat,
optimizer="efs")
## plot and look at residuals...
plot(b,pages=1,scale=0)
summary(b)
......
File mode changed from 100644 to 100755
......@@ -147,10 +147,21 @@ dat <- gamSim(1,n=n,dist="normal",scale=2)
## regular gam fit for comparison...
b0 <- gam(y~s(x0)+s(x1) + s(x2)+s(x3),data=dat,method="REML")
## Set directory and file name for file containing jags code.
## In real use you would *never* use tempdir() for this. It is
## only done here to keep CRAN happy, and avoid any chance of
## an accidental overwrite. Instead you would use
## setwd() to set an appropriate working directory in which
## to write the file, and just set the file name to what you
## want to call it (e.g. "test.jags" here).
jags.file <- paste(tempdir(),"/test.jags",sep="")
## Set up JAGS code and data. In this one might want to diagonalize
## to use conjugate samplers. Usually call 'setwd' first, to set
## directory in which model file ("test.jags") will be written.
jd <- jagam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,file="test.jags",
jd <- jagam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,file=jags.file,
sp.prior="gamma",diagonalize=TRUE)
## In normal use the model in "test.jags" would now be edited to add
......@@ -159,7 +170,7 @@ jd <- jagam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,file="test.jags",
\dontrun{
require(rjags)
load.module("glm") ## improved samplers for GLMs often worth loading
jm <-jags.model("test.jags",data=jd$jags.data,inits=jd$jags.ini,n.chains=1)
jm <-jags.model(jags.file,data=jd$jags.data,inits=jd$jags.ini,n.chains=1)
list.samplers(jm)
sam <- jags.samples(jm,c("b","rho","scale"),n.iter=10000,thin=10)
jam <- sim2jam(sam,jd$pregam)
......@@ -178,7 +189,7 @@ dat <- gamSim(1,n=n,dist="normal",scale=2)
scale <- .5; Ey <- exp(dat$f/2)
dat$y <- rgamma(n,shape=1/scale,scale=Ey*scale)
jd <- jagam(y~s(x0)+te(x1,x2)+s(x3),data=dat,family=Gamma(link=log),
file="test.jags",sp.prior="log.uniform")
file=jags.file,sp.prior="log.uniform")
## In normal use the model in "test.jags" would now be edited to add
## the non-standard stochastic elements that require use of JAGS....
......@@ -189,7 +200,7 @@ require(rjags)
## models are still not fully repeatable (JAGS 4 should fix this)
jd$jags.ini$.RNG.name <- "base::Mersenne-Twister" ## setting RNG
jd$jags.ini$.RNG.seed <- 6 ## how to set RNG seed
jm <-jags.model("test.jags",data=jd$jags.data,inits=jd$jags.ini,n.chains=1)
jm <-jags.model(jags.file,data=jd$jags.data,inits=jd$jags.ini,n.chains=1)
list.samplers(jm)
sam <- jags.samples(jm,c("b","rho","scale","mu"),n.iter=10000,thin=10)
jam <- sim2jam(sam,jd$pregam)
......
\name{k.check}
\alias{k.check}
\title{Checking smooth basis dimension }
\description{ Takes a fitted \code{gam} object produced by \code{gam()} and runs
diagnostic tests of whether the basis dimension choises are adequate.
}
\usage{
k.check(b, subsample=5000, n.rep=400)
}
\arguments{
\item{b}{a fitted \code{gam} object as produced by \code{\link{gam}()}.}
\item{subsample}{above this number of data, testing uses a random sub-sample of data of this size.}
\item{n.rep}{how many re-shuffles to do to get p-value for k testing.}
}
\value{A matrix contaning the output of the tests described above.}
\details{
The test of whether the basis dimension for a smooth is adequate (Wood, 2017, section 5.9) is based on computing an estimate of the residual variance based on differencing residuals that are near neighbours according to the (numeric) covariates of the smooth. This estimate divided by the residual variance is the \code{k-index} reported. The further below 1 this is, the more likely it is that there is missed pattern left in the residuals. The \code{p-value} is computed by simulation: the residuals are randomly re-shuffled \code{n.rep} times to obtain the null distribution of the differencing variance estimator, if there is no pattern in the residuals. For models fitted to more than \code{subsample} data, the tests are based of \code{subsample} randomly sampled data. Low p-values may indicate that the basis dimension, \code{k}, has been set too low, especially if the reported \code{edf} is close to \code{k\'}, the maximum possible EDF for the term. Note the disconcerting fact that if the test statistic itself is based on random resampling and the null is true, then the associated p-values will of course vary widely from one replicate to the next. Currently smooths of factor variables are not supported and will give an \code{NA} p-value.
Doubling a suspect \code{k} and re-fitting is sensible: if the reported \code{edf} increases substantially then you may have been missing something in the first fit. Of course p-values can be low for reasons other than a too low \code{k}. See \code{\link{choose.k}} for fuller discussion.
}
\references{
Wood S.N. (2017) Generalized Additive Models: An Introduction with R (2nd edition). Chapman
and Hall/CRC Press.
\url{http://www.maths.bris.ac.uk/~sw15190/}
}
\author{ Simon N. Wood \email{simon.wood@r-project.org}}
\seealso{ \code{\link{choose.k}}, \code{\link{gam}}, \code{\link{gam.check}}}
\examples{
library(mgcv)
set.seed(0)
dat <- gamSim(1,n=200)
b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
plot(b,pages=1)
k.check(b)
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ...
......@@ -31,11 +31,12 @@ variable then the smooth will automatically be centred in order to ensure identi
will not be. Note also that for centred smooths it can be worth replacing the constant term in the model with \code{rowSums(L)}
in order to ensure that predictions are automatically on the right scale.
Note that \code{\link{predict.gam}} can accept matrix predictors for prediction with such terms, in which case its
\code{\link{predict.gam}} can accept matrix predictors for prediction with such terms, in which case its
\code{newdata} argument will need to be a list. However when