transform.hazards

Transforms cumulative hazard estimates to estimate survival analysis parameters specified by differential equations. We demonstrate how the method works om some commonly studied parameters.

Parameters

We are interested in assessing parameters that solve differential equations driven by cumulative hazards \(A\), i.e. \[\begin{equation}X_t = X_0 + \int_0^t F(X_s) dA_s,\end{equation}\] where \(F = (F_1,F_2,\cdots)\) is Lipschitz, two times continuously differentiable, and satisfies a linear growth bound. Several examples of such parameters can be found in ¹ ². The main function pluginEstimate in this package estimate such parameters. Among other things, pluginEstimate has the following input

The increments of the cumulative hazard estimates
The function \(F\)
The Jacobian matrices of the columns of \(F\), i.e. \(J_{F_1},J_{F_2},\cdots\), as a list

We illustrate how to provide the correct inputs by use of examples below.

Survival

The survival function \(S\) solves a differential equation with initial value 1; \[\begin{equation} S_t = 1 - \int_0^t S_{s} dA_s.\end{equation}\]

We generate a sample of exponentially distributed survival times with independent censoring, and calculate cumulative hazard estimates

n1 <- 300
Samp1 <- rexp(n1,1)
cTimes <- runif(n1)*4
Samp1 <- pmin(Samp1,cTimes)
isNotCensored <- cTimes != Samp1
startTimes1 <- rep(0,n1)

aaMod1 <- aalen(Surv(startTimes1,Samp1,isNotCensored)~1)

We want to estimate the survival function. The cumulative hazard estimates are inserted as a time-ordered matrix where each column is an increment of the vector of cumulatie hazard estimates;

dA_est1 <- diff(c(0,aaMod1$cum[,2]))
hazMatrix <- matrix(dA_est1,nrow=1)

In this case, \(F\) takes the simple form \(F(x) = -x\), and its Jacobian is therefore \(J_{F}(x) = -1\). These must be provided as matrix-valued functions;

F_survival <- function(X)matrix(-X,nrow=1,ncol=1)
F_survival_JacobianList <- list(function(X)matrix(-1,nrow=1,ncol=1))

We obtain plugin estimates of \(S\) and its covariance using the call

param <- pluginEstimate(n1, hazMatrix, F_survival, F_survival_JacobianList,matrix(1,nrow=1,ncol=1),matrix(0,nrow=1,ncol=1))

Finally, we may plot the results with approximate 95 % confidence intervals

times1 <- aaMod1$cum[,1]
plot(times1,param$X,type="s",xlim=c(0,4),xlab="t",ylab="",main="Survival")
lines(times1,param$X + 1.96 * sqrt(param$covariance[1,1,]),type="s",lty=2)
lines(times1,param$X - 1.96 * sqrt(param$covariance[1,1,]),type="s",lty=2)
lines(times1,exp(-times1),col=2)
legend("topright",c("estimated","exact"),lty=1,col=c(1,2),bty="n")

Restricted mean survival

The restricted mean survival function \(R\) solves the system \[\begin{equation} \begin{pmatrix} S_t \\ R_t \end{pmatrix} = \begin{pmatrix}1 \\ 0 \end{pmatrix} + \int_0^t \begin{pmatrix} -S_{s} & 0 \\ 0 & S_s \end{pmatrix}d \begin{pmatrix}A_s \\ s \end{pmatrix},\end{equation}\]

where \(S\) is the survival function. Here, the colums of \(F\), \(F_1\) and \(F_2\), are given by \[\begin{equation} F_1(x_1,x_2) = \begin{pmatrix} -x_1 \\ 0 \end{pmatrix}, F_2(x_1,x_2) = \begin{pmatrix} 0 \\ x_1 \end{pmatrix}.\end{equation}\] The Jacobian matrices are therefore \[\begin{equation} J_{F_1}(x_1,x_2) = \begin{pmatrix} -1 & 0 \\ 0 & 0 \end{pmatrix}, J_{F_2}(x_1,x_2) = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}. \end{equation}\] We define these along with initial values below

F_restrict <- function(X)matrix(c(-X[1],0,0,X[1]),nrow=2)
F_restrict_JacobianList <- list(function(X)matrix(c(-1,0,0,0),nrow=2),
                                function(X)matrix(c(0,0,1,0),nrow=2,byrow=T))

X0_restrict <- matrix(c(1,0),nrow=2)
V0_restrict <- matrix(0,nrow=2,ncol=2)

The restricted mean survival is a ‘regular’ (i.e. Lebesgue) integral, and we must therefore provide the time increments. We choose the time interval \([0,4]\) in \(10^4\) increments:

fimeTimes <- seq(0,4,length.out = 1e4+1)

tms <- sort(unique(c(fimeTimes,times1)))

hazMatrix <- matrix(0,nrow=2,ncol=length(tms))
hazMatrix[1,match(times1,tms)] <- dA_est1
hazMatrix[2,] <- diff(c(0,tms))

We obtain plugin estimates using the call (note the last argument that in this example can be used to improve efficiency);

param <- pluginEstimate(n1, hazMatrix, F_restrict, F_restrict_JacobianList,X0_restrict,V0_restrict,isLebesgue = 2)

We plot the results

plot(tms,param$X[2,],type="s",xlim=c(0,4),ylim=c(0,1.5),xlab="t",ylab="",main="Restricted mean survival")
lines(tms,param$X[2,] + 1.96 * sqrt(param$covariance[2,2,]),type="s",lty=2)
lines(tms,param$X[2,] - 1.96 * sqrt(param$covariance[2,2,]),type="s",lty=2)
lines(tms,1 - exp(-tms),col=2)
legend("topleft",c("estimated","exact"),lty=1,col=c(1,2),bty="n")

Relative survival

We generate another set of exponentially distributed survival times and compare the two groups. A new matrix of cumulative hazard increments must be created;

n2 <- 200
Samp2 <- rexp(n2,1.3)
censTimes2 <- runif(n2)*3
Samp2 <- pmin(Samp2,censTimes2)
isNotCensored2 <- censTimes2 != Samp2
startTimes2 <- rep(0,n2)

aaMod2 <- aalen(Surv(startTimes2,Samp2,isNotCensored2)~1)
dA_est2 <- diff(c(0,aaMod2$cum[,2]))

times2 <- aaMod2$cum[,1]
times <- sort(unique(c(times1,times2)))

hazMatrix <- matrix(0,nrow=2,ncol=length(times))
hazMatrix[1,match(times1,times)] <- dA_est1
hazMatrix[2,match(times2,times)] <- dA_est2

The relative survival function \(RS\) between groups 1 and 2 solve the equation \[\begin{equation} RS_t = 1 + \int_0^t \begin{pmatrix} -RS_s & RS_s\end{pmatrix} d\begin{pmatrix} A_s^1 \\ A_s^2 \end{pmatrix}, \end{equation}\]

i.e. \(F(x) = (-x,x)\). The Jacobians of the columns of \(F\) are \(J_{F_1}(x) = -1\), and \(J_{F_2}(x) = 1\). We specify these functions as follows:

F_relsurv <- function(X)matrix(c(-X,X),nrow=1)
F_relsurv_JacobianList <- list(function(X)matrix(-1,nrow=1),
                               function(X)matrix(1,nrow=1))

The estimates are then obtained by the call

param <- pluginEstimate(n1+n2, hazMatrix, F_relsurv, F_relsurv_JacobianList,matrix(1,nrow=1,ncol=1),matrix(0,nrow=1,ncol=1))

We may plot the results with approximate 95% confidence intervals

plot(times,param$X,type="s",xlim=c(0,4),ylim=c(-0.5,4),xlab="t",ylab="",main="Relative survival")
lines(times,param$X + 1.96 * sqrt(param$covariance[1,1,]),type="s",lty=2)
lines(times,param$X - 1.96 * sqrt(param$covariance[1,1,]),type="s",lty=2)
lines(times,exp(-times*(1/1.3 - 1)),col=2)
legend("topleft",c("estimated","exact"),lty=1,col=c(1,2),bty="n")

Cumulative incidence

We may be in a situation with two competing risks. The cumulative incidences \(C^1,C^2\) solve the system \[\begin{equation} \begin{pmatrix} S_t \\ C_t^1 \\ C_t^2 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} + \int_0^t \begin{pmatrix} -S_s & -S_s \\ S_s & 0 \\ 0 & S_s \end{pmatrix} d \begin{pmatrix} A^1_s \\ A^2_s \end{pmatrix}, \end{equation}\] where \(S\) is the survival. The first columns of \(F\) are

\[\begin{equation}F_1(x_1,x_2,x_2) = \begin{pmatrix} -x_1 \\ x_1 \\ 0 \end{pmatrix}, F_2(x_1,x_2,x_2) = \begin{pmatrix} -x_1 \\ 0 \\ x_1 \end{pmatrix} \end{equation}\]

The reader may verify that the Jacobian matrix of \(F_1\) and \(F_1\) are \[\begin{equation} J_{F_1}(x_1,x_2,x_3) = \begin{pmatrix} -1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} , J_{F_2}(x_1,x_2,x_3) = \begin{pmatrix} -1 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix} \end{equation}\]

We specify the function \(F\) and the list of the Jacobian matrices below, along with the initial values:

F_cuminc <- function(X)matrix(c(-X[1],-X[1],X[1],0,0,X[1]),nrow=3,byrow=T)
F_cuminc_JacobianList <- list(function(X)matrix(c(-1,0,0,1,0,0,0,0,0),nrow=3,byrow=T),
                               function(X)matrix(c(-1,0,0,0,0,0,1,0,0),nrow=3,byrow=T))
X0_cuminc <- matrix(c(1,0,0),nrow=3)
v0_cuminc <- matrix(0,nrow=3,ncol=3)

Now for obtaining hazard estimates. We generate survival times first, before sampling two causes of death (and censoring). We then estimate the cause-specific cumulative hazards, and transform the estimates:

# Generate the data
dfr <- data.frame(from=0,to=rexp(n1+n2,1),from.state = 1)

# Adding causes of death (to.state = 2 or 3) and censoring (to.state = 0)
dfr$to.state <- sample(c(0,2,3),n1+n2,replace=T,prob=c(0.2,0.2,0.6))


# Obtaining cause-specific cumulative hazard estimates and extracting the increments
aamod22 <- aalen(Surv(from,to,to.state %in% c(2)) ~1, data=dfr)
aamod33 <- aalen(Surv(from,to,to.state %in% c(3)) ~1, data=dfr)
times = sort(unique(c(aamod22$cum[,1],aamod33$cum[,1])))
hazMatrix = matrix(0,nrow=2,ncol=length(times))
mt1 = match(aamod22$cum[,1],times)
mt2 = match(aamod33$cum[,1],times)
hazMatrix[1,mt1] = diff(c(0, aamod22$cum[,2] ))
hazMatrix[2,mt2] = diff(c(0, aamod33$cum[,2] ))

# Transforming the estimates
param <- pluginEstimate(n1+n2,hazMatrix,F_cuminc,F_cuminc_JacobianList,X0_cuminc,v0_cuminc)

Finally, we plot the estimates of the cumulative incidences \(C^1\) and \(C^2\) with confidence intervals:

plot(times,param$X[2,],type="s",xlim=c(0,4),ylim=c(0,1),xlab="t",ylab="",main="Cumulative incidence")
lines(times,param$X[2,] + 1.96 * sqrt(param$covariance[2,2,]),type="s",lty=2)
lines(times,param$X[2,] - 1.96 * sqrt(param$covariance[2,2,]),type="s",lty=2)
lines(times,param$X[3,],type="s",col=3)
lines(times,param$X[3,] + 1.96 * sqrt(param$covariance[3,3,]),type="s",lty=2,col=3)
lines(times,param$X[3,] - 1.96 * sqrt(param$covariance[3,3,]),type="s",lty=2,col=3)
legend("topleft",c(expression(C^1),expression(C^2)),lty=1,col=c(1,3),bty="n")

[Transforming cumulative hazard estimates] (Bioaetrika, 2018).↩︎
[On null hypothesis in survival analysis] (Biometrics, 2019)↩︎