Commit 86e4d7fd authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-0

parent e62c8468
** denotes quite substantial/important changes
*** denotes really big changes
1.7-29
1.8-0
* Single character change to Makevars file so that openMP multi-threading
actually works.
*** Cox Proportional Hazard family 'cox.ph' added as example of general
penalized likelihood families now useable with 'gam'.
*** 'ocat', 'tw', 'nb', 'betar', 'ziP' and 'scat' families added for
ordered categorical data, Tweedie with estimation of 'p', negative binomial
with (fast) estimation of 'theta', beta regression for proportions, simple
zero inflated Poisson regression and heavy tailed regression with scaled t
distribution. These are all examples of 'extended families' now useable
with 'gam'.
*** 'gaulss' and 'ziplss' families, implementing models with multiple linear
predictors. For gaulss there is a linear predictor for the Gaussian mean
and another for the standard deviation. For ziplss there is a linear
predictor controlling `presence' and another controlling
the Poisson parameter, given presence.
*** 'mvn' family for multivariate normal additive models.
** AIC computation changed for bam and gam models estimated by REML/ML
to account for smoothing parameter uncertainty in degrees of freedom
term.
* With REML/ML smoothness selection in gam/bam an extra covariance matrix 'Vc'
is now computed which allows for smoothing parameter uncertainty. See
the 'unconditional' arguments to 'predict.gam' and 'plot.gam' to use this.
* 'gam.vcomp' bug fix. Computed intervals for families with fixed scale
parameter were too wide.
* gam now defaults to the Pearson estimator of the scale parameter to avoid
poor scale estimates in the quasipoisson case with low counts (and possibly
elsewhere). Gaussian, Poisson and binomial inference invariant to change.
Thanks to Greg Dropkin, for reporting the issue.
* Polish translation added thanks to Lukasz Daniel.
* gam.fit3 now forces eta and mu to be consistent with coef and valid on
return (previously could happen that if step halving was used in final
iteration then eta or mu could be invalid, e.g. when using identity link
with non-negative data)
* gam.fit3 now bases its convergence criteria on grad deviance w.r.t. model
coefs, rather than changes in model coefs. This prevents problems when
there is rank deficiency but different coefs get dropped at different
iterations. Thanks to Kristynn Sullivan.
* If mgcv is not on the search path then interpret.gam now tries to
evaluate in namespace of mgcv with environment of formula as enclosing
environment, if evaluation in the environment of the formula fails.
* bug fix to sos plotting method so that it now works with 'by' variables.
* 'plot.gam' now weights partial residuals by *normalized* square root
iterative weights so that the average weight is 1 and the residuals
should have constant variance if all is ok.
* 'pcls' now reports if the initial point is not feasible.
* 'print.gam' and 'summary.gam' now report the rank of the model if it is
rank deficient. 'gam.check' reports the model rank whenever it is
available.
* fix of bug in 'k.check' called by 'gam.check' that gave an error for
smooths with by variables.
* predict.gam now checks that factors in newdata do not contain more
levels than those used in fitting.
* predict.gam could fail for type "terms" with no intercept - fixed.
* 'bfgs' now uses a finite difference approximation for the initial inverse
Hessian.
1.7-28
......
Package: mgcv
Version: 1.7-29
Version: 1.8-0
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
......@@ -10,11 +10,11 @@ Description: Routines for GAMs and other generalized ridge regression
Priority: recommended
Depends: R (>= 2.14.0), nlme (>= 3.1-64)
Imports: methods, stats, graphics, Matrix
Suggests: splines, parallel
Suggests: splines, parallel, survival, MASS
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2014-03-30 21:09:39 UTC; sw283
Packaged: 2014-06-25 05:49:03 UTC; sw283
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2014-03-31 08:57:05
Date/Publication: 2014-06-25 09:05:33
This diff is collapsed.
useDynLib(mgcv, .registration = TRUE, .fixes = "C_")
export(anova.gam, bam, bam.update, concurvity, cSplineDes,
export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
exclude.too.far,extract.lme.cov, extract.lme.cov2,
formXtViX, full.score, formula.gam,fixDependence,fix.family.link,
fix.family.var, fix.family.ls, fix.family.qf,fix.family.rd,
fs.test,fs.boundary,gam, gam2derivative,
gam2objective,
gamm, gam.check, gam.control,gam.fit3,
gam.fit, gam.outer,gam.vcomp, gamSim , influence.gam,
gam.fit, gam.outer,gam.vcomp, gamSim ,
gaulss,
influence.gam,
in.out,inSide,interpret.gam,
gam.side,
get.var,ldTweedie,
initial.sp,logLik.gam,ls.size,
magic, magic.post.proc, model.matrix.gam,
mono.con, mroot, negbin, new.name,
mono.con, mroot, mvn, nb, negbin, new.name,
notExp,notExp2,notLog,notLog2,pcls,null.space.dimension,
ocat,
pen.edf,pdIdnot,pdTens,
place.knots, plot.gam, polys.plot,print.anova.gam,
print.gam,print.summary.gam,predict.gam,
......@@ -36,7 +39,7 @@ export(anova.gam, bam, bam.update, concurvity, cSplineDes,
Predict.matrix.t2.smooth,
qq.gam,
residuals.gam,rig,rTweedie,
Rrank,s,
Rrank,s,scat,
slanczos,
smoothCon,smooth.construct,smooth.construct2,
smooth.construct.cc.smooth.spec,
......@@ -59,25 +62,16 @@ export(anova.gam, bam, bam.update, concurvity, cSplineDes,
summary.gam,sp.vcov,
spasm.construct,spasm.sp,spasm.smooth,
t2,te,ti,tensor.prod.model.matrix,tensor.prod.penalties,
Tweedie,uniquecombs, vcov.gam, vis.gam)
Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
importFrom(grDevices,cm.colors,gray,heat.colors,terrain.colors,topo.colors)
importFrom(graphics,axis,box,contour,hist,lines,mtext, par, persp,plot,points,
polygon,strheight,strwidth,text)
#importFrom(stats,.checkMFClasses,.getXlevels, as.formula, anova,anova.glmlist,
# coef,cooks.distance,cor,delete.response,family,fitted,formula,
# gaussian,glm,influence,lm,logLik,median,model.frame,model.matrix,
# model.offset,na.pass,napredict,naresid,nlm,optim,pchisq,pf,pnorm,pt,
# predict,printCoefmat ,quantile,qqnorm, reformulate,residuals,rnorm,
# 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,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)
......@@ -111,24 +105,44 @@ S3method(pdMatrix,pdIdnot)
S3method(solve,pdIdnot)
S3method(summary,pdIdnot)
S3method(fix.family.link,family)
S3method(fix.family.link,extended.family)
S3method(fix.family.link,general.family)
S3method(smooth.construct,ad.smooth.spec)
S3method(smooth.construct,ps.smooth.spec)
S3method(smooth.construct,cp.smooth.spec)
S3method(smooth.construct, cc.smooth.spec)
S3method(smooth.construct,cp.smooth.spec)
S3method(smooth.construct, cr.smooth.spec)
S3method(smooth.construct, cs.smooth.spec)
S3method(smooth.construct,ds.smooth.spec)
S3method(smooth.construct, fs.smooth.spec)
S3method(smooth.construct, mrf.smooth.spec)
S3method(smooth.construct,ps.smooth.spec)
S3method(smooth.construct, re.smooth.spec)
S3method(smooth.construct,so.smooth.spec)
S3method(smooth.construct,sos.smooth.spec)
S3method(smooth.construct, tp.smooth.spec)
S3method(smooth.construct, tensor.smooth.spec)
S3method(smooth.construct, cs.smooth.spec)
S3method(smooth.construct, t2.smooth.spec)
S3method(smooth.construct, ts.smooth.spec)
S3method(Predict.matrix,cs.smooth)
S3method(Predict.matrix,ts.smooth)
S3method(Predict.matrix,pspline.smooth)
S3method(Predict.matrix,cr.smooth)
S3method(Predict.matrix,cs.smooth)
S3method(Predict.matrix,cyclic.smooth)
S3method(Predict.matrix,tensor.smooth)
S3method(Predict.matrix,cpspline.smooth)
S3method(Predict.matrix,duchon.spline)
S3method(Predict.matrix,fs.interaction)
S3method(Predict.matrix,mrf.smooth)
S3method(Predict.matrix,pspline.smooth)
S3method(Predict.matrix,random.effect)
S3method(Predict.matrix,tprs.smooth)
S3method(Predict.matrix,ts.smooth)
S3method(Predict.matrix,tensor.smooth)
S3method(Predict.matrix,t2.smooth)
S3method(Predict.matrix,soap.film)
S3method(Predict.matrix,sf)
S3method(Predict.matrix,sw)
S3method(Predict.matrix,sos.smooth)
S3method(spasm.construct,cus)
S3method(spasm.sp,cus)
......
......@@ -32,10 +32,6 @@ chol2qr <- function(XX,Xy) {
## equivalent to qr update. Uses simple
## regularization if XX not +ve def.
#d <- diag(XX)
#ind <- d>0
#d[!ind] <- 1;d[ind] <- 1/sqrt(d[ind])
#XX <- t(d*t(d*XX)) ## diagonal pre-conditioning
XXeps <- norm(XX)*.Machine$double.eps^.9
n <- ncol(XX)
......@@ -44,13 +40,11 @@ chol2qr <- function(XX,Xy) {
R <- try(chol(XX+diag(n)*XXeps*i,pivot=TRUE),silent=TRUE) ## R'R = X'X
if (inherits(R,"try-error")) ok <- FALSE else {
ipiv <- piv <- attr(R,"pivot")
#f <- try(forwardsolve(t(R),(Xy*d)[piv]),silent=TRUE)
f <- try(forwardsolve(t(R),Xy[piv]),silent=TRUE)
if (inherits(f,"try-error")) ok <- FALSE
}
if (ok) {
ipiv[piv] <- 1:ncol(R)
#R <- t(t(R[,ipiv])/d)
R <- R[,ipiv]
break; ## success
}
......@@ -97,7 +91,6 @@ qr.up <- function(arg) {
dev <- 0
for (b in 1:arg$n.block) {
ind <- arg$start[b]:arg$stop[b]
##arg$G$model <- arg$mf[ind,]
X <- predict(arg$G,newdata=arg$mf[ind,],type="lpmatrix",newdata.guaranteed=TRUE,block.size=length(ind))
rownames(X) <- NULL
if (is.null(arg$coef)) eta1 <- arg$eta[ind] else eta1 <- drop(X%*%arg$coef) + arg$offset[ind]
......@@ -315,9 +308,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
dev <- 0
if (n.threads == 1) { ## use original serial update code
for (b in 1:n.block) {
ind <- start[b]:stop[b]
##G$model <- mf[ind,]
X <- predict(G,newdata=mf[ind,],type="lpmatrix",newdata.guaranteed=TRUE,block.size=length(ind))
rownames(X) <- NULL
if (is.null(coef)) eta1 <- eta[ind] else eta1 <- drop(X%*%coef) + offset[ind]
......@@ -386,7 +377,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
rss.extra <- qrx$y.norm2 - sum(qrx$f^2)
if (control$trace)
cat("Deviance =", dev, "Iterations -", iter,"\n")
gettextf("Deviance = %s Iterations - %d", dev, iter, domain = "R-mgcv")
if (!is.finite(dev)) stop("Non-finite deviance")
......@@ -420,7 +411,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
log.phi=log.phi,phi.fixed=scale>0,rss.extra=rss.extra,
nobs =nobs+nobs.extra,Mp=um$Mp)
res <- Sl.postproc(Sl,fit,um$undrop,qrx$R,cov=FALSE)
object <- list(coefficients=res$beta,
object <- list(coefficients=res$beta,db.drho=fit$d1b,
gcv.ubre=fit$reml,mgcv.conv=list(iter=fit$iter,
message=fit$conv),rank=ncol(um$X),
Ve=NULL,scale.estimated = scale<=0,outer.info=fit$outer.info,
......@@ -457,7 +448,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
object$coefficients <- fit$b
object$edf <- post$edf
object$edf1 <- post$edf1
object$F <- post$F
##object$F <- post$F
object$full.sp <- fit$sp.full
object$gcv.ubre <- fit$score
object$hat <- post$hat
......@@ -476,8 +467,8 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
if (any(!is.finite(coef))) {
conv <- FALSE
warning("non-finite coefficients at iteration ",
iter)
warning(gettextf("non-finite coefficients at iteration %d",
iter))
break
}
} ## end fitting iteration
......@@ -486,10 +477,12 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
res <- Sl.postproc(Sl,fit,um$undrop,qrx$R,cov=TRUE,scale=scale)
object$edf <- res$edf
object$edf1 <- res$edf1
object$F <- res$F
object$edf2 <- res$edf2
##object$F <- res$F
object$hat <- res$hat
object$Vp <- res$Vp
object$Ve <- res$Ve
object$Vc <- res$Vc
}
if (!conv)
......@@ -596,7 +589,7 @@ bgam.fit2 <- function (G, mf, chunk.size, gp ,scale ,gamma,method, etastart = NU
rss.extra <- qrx$y.norm2 - sum(qrx$f^2)
if (control$trace)
cat("Deviance =", dev, "Iterations -", iter,"\n")
gettextf("Deviance = %s Iterations - %d", dev, iter, domain = "R-mgcv")
if (!is.finite(dev)) stop("Non-finite deviance")
......@@ -634,7 +627,7 @@ bgam.fit2 <- function (G, mf, chunk.size, gp ,scale ,gamma,method, etastart = NU
object$coefficients <- fit$b
object$edf <- post$edf
object$edf1 <- post$edf1
object$F <- post$F
#object$F <- post$F
object$full.sp <- fit$sp.full
object$gcv.ubre <- fit$score
object$hat <- post$hat
......@@ -1002,10 +995,12 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
log.phi=log.phi,phi.fixed=scale>0,rss.extra=rss.extra,
nobs =n,Mp=um$Mp)
res <- Sl.postproc(Sl,fit,um$undrop,qrx$R,cov=TRUE,scale=scale)
object <- list(coefficients=res$beta,edf=res$edf,edf1=res$edf1,F=res$F,
object <- list(coefficients=res$beta,edf=res$edf,edf1=res$edf1,edf2=res$edf2,##F=res$F,
db.drho=fit$d1b,
gcv.ubre=fit$reml,hat=res$hat,mgcv.conv=list(iter=fit$iter,
message=fit$conv),rank=ncol(um$X),
Ve=res$Ve,Vp=res$Vp,scale.estimated = scale<=0,outer.info=fit$outer.info,
Ve=res$Ve,Vp=res$Vp,Vc=res$Vc,
scale.estimated = scale<=0,outer.info=fit$outer.info,
optimizer=c("perf","newton"))
if (scale<=0) { ## get sp's and scale estimate
nsp <- length(fit$rho)
......@@ -1050,7 +1045,7 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
object$coefficients <- fit$b
object$edf <- post$edf
object$edf1 <- post$edf1
object$F <- post$F
##object$F <- post$F
object$full.sp <- fit$sp.full
object$gcv.ubre <- fit$score
object$hat <- post$hat
......@@ -1201,7 +1196,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
G$family <- family
G$terms<-terms;##G$pterms<-pterms
pvars <- all.vars(delete.response(terms))
G$pred.formula <- if (length(pvars)>0) reformulate(pvars) else NULL
G$pred.formula <- gp$pred.formula # if (length(pvars)>0) reformulate(pvars) else NULL
n <- nrow(mf)
......@@ -1302,7 +1297,7 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
object$offset <- G$offset
object$prior.weights <- G$w
object$pterms <- G$pterms
object$pred.formula <- G$pred.formula
object$pred.formula <- G$pred.formula
object$smooth <- G$smooth
object$terms <- G$terms
......@@ -1452,10 +1447,10 @@ bam.update <- function(b,data,chunk.size=10000) {
res <- Sl.postproc(b$Sl,fit,um$undrop,b$qrx$R,cov=TRUE,scale=scale)
object <- list(coefficients=res$beta,edf=res$edf,edf1=res$edf1,F=res$F,
object <- list(coefficients=res$beta,edf=res$edf,edf1=res$edf1,edf2=res$edf2,##F=res$F,
gcv.ubre=fit$reml,hat=res$hat,outer.info=list(iter=fit$iter,
message=fit$conv),optimizer="fast-REML",rank=ncol(um$X),
Ve=NULL,Vp=res$V,scale.estimated = scale<=0)
Ve=NULL,Vp=res$V,Vc=res$Vc,db.drho=fit$d1b,scale.estimated = scale<=0)
if (scale<=0) { ## get sp's and scale estimate
nsp <- length(fit$rho)
object$sig2 <- object$scale <- exp(fit$rho[nsp])
......@@ -1504,7 +1499,7 @@ bam.update <- function(b,data,chunk.size=10000) {
b$coefficients <- fit$b
b$edf <- post$edf
b$edf1 <- post$edf1
b$F <- post$F
##b$F <- post$F
b$full.sp <- fit$sp.full
b$gcv.ubre <- fit$score
b$hat <- post$hat
......@@ -1520,7 +1515,7 @@ bam.update <- function(b,data,chunk.size=10000) {
b$coefficients <- object$coefficients
b$edf <- object$edf
b$edf1 <- object$edf1
b$F <- object$F
##b$F <- object$F
b$full.sp <- object$sp.full
b$gcv.ubre <- object$gcv.ubre
b$hat <- object$hat
......
## (c) Simon N. Wood (2013, 2014) coxph model extended family.
## Released under GPL2 ...
cox.ph <- function (link = "identity") {
## Extended family object for Cox PH.
linktemp <- substitute(link)
if (!is.character(linktemp)) linktemp <- deparse(linktemp)
if (linktemp %in% c("identity")) stats <- make.link(linktemp)
else if (is.character(link)) {
stats <- make.link(link)
linktemp <- link
} else {
if (inherits(link, "link-glm")) {
stats <- link
if (!is.null(stats$name))
linktemp <- stats$name
}
else stop(linktemp, " link not available for coxph family; available link is \"identity\" ")
}
env <- new.env(parent = .GlobalEnv)
validmu <- function(mu) all(is.finite(mu))
# aic <- function(y, mu, theta=NULL, wt, dev) {
# ## this needs to call coxlpl - not really enough info
# ## use store and retrieve approach, actually this is not needed (gam.fit5 computes from likelihood)
# get(".log.partial.likelihood")
# }
## initialization is tough here... need data frame in reverse time order,
## and intercept removed from X...
preinitialize <- expression({
## code to evaluate in estimate.gam...
## sort y (time) into decending order, and
## re-order weights and rows of X accordingly
G$family$data <- list()
y.order <- order(G$y,decreasing=TRUE)
G$family$data$y.order <- y.order
G$y <- G$y[y.order]
G$X <- G$X[y.order,,drop=FALSE]
G$w <- G$w[y.order]
})
postproc <- expression({
## code to evaluate in estimate.gam, to do with data ordering and
## baseline hazard estimation...
## first get the estimated hazard and prediction information...
object$family$data <- G$family$hazard(G$y,G$X,object$coefficients,G$w)
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
## now put the survivor function in object$fitted
object$fitted.values <- exp(-object$family$data$h[object$family$data$r]*exp(object$linear.predictors))
## compute the null deviance...
s.base <- exp(-object$family$data$h[object$family$data$r]) ## baseline survival
s.base[s.base >= 1] <- 1 - 2*.Machine$double.eps ## avoid NA later
object$null.deviance <- ## sum of squares of null deviance residuals
2*sum(abs((object$prior.weights + log(s0.base) + object$prior.weights*(log(-log(s0.base))))))
## and undo the re-ordering...
object$linear.predictors[y.order] <- object$linear.predictors
object$fitted.values[y.order] <- object$fitted.values
object$y[y.order] <- object$y
object$prior.weights[y.order] <- object$prior.weights
})
initialize <- expression({
n <- rep(1, nobs)
if (is.null(start)) start <- rep(0,ncol(x))
})
hazard <- function(y, X,beta,wt) {
## 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)
tr <- unique(y);r <- match(y,tr);nt <- length(tr)
oo <- .C("coxpp",as.double(X%*%beta),A=as.double(X),as.integer(r),d=as.integer(wt),
h=as.double(rep(0,nt)),q=as.double(rep(0,nt)),km=as.double(rep(0,nt)),n=as.integer(nrow(X)),p=as.integer(ncol(X)),
nt=as.integer(nt),PACKAGE="mgcv")
p <- ncol(X)
list(tr=tr,h=oo$h,q=oo$q,a=matrix(oo$A[p*nt],p,nt),nt=nt,r=r,km=oo$km)
}
residuals <- function(object,type=c("deviance","martingale")) {
type <- match.arg(type)
w <- object$prior.weights;log.s <- log(object$fitted.values)
res <- w + log.s ## martingale residuals
if (type=="deviance") {
log.s[log.s>-1e-50] <- -1e-50
res <- sign(res)*sqrt(-2*(res + w * log(-log.s)))
}
res
}
predict <- function(family,se=FALSE,eta=NULL,y=NULL,
X=NULL,beta=NULL,off=NULL,Vb=NULL) {
## prediction function.
ii <- order(y,decreasing=TRUE) ## C code expects non-increasing
n <- nrow(X)
oo <- .C("coxpred",as.double(X[ii,]),t=as.double(y[ii]),as.double(beta),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),
s=as.double(rep(0,n)),se=as.double(rep(0,n)),PACKAGE="mgcv")
s <- sef <- oo$s
s[ii] <- oo$s
sef[ii] <- oo$se
if (se) return(list(fit=s,se.fit=sef)) else return(list(fit=s))
}
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) {
## function defining the cox model log lik.
## Calls C code "coxlpl"
## deriv codes: 0 - eval; 1 - grad and Hessian
## 2 - d1H (diagonal only)
## 3 - d1H; 4 d2H (diag)
## Hp is the preconditioned penalized hessian of the log lik
## which is of rank 'rank'.
## fh is a factorization of Hp - either its eigen decomposition
## 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)
##tr <- sort(unique(y),decreasing=TRUE)
tr <- unique(y)
r <- match(y,tr)
p <- ncol(X)
deriv <- deriv - 1
mu <- X%*%coef
g <- rep(0,p);H <- rep(0,p*p)
if (deriv > 0) {
M <- ncol(d1b)
d1H <- if (deriv==1) rep(0,p*M) else rep(0,p*p*M)
} else M <- d1H <- 0
if (deriv > 2) {
d2H <- rep(0,p*M*(M+1)/2)
#X <- t(forwardsolve(t(L),t(X)))
#d1b <- L %*% d1b; d2b <- L %*% d2b
if (is.list(fh)) {
ev <- fh
} else { ## need to compute eigen here
ev <- eigen(Hp,symmetric=TRUE)
if (rank < p) ev$values[(rank+1):p] <- 0
}
X <- X%*%(ev$vectors*D)
d1b <- t(ev$vectors)%*%(d1b/D); d2b <- t(ev$vectors)%*%(d2b/D)
} else trHid2H <- d2H <- 0
## note that the following call can not use .C(C_coxlpl,...) since the ll
## function is not in the mgcv namespace.
oo <- .C("coxlpl",as.double(mu),as.double(X),as.integer(r),as.integer(wt),
as.double(tr),n=as.integer(length(y)),p=as.integer(p),nt=as.integer(length(tr)),
lp=as.double(0),g=as.double(g),H=as.double(H),
d1b=as.double(d1b),d1H=as.double(d1H),d2b=as.double(d2b),d2H=as.double(d2H),
n.sp=as.integer(M),deriv=as.integer(deriv),PACKAGE="mgcv");
if (deriv==1) d1H <- matrix(oo$d1H,p,M) else
if (deriv>1) {
ind <- 1:(p^2)
d1H <- list()
for (i in 1:M) {
d1H[[i]] <- matrix(oo$d1H[ind],p,p)
ind <- ind + p^2
}
}
if (deriv>2) {
d2H <- matrix(oo$d2H,p,M*(M+1)/2)
#trHid2H <- colSums(d2H)
d <- ev$values
d[d>0] <- 1/d[d>0];d[d<=0] <- 0
trHid2H <- colSums(d2H*d)
}
assign(".log.partial.likelihood", oo$lp, envir=environment(sys.function()))
list(l=oo$lp,lb=oo$g,lbb=matrix(oo$H,p,p),d1H=d1H,d2H=d2H,trHid2H=trHid2H)
}
# environment(dev.resids) <- environment(aic) <- environment(getTheta) <-
# environment(rd)<- environment(qf)<- environment(variance) <- environment(putTheta)
#environment(aic) <-
environment(ll) <- env
structure(list(family = "Cox PH", link = linktemp, linkfun = stats$linkfun,
linkinv = stats$linkinv, ll=ll, ## aic = aic,
mu.eta = stats$mu.eta,
initialize = initialize,preinitialize=preinitialize,postproc=postproc,
hazard=hazard,predict=predict,residuals=residuals,
validmu = validmu, valideta = stats$valideta,
rd=rd,qf=qf,drop.intercept = TRUE,
ls=1, ## signal ls not needed
available.derivs = 2 ## can use full Newton here
),
class = c("general.family","extended.family","family"))
} ## cox.ph
This diff is collapsed.
## code for fast REML computation. key feature is that first and
## second derivatives come at no increase in leading order
## computational cost, relative to evaluation!
## (c) Simon N. Wood, 2010-2012
## (c) Simon N. Wood, 2010-2014
Sl.setup <- function(G) {
## Sets up a list representing a block diagonal penalty matrix.
......@@ -58,20 +58,11 @@ Sl.setup <- function(G) {
} ## finished paraPen
## now work through the smooths....
for (i in 1:length(G$smooth)) {
if (length(G$smooth)) for (i in 1:length(G$smooth)) {
if (!is.null(G$smooth[[i]]$fixed)&&G$smooth[[i]]$fixed) m <- 0 else
m <- length(G$smooth[[i]]$S)
# if (m>0&&!is.null(G$smooth[[i]]$fx)) { ## drop any fixed
# for (j in m:1) {
# if (G$smooth[[i]]$fx[j]==TRUE) {
# G$smooth[[i]]$S[[j]] <- NULL
# m <- m - 1
# }
# }
# }
if (m>0) {
Sl[[b]] <- list()
Sl[[b]]$start <- G$smooth[[i]]$first.para
......@@ -116,7 +107,6 @@ Sl.setup <- function(G) {
b <- b + 1
}
} else { ## not possible to split
#Sl[[b]]$rank <- nrow(G$smooth[[i]]$S[[1]]) - G$smooth[[i]]$null.space.dim+1
Sl[[b]]$S <- G$smooth[[i]]$S
} ## additive block finished
......@@ -129,6 +119,8 @@ Sl.setup <- function(G) {
## Singletons need to be transformed to identity penalties, while
## multiples need to be projected into total penalty range space.
if (length(Sl)==0) return(Sl) ## nothing to do
np <- ncol(G$X)
E <- matrix(0,np,np) ## well scaled square root penalty
lambda <- rep(0,0)
......@@ -155,7 +147,6 @@ Sl.setup <- function(G) {
Sl[[b]]$D <- t(D*t(U)) ## D <- U%*%diag(D)
Sl[[b]]$Di <- t(U)/D
## so if X is smooth model matrix X%*%D is re-parameterized form
## Sl[[b]]$D <- D;
Sl[[b]]$ind <- ind
}
## add penalty square root into E
......@@ -179,8 +170,7 @@ Sl.setup <- function(G) {
}
Sl[[b]]$ind <- rep(FALSE,ncol(U))