Nothing
###{{{ phreg0
phreg0 <- function(X,entry,exit,status,id=NULL,strata=NULL,beta,stderr=TRUE,method="NR",...) {
p <- ncol(X)
if (missing(beta)) beta <- rep(0,p)
if (p==0) X <- cbind(rep(0,length(exit)))
if (!is.null(strata)) {
stratalev <- levels(strata)
strataidx <- lapply(stratalev,function(x) which(strata==x))
if (!all(unlist(lapply(strataidx,function(x) length(x)>0))))
stop("Strata without any observation")
dd <- lapply(strataidx, function(ii)
.Call("FastCoxPrep",
entry[ii],exit[ii],status[ii],
as.matrix(X)[ii,,drop=FALSE],
id[ii],
!is.null(id),
!is.null(entry),
package="mets"))
if (!is.null(id))
id <- unlist(lapply(dd,function(x) x$id[x$jumps+1]))
obj <- function(pp,U=FALSE,all=FALSE) {
val <- lapply(dd,function(d)
with(d,
.Call("FastCoxPL",pp,X,XX,sign,jumps,package="mets")))
ploglik <- Reduce("+",lapply(val,function(x) x$ploglik))
gradient <- Reduce("+",lapply(val,function(x) x$gradient))
hessian <- Reduce("+",lapply(val,function(x) x$hessian))
if (all) {
U <- do.call("rbind",lapply(val,function(x) x$U))
time <- lapply(dd,function(x) x$time[x$ord+1])
ord <- lapply(dd,function(x) x$ord+1)
jumps <- lapply(dd,function(x) x$jumps+1)
jumptimes <- lapply(dd,function(x) x$time[x$ord+1][x$jumps+1])
S0 <- lapply(val,function(x) x$S0)
nevent <- unlist(lapply(S0,length))
return(list(ploglik=ploglik,gradient=gradient,hessian=hessian,
U=U,S0=S0,nevent=nevent,
ord=ord,time=time,jumps=jumps,jumptimes=jumptimes))
}
structure(-ploglik,gradient=-gradient,hessian=-hessian)
}
} else {
browser()
system.time(dd <- .Call("mets_FastCoxPrep",
entry,exit,status,X,
as.integer(seq_along(entry)),
is.null(id),
!is.null(entry),
package="mets"))
if (!is.null(id))
id <- dd$id[dd$jumps+1]
obj <- function(pp,U=FALSE,all=FALSE) {
val <- with(dd,
.Call("FastCoxPL",pp,X,XX,sign,jumps,package="mets"))
if (all) {
val$time <- dd$time[dd$ord+1]
val$ord <- dd$ord+1
val$jumps <- dd$jumps+1
val$jumptimes <- val$time[val$jumps]
val$nevent <- length(val$S0)
return(val)
}
with(val, structure(-ploglik,gradient=-gradient,hessian=-hessian))
}
}
opt <- NULL
if (p>0) {
if (tolower(method)=="nr") {
opt <- lava::NR(beta,obj,...)
opt$estimate <- opt$par
} else {
opt <- nlm(obj,beta,...)
opt$method <- "nlm"
}
cc <- opt$estimate; names(cc) <- colnames(X)
if (!stderr) return(cc)
val <- c(list(coef=cc),obj(opt$estimate,all=TRUE))
} else {
val <- obj(0,all=TRUE)
val[c("ploglik","gradient","hessian","U")] <- NULL
}
res <- c(val,
list(strata=strata,
entry=entry,
exit=exit,
status=status,
p=p,
X=X,
id=id, opt=opt))
class(res) <- "phreg"
res
}
###}}} phreg0
###{{{ simcox
##' @export
simCox <- function(n=1000, seed=1, beta=c(1,1), entry=TRUE) {
if (!is.null(seed))
set.seed(seed)
m <- lvm()
regression(m,T~X1+X2) <- beta
distribution(m,~T+C) <- coxWeibull.lvm(scale=1/100)
distribution(m,~entry) <- coxWeibull.lvm(scale=1/10)
m <- eventTime(m,time~min(T,C=0),"status")
d <- sim(m,n);
if (!entry) d$entry <- 0
else d <- subset(d, time>entry,select=-c(T,C))
return(d)
}
###}}} simcox
###{{{ phreg
##' Fast Cox PH regression
##'
##' Fast Cox PH regression
##' @param formula formula with 'Surv' outcome (see \code{coxph})
##' @param data data frame
##' @param ... Additional arguments to lower level funtions
##' @author Klaus K. Holst
##' @export
##' @aliases phreg phreg.par
##' @examples
##' simcox <- function(n=1000, seed=1, beta=c(1,1), entry=TRUE) {
##' if (!is.null(seed))
##' set.seed(seed)
##' library(lava)
##' m <- lvm()
##' regression(m,T~X1+X2) <- beta
##' distribution(m,~T+C) <- coxWeibull.lvm(scale=1/100)
##' distribution(m,~entry) <- coxWeibull.lvm(scale=1/10)
##' m <- eventTime(m,time~min(T,C=0),"status")
##' d <- sim(m,n);
##' if (!entry) d$entry <- 0
##' else d <- subset(d, time>entry,select=-c(T,C))
##' return(d)
##' }
##' \dontrun{
##' n <- 1e3;
##' d <- mets:::simCox(n); d$id <- seq(nrow(d)); d$group <- factor(rbinom(nrow(d),1,0.5))
##'
##' (m1 <- phreg(Surv(entry,time,status)~X1+X2,data=d))
##' (m2 <- coxph(Surv(entry,time,status)~X1+X2+cluster(id),data=d))
##' (coef(m3 <-cox.aalen(Surv(entry,time,status)~prop(X1)+prop(X2),data=d)))
##'
##'
##' (m1b <- phreg(Surv(entry,time,status)~X1+X2+strata(group),data=d))
##' (m2b <- coxph(Surv(entry,time,status)~X1+X2+cluster(id)+strata(group),data=d))
##' (coef(m3b <-cox.aalen(Surv(entry,time,status)~-1+group+prop(X1)+prop(X2),data=d)))
##'
##' m <- phreg(Surv(entry,time,status)~X1*X2+strata(group)+cluster(id),data=d)
##' m
##' plot(m,ylim=c(0,1))
##' }
phreg <- function(formula,data,...) {
cl <- match.call()
m <- match.call(expand.dots = TRUE)[1:3]
special <- c("strata", "cluster")
Terms <- terms(formula, special, data = data)
m$formula <- Terms
m[[1]] <- as.name("model.frame")
m <- eval(m, parent.frame())
Y <- model.extract(m, "response")
if (!is.Surv(Y)) stop("Expected a 'Surv'-object")
if (ncol(Y)==2) {
exit <- Y[,1]
entry <- NULL ## rep(0,nrow(Y))
status <- Y[,2]
} else {
entry <- Y[,1]
exit <- Y[,2]
status <- Y[,3]
}
id <- strata <- NULL
if (!is.null(attributes(Terms)$specials$cluster)) {
ts <- survival::untangle.specials(Terms, "cluster")
Terms <- Terms[-ts$terms]
id <- m[[ts$vars]]
}
if (!is.null(stratapos <- attributes(Terms)$specials$strata)) {
ts <- survival::untangle.specials(Terms, "strata")
Terms <- Terms[-ts$terms]
strata <- m[[ts$vars]]
}
X <- model.matrix(Terms, m)
if (!is.null(intpos <- attributes(Terms)$intercept))
X <- X[,-intpos,drop=FALSE]
if (ncol(X)==0) X <- matrix(nrow=0,ncol=0)
res <- c(phreg0(X,entry,exit,status,id,strata,...),list(call=cl,model.frame=m))
class(res) <- "phreg"
res
}
###}}} phreg
###{{{ vcov
##' @export
vcov.phreg <- function(object,...) {
res <- crossprod(ii <- iid(object,...))
attributes(res)$ncluster <- attributes(ii)$ncluster
attributes(res)$invhess <- attributes(ii)$invhess
colnames(res) <- rownames(res) <- names(coef(object))
res
}
###}}} vcov
###{{{ coef
##' @export
coef.phreg <- function(object,...) {
object$coef
}
###}}} coef
###{{{ iid
##' @export
iid.phreg <- function(x,...) {
invhess <- solve(x$hessian)
ncluster <- NULL
if (!is.null(x$id)) {
ii <- mets::cluster.index(x$id)
UU <- matrix(nrow=ii$uniqueclust,ncol=ncol(invhess))
for (i in seq(ii$uniqueclust)) {
UU[i,] <- colSums(x$U[ii$idclustmat[i,seq(ii$cluster.size[i])]+1,,drop=FALSE])
}
ncluster <- nrow(UU)
} else {
UU <- x$U
}
structure(UU%*%invhess,invhess=invhess,ncluster=ncluster)
}
###}}}
###{{{ summary
##' @export
summary.phreg <- function(object,se="robust",...) {
cc <- ncluster <- NULL
if (length(object$p)>0) {
I <- -solve(object$hessian)
V <- vcov(object)
cc <- cbind(coef(object),diag(V)^0.5,diag(I)^0.5)
cc <- cbind(cc,2*(pnorm(abs(cc[,1]/cc[,2]),lower.tail=FALSE)))
colnames(cc) <- c("Estimate","S.E.","dU^-1/2","P-value")
if (!is.null(ncluster <- attributes(V)$ncluster))
rownames(cc) <- names(coef(object))
}
Strata <- levels(object$strata)
if (!is.null(Strata)) {
n <- unlist(lapply(object$time,length))
} else {
n <- length(object$time)
}
res <- list(coef=cc,n=n,nevent=object$nevent,
strata=Strata,ncluster=ncluster)
class(res) <- "summary.phreg"
res
}
###}}} summary
###{{{ print.summary
##' @export
print.summary.phreg <- function(x,max.strata=5,...) {
cat("\n")
nn <- cbind(x$n, x$nevent)
rownames(nn) <- levels(x$strata); colnames(nn) <- c("n","events")
if (is.null(rownames(nn))) rownames(nn) <- rep("",NROW(nn))
if (length(x$strata)>max.strata) {
nn <- rbind(c(colSums(nn),length(x$strata)));
colnames(nn) <- c("n","events","stratas")
rownames(nn) <- ""
}
print(nn,quote=FALSE)
if (!is.null(x$ncluster)) cat("\n ", x$ncluster, " clusters\n",sep="")
if (!is.null(x$coef)) {
cat("\n")
printCoefmat(x$coef,...)
}
cat("\n")
}
###}}} print.summary
###{{{ predict
predictPhreg <- function(jumptimes,S0,beta,time=NULL,X=NULL,surv=FALSE,...) {
## Brewslow estimator
chaz <- cbind(jumptimes,cumsum(1/S0))
if (!is.null(time)) {
chaz <- Cpred(chaz,time)
}
colnames(chaz) <- c("time","chaz")
if (!is.null(X)) {
H <- exp(X%*%beta)
if (nrow(chaz)==length(H)) {
chaz[,2] <- chaz[,2]*H
} else {
chaz2 <- c()
X <- rbind(X)
for (i in seq(nrow(X)))
chaz2 <- rbind(chaz2,
cbind(chaz[,1],chaz[,2]*H[i],
rep(1,nrow(chaz))%x%X[i,,drop=FALSE]))
chaz <- chaz2;
nn <- c("time","chaz",names(beta))
colnames(chaz) <- nn
}
}
if (surv) {
chaz[,2] <- exp(-chaz[,2])
colnames(chaz)[2] <- "surv"
}
return(chaz)
}
##' @export
predict.phreg <- function(object,data,surv=FALSE,time=object$exit,X=object$X,strata=object$strata,...) {
if (!is.null(object$strata)) {
lev <- levels(object$strata)
if (!is.null(object$strata) &&
!(is.list(time) & !is.data.frame(time)) &&
!(is.list(X) & !is.data.frame(X))) {
X0 <- X
time0 <- time
X <- time <- c()
for (i in seq(length(lev))) {
idx <- which(strata==lev[i])
X <- c(X,list(X0[idx,,drop=FALSE]))
time <- c(time,list(time0[idx]))
}
}
chaz <- c()
for (i in seq(length(lev)))
chaz <- c(chaz,list(predictPhreg(object$jumptimes[[i]],
object$S0[[i]],
coef(object),
time[[i]],X[[i]],surv)))
names(chaz) <- lev
} else {
chaz <- predictPhreg(object$jumptimes,object$S0,coef(object),time,X,surv)
}
return(chaz)
}
###}}} predict
###{{{ plot
##' @export
plot.phreg <- function(x,surv=TRUE,X=NULL,time=NULL,add=FALSE,...) {
if (!is.null(X) && nrow(X)>1) {
P <- lapply(split(X,seq(nrow(X))),function(xx) predict(x,X=xx,time=time,surv=surv))
} else {
P <- predict(x,X=X,time=time,surv=surv)
}
if (!is.list(P)) {
if (add) {
lines(P,type="s",...)
} else {
plot(P,type="s",...)
}
return(invisible(P))
}
if (add) {
lines(P[[1]][,1:2],type="s",lty=1,col=1,...)
} else {
plot((P[[1]])[,1:2],type="s",lty=1,col=1,...)
}
for (i in seq_len(length(P)-1)+1) {
lines(P[[i]][,1:2],type="s",lty=i,col=i,...)
}
return(invisible(P))
}
###}}} plot
###{{{ print
##' @export
print.phreg <- function(x,...) {
cat("Call:\n")
dput(x$call)
print(summary(x),...)
}
###}}} print
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.