##' GOF for Cox PH regression
##'
##' Cumulative score process residuals for Cox PH regression
##' p-values based on Lin, Wei, Ying resampling.
##' @param object is phreg object
##' @param n.sim number of simulations for score processes
##' @param silent to show timing estimate will be produced for longer jobs
##' @param robust to control wether robust dM_i(t) or dN_i are used for simulations
##' @param ... Additional arguments to lower level funtions
##' @author Thomas Scheike and Klaus K. Holst
##' @export
##' @aliases gof.phreg
##' @examples
##' library(mets)
##' data(sTRACE)
##'
##' m1 <- phreg(Surv(time,status==9)~vf+chf+diabetes,data=sTRACE)
##' gg <- gof(m1)
##' gg
##' par(mfrow=c(1,3))
##' plot(gg)
##'
##' m1 <- phreg(Surv(time,status==9)~strata(vf)+chf+diabetes,data=sTRACE)
##' ## to get Martingale ~ dN based simulations
##' gg <- gof(m1)
##' gg
##'
##' ## to get Martingale robust simulations, specify cluster in call
##' sTRACE$id <- 1:500
##' m1 <- phreg(Surv(time,status==9)~vf+chf+diabetes+cluster(id),data=sTRACE)
##' gg <- gof(m1)
##' gg
##'
##' m1 <- phreg(Surv(time,status==9)~strata(vf)+chf+diabetes+cluster(id),data=sTRACE)
##' gg <- gof(m1)
##' gg
##' @export
gof.phreg <- function(object,n.sim=1000,silent=1,robust=NULL,...)
{# {{{
### test for proportionality
p <- length(object$coef)
nnames <- names(object$coef)
ii <- solve(object$hessian)
jumptimes <- object$jumptimes
Pt <- object$hessianttime
U <- object$U
Pt <- apply(Pt,2,cumsum)
Ut <- apply(U,2,cumsum)
nd <- nrow(object$U)
Pt <- .Call("XXMatFULL",Pt,p,PACKAGE="mets")$XXf
Pt <- .Call("CubeMat",Pt,ii,PACKAGE="mets")$XXX
sup <- matrix(0,n.sim,nrow(ii))
hatti <- matrix(0,nd,nrow(ii))
obs <- apply(abs(Ut),2,max)
if (is.null(robust)) {
if (!is.null(object$call.id)) robust <- TRUE else robust <- FALSE
if (inherits(object,c("cifreg"))) {
robust <- FALSE ## approximative test using MG framework
### if (is.null(object$cox.prep)) stop("Must be called with cox.prep=TRUE\n")
}
}
### cluster call or robust \hat M_i(t) based
if (robust) {# {{{
xx <- object$cox.prep
Z <- xx$X
rrw <- c(xx$sign * exp(Z %*% coef(object) + xx$offset)*xx$weights)
UdN <- xx$weights[xx$jumps+1]*object$U
### also weights
nn <- nrow(Z)
tt <- system.time(simcox1<- .Call("PropTestCoxClust",UdN,Pt,rrw,xx$X,
object$S0,object$E,
10,obs,nn,xx$id,xx$strata,xx$nstrata,object$strata.jumps,xx$jumps))
prt <- n.sim*tt[3]/(10*60)
if (prt>1 & silent==0) cat(paste("Predicted time minutes",signif(prt,2),"\n"))
simcox <- .Call("PropTestCoxClust",UdN,Pt,rrw,xx$X,
object$S0,object$E,
n.sim,obs,nn,xx$id,xx$strata,xx$nstrata,object$strata.jumps,xx$jumps)
} else {# }}}
### or dN_i based # {{{
tt <- system.time(simcox1<-.Call("PropTestCox",U,Pt,10,obs,PACKAGE="mets"))
prt <- n.sim*tt[3]/(10*60)
if (prt>1 & silent==0) cat(paste("Predicted time minutes",signif(prt,2),"\n"))
simcox <- .Call("PropTestCox",U,Pt,n.sim,obs,PACKAGE="mets")
}# }}}
sup <- simcox$supUsim
res <- cbind(obs,simcox$pval)
colnames(res) <- c("Sup|U(t)|","pval")
rownames(res) <- nnames
if (silent==0) {
cat("Cumulative score process test for Proportionality:\n")
prmatrix(round(res,digits=2))
}
out <- list(jumptimes=object$jumptimes,supUsim=sup,res=res,supU=obs,
pvals=simcox$pval,score=Ut,simUt=simcox$simUt,type="prop",robust=robust)
class(out) <- "gof.phreg"
return(out)
}# }}}
##' GOF for Cox covariates in PH regression
##'
##' Cumulative residuals after model matrix for Cox PH regression
##' p-values based on Lin, Wei, Ying resampling.
##'
##' That is, computes
##' \deqn{
##' U(t) = \int_0^t M^t d \hat M
##' }
##' and resamples its asymptotic distribution.
##'
##' This will show if the residuals are consistent with the model. Typically,
##' M will be a design matrix for the continous covariates that gives for example
##' the quartiles, and then the plot will show if for the different quartiles of the covariate the risk
##' prediction is consistent over time (time x covariate interaction).
##'
##' @param formula formula for cox regression
##' @param data data for model
##' @param offset offset
##' @param weights weights
##' @param modelmatrix matrix for cumulating residuals
##' @param n.sim number of simulations for score processes
##' @param silent to keep it absolutely silent, otherwise timing estimate will be prduced for longer jobs.
##' @param ... Additional arguments to lower level funtions
##' @author Thomas Scheike and Klaus K. Holst
##' @export
##' @examples
##' library(mets)
##' data(TRACE)
##' set.seed(1)
##' TRACEsam <- blocksample(TRACE,idvar="id",replace=FALSE,100)
##'
##' dcut(TRACEsam) <- ~.
##' mm <- model.matrix(~-1+factor(wmicat.4),data=TRACEsam)
##' m1 <- gofM.phreg(Surv(time,status==9)~vf+chf+wmi,data=TRACEsam,modelmatrix=mm)
##' summary(m1)
##' if (interactive()) {
##' par(mfrow=c(2,2))
##' plot(m1)
##' }
##'
##' m1 <- gofM.phreg(Surv(time,status==9)~strata(vf)+chf+wmi,data=TRACEsam,modelmatrix=mm)
##' summary(m1)
##'
##' ## cumulative sums in covariates, via design matrix mm
##' mm <- cumContr(TRACEsam$wmi,breaks=10,equi=TRUE)
##' m1 <- gofM.phreg(Surv(time,status==9)~strata(vf)+chf+wmi,data=TRACEsam,
##' modelmatrix=mm,silent=0)
##' summary(m1)
##'
##' @export
gofM.phreg <- function(formula,data,offset=NULL,weights=NULL,modelmatrix=NULL,
n.sim=1000,silent=1,...)
{# {{{
if (is.null(modelmatrix)) stop(" must give matrix for cumulating residuals\n");
cox1 <- phreg(formula,data,offset=NULL,weights=NULL,Z=modelmatrix,cumhaz=FALSE,...)
offsets <- cox1$X %*% cox1$coef
if (!is.null(offset)) offsets <- offsets*offset
if (!is.null(cox1$strata.name))
coxM <- phreg(cox1$model.frame[,1]~modelmatrix+strata(cox1$strata),data,offset=offsets,weights=weights,no.opt=TRUE,cumhaz=FALSE,no.var=1,...)
else coxM <- phreg(cox1$model.frame[,1]~modelmatrix,data,offset=offsets,weights=weights,no.opt=TRUE,cumhaz=FALSE,no.var=1,...)
nnames <- colnames(modelmatrix)
Ut <- apply(coxM$U,2,cumsum)
jumptimes <- coxM$jumptimes
U <- coxM$U
Ubeta <- cox1$U
ii <- -solve(cox1$hessian)
EE <- .Call("vecMatMat",coxM$E,cox1$E,PACKAGE="mets")$vXZ;
Pt <- cox1$ZX - EE
Pt <- apply(Pt,2,cumsum)
betaiid <- t(ii %*% t(Ubeta))
obs <- apply(abs(Ut),2,max)
simcox <- .Call("ModelMatrixTestCox",U,Pt,betaiid,n.sim,obs,PACKAGE="mets")
sup <- simcox$supUsim
res <- cbind(obs,simcox$pval)
colnames(res) <- c("Sup_t |U(t)|","pval")
rownames(res) <- nnames
if (silent==0) {
cat("Cumulative score process test for modelmatrix:\n")
prmatrix(round(res,digits=2))
}
## pvals efter z i model.matrix sup_z | M(z,tau) |
Utlast <- max(abs(tail(Ut,1)))
maxlast <- apply(abs(simcox$last),1,max)
pval.last <- mean(maxlast>=Utlast)
res.last <- matrix(c(Utlast,pval.last),1,2)
colnames(res.last) <- c("Sup_z |U(tau,z)|","pval")
rownames(res.last) <- "matrixZ"
out <- list(jumptimes=jumptimes,supUsim=simcox$supUsim,res=res,supU=obs,
pvals=simcox$pval,score=Ut,simUt=simcox$simUt,
simUtlast=simcox$last,Utlast=Utlast,pval.last=pval.last,
res.last=res.last, type="modelmatrix")
class(out) <- "gof.phreg"
return(out)
}# }}}
##' GOF for Cox covariates in PH regression
##'
##' That is, computes
##' \deqn{
##' U(z,\tau) = \int_0^\tau M(z)^t d \hat M
##' }
##' and resamples its asymptotic distribution.
##'
##' This will show if the residuals are consistent with the model evaulated in the z covariate.
##' M is here chosen based on a grid (z_1, ..., z_m) and the different columns are \eqn{I(Z_i \leq z_l)}.
##' for \eqn{l=1,...,m}.
##' The process in z is resampled to find extreme values. The time-points of evuluation is by default
##' 50 points, chosen as 2%,4%,..., percentiles of the covariates.
##'
##' The p-value is valid but depends on the chosen grid. When the number of break points are high
##' this will give the orginal test of Lin, Wei and Ying for linearity, that is also computed in
##' the timereg package.
##'
##' @param formula formula for cox regression
##' @param data data for model
##' @param vars which variables to test for linearity
##' @param offset offset
##' @param weights weights
##' @param breaks number of breaks for cumulatives in covarirate direction
##' @param equi equidistant breaks or not
##' @param n.sim number of simulations for score processes
##' @param silent to keep it absolutely silent, otherwise timing estimate will be prduced for longer jobs.
##' @param ... Additional arguments to lower level funtions
##' @author Thomas Scheike and Klaus K. Holst
##' @export
##' @examples
##' library(mets)
##' data(TRACE)
##' set.seed(1)
##' TRACEsam <- blocksample(TRACE,idvar="id",replace=FALSE,100)
##'
##' ## cumulative sums in covariates, via design matrix mm
##' \donttest{ ## Reduce Ex.Timings
##' m1 <- gofZ.phreg(Surv(time,status==9)~strata(vf)+chf+wmi+age,data=TRACEsam)
##' summary(m1)
##' plot(m1,type="z")
##' }
##' @aliases cumContr
##' @export
gofZ.phreg <- function(formula,data,vars=NULL,offset=NULL,weights=NULL,breaks=50,equi=FALSE,
n.sim=1000,silent=1,...)
{# {{{
if (is.null(vars)) {
vars <- NULL
## find strata var's
yxzf <- procform(formula,x=NULL,z=NULL,data=data,do.filter=FALSE)
avars <- all.vars(formula[-2])
svar <- grep("strata",yxzf$predictor)
if (length(svar)>=1) {
avars <- avars[-svar]
}
## check that it is not a factor and that there are more than 2 levels
for (vv in avars) {
if (length(unique(data[,vv]))>2 & !is.factor(data[,vv]))
vars <- c(vars,vv)
}
}
gres <- list()
res <- matrix(0,length(vars),2)
colnames(res) <- c("Sup_z |U(tau,z)|","pval")
rownames(res) <- vars
i <- 1
for (vv in vars) {
modelmatrix <- cumContr(data[,vv],breaks=breaks,equi=equi)
lres <- gofM.phreg(formula,data,modelmatrix=modelmatrix)
lres$xaxs <- attr(modelmatrix,"breaks")
res[i,] <- c(lres$Utlast,lres$pval.last)
i <- i+1
lres <- list(lres)
names(lres) <- vv
gres <- c(gres,lres)
}
out <- list(res=res,Zres=gres,type="Zmodelmatrix")
class(out) <- c("gof.phreg")
return(out)
}# }}}
##' @export
cumContr <- function(data,breaks=4,probs=NULL,equi=TRUE,na.rm=TRUE,unique.breaks=TRUE,...)
{# {{{
if (is.vector(data)) {
if (is.list(breaks))
breaks <- unlist(breaks)
if (length(breaks) == 1) {
if (!is.null(probs)) {
breaks <- quantile(data, probs, na.rm = na.rm,...)
breaks <- breaks[-1]
}
else {
if (!equi) {
probs <- seq(0, 1, length.out = breaks + 1)
breaks <- quantile(data, probs, na.rm = na.rm, ...)
if (unique.breaks) breaks <- unique(breaks)
breaks <- breaks[-1]
}
if (equi) {
rr <- range(data, na.rm = na.rm)
breaks <- seq(rr[1], rr[2], length.out = breaks + 1)
breaks <- breaks[-1]
}
}
}
if (sum(duplicated(breaks)) == 0) {
i <- 0;
gm <- matrix(0,length(data),length(breaks))
for (bb in breaks) {
i <- i+1
gm[,i] <- (data <= bb)*1
}
} else {
wd <- which(duplicated(breaks))
mb <- min(diff(breaks[-wd]))
breaks[wd] <- breaks[wd] + (mb/2) * seq(length(wd))/length(wd)
i <- 0; gm <- matrix(0,length(data),length(breaks))
for (bb in breaks) {
i <- i+1
gm[,i] <- (data <= bb)*1
}
warning(paste("breaks duplicated"))
}
colnames(gm) <- paste("<=",breaks,sep="")
attr(gm,"breaks") <- breaks
return(gm)
}
}# }}}
##' Stratified baseline graphical GOF test for Cox covariates in PH regression
##'
##' Looks at stratified baseline in Cox model and plots all baselines versus each
##' other to see if lines are straight, with 50 resample versions under the
##' assumptiosn that the stratified Cox is correct
##'
##' @param x phreg object
##' @param sim to simulate som variation from cox model to put on graph
##' @param silent to keep it absolutely silent
##' @param lm addd line to plot, regressing the cumulatives on each other
##' @param ... Additional arguments to lower level funtions
##' @author Thomas Scheike and Klaus K. Holst
##' @export
##' @examples
##' data(tTRACE)
##'
##' m1 <- phreg(Surv(time,status==9)~strata(vf)+chf+wmi,data=tTRACE)
##' m2 <- phreg(Surv(time,status==9)~vf+strata(chf)+wmi,data=tTRACE)
##' par(mfrow=c(2,2))
##'
##' gofG.phreg(m1)
##' gofG.phreg(m2)
##'
##' bplot(m1,log="y")
##' bplot(m2,log="y")
##' @export
gofG.phreg <- function(x,sim=0,silent=1,lm=TRUE,...)
{# {{{
p <- length(x$coef)
nnames <- names(x$coef)
strata <- x$strata[x$jumps]
nstrata <- x$nstrata
jumptimes <- x$jumptimes
cumhaz <- x$cumhaz
ms <- match(x$strata.name,names(x$model.frame))
lstrata <- levels(x$model.frame[,ms])
stratn <- substring(x$strata.name,8,nchar(x$strata.name)-1)
stratnames <- paste(stratn,lstrata,sep=":")
if (is.null(cumhaz)) stop("Must run phreg with cumhaz=TRUE (default)");
if (nstrata==1) stop("Stratified Cox to look at baselines");
if ((x$no.opt) | is.null(x$coef)) fixbeta<- 1 else fixbeta <- 0
for (i in 0:(nstrata-2))
for (j in (i+1):(nstrata-1)) {
iij <- which(strata %in% c(i,j))
ii <- which(strata %in% i)
ij <- which(strata %in% j)
dijjumps <- jumptimes[iij]
cumhazi <- cpred(cumhaz[ii,],dijjumps)
cumhazj <- cpred(cumhaz[ij,],dijjumps)
plot(cumhazj[,2],cumhazi[,2],type="s",lwd=2,xlab=stratnames[j+1],ylab=stratnames[i+1])
graphics::title(paste("Stratified baselines for",stratn))
if ((fixbeta==0 | sim==0) & lm )
graphics::legend("topleft",c("Nonparametric","lm"),lty=1,col=1:2)
ab <- lm(cumhazi[,2]~-1+cumhazj[,2])
if (sim==1 & fixbeta==0) {
Pt <- DLambeta.t <- apply(x$E/c(x$S0),2,cumsumstrata,strata,nstrata)
II <- -solve(x$hessian)
betaiid <- t(II %*% t(x$U))
simband <- .Call("simBandCumHazCox",1/x$S0,Pt,betaiid,50,rep(1,nrow(Pt)),PACKAGE="mets")
simU <-simband$simUt
for (k in 1:50)
{
di <- cpred(cbind(jumptimes[ii],simU[ii,k]),dijjumps)[,2]
dj <- cpred(cbind(jumptimes[ij],simU[ij,k]),dijjumps)[,2]
lines(cumhazj[,2]+dj,cumhazi[,2]+di,type="s",lwd=0.1,col=3)
}
}
lines(cumhazj[,2],cumhazi[,2],type="s",lwd=2,col=1)
if (lm==TRUE) abline(c(0,coef(ab)),col=2,lwd=2)
}
}# }}}
##' @export
plot.gof.phreg <- function(x,col=3,type=NULL,...)
{# {{{
if (is.null(type)) {
if (x$type=="prop") type <- "time"
if (x$type=="modelmatrix" ) type <- "modelmatrix"
if (x$type=="Zmodelmatrix") type <- "z"
}
if (type=="time" || type=="modelmatrix") {
p <- ncol(x$score)
for (i in 1:p)
{
simU <- x$simUt[,(0:49)*p+i]
rsU <- max(abs(simU))
rsU <- max(rsU,abs(x$score[,i]))
plot(x$jumptimes,x$score[,i],type="s",ylim=c(-rsU,rsU),xlab="",ylab="")
title(main=rownames(x$res)[i])
matlines(x$jumptimes,simU,type="s",lwd=0.3,col=col)
lines(x$jumptimes,x$score[,i],type="s",lwd=1.5)
}
} else {
if (type=="modelmatrix") {
obsz <- c(tail(x$score,1))
times <- 1:length(obsz)
rsU <- max(max(abs(obsz)),max(abs(x$simUtlast[1:50,])))
plot(times,obsz,type="l",ylim=c(-rsU,rsU),xlab="",ylab="")
matlines(times,t(x$simUtlast[1:50,]),type="l",lwd=0.3,col=col)
## redraw with thick to make observed clear
lines(times,obsz,lwd=2,col=1)
} else {
for (i in 1:length(x$Zres))
{
xr <- x$Zres[[i]]
obsz <- c(tail(xr$score,1))
times <- xr$xaxs
rsU <- max(max(abs(obsz)),max(abs(xr$simUtlast[1:50,])))
plot(times,obsz,type="l",ylim=c(-rsU,rsU),xlab="",ylab="")
title(main=rownames(x$res)[i])
matlines(times,t(xr$simUtlast[1:50,]),type="l",lwd=0.3,col=col)
## redraw with thick to make observed clear
lines(times,obsz,lwd=2,col=1)
}
}
}
}# }}}
##' @export
summary.gof.phreg <- function(object,...)
{# {{{
if (object$type=="prop")
cat("Cumulative score process test for Proportionality:\n")
else cat("Cumulative residuals versus modelmatrix :\n")
print(object$res)
if (!is.null(object$res.last)) {
cat("\n")
cat("Cumulative score process versus covariates (discrete z via model.matrix):\n")
print(object$res.last)
}
} # }}}
##' @export
print.gof.phreg <- function(x,...)
{# {{{
if (x$type=="prop")
cat("Cumulative score process test for Proportionality:\n")
else cat("Cumulative residuals versus modelmatrix :\n")
print(x$res)
if (!is.null(x$res.last)) {
cat("\n")
cat("Cumulative score process versus covariates (discrete z via model.matrix):\n")
print(x$res.last)
}
} # }}}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.