Package {HQM}

Time-dependent ROC curve estimation

Description

Inverse Probability of Censoring Weighting (IPCW) estimation of Cumulative/Dynamic time-dependent ROC curve. The function works in the usual survival setting as well as in the competing risks setting. Computation of the iid-representation of areas under time-dependent ROC curves is implemented. This enables computation of inference procedures: Confidence intervals and tests for comparing two AUCs of two different markers measured on the same subjects.

Usage

AUC.hat(T.i, delta.i, pi.i, other.markers = NULL, cause,
	    weighting = "marginal", t, ROC = TRUE, iid = FALSE)

Arguments

Details

This function computes Inverse Probability of Censoring Weighting (IPCW) estimates of Cumulative/Dynamic time-dependent ROC curve. By definition, time-dependent ROC curve intrinsically depends on the definitions of time-dependent cases and controls.
Let T_i denote the event time of the subject i.

Without competing risks : A case is defined as a subject i with T_i \leq t. A control is defined as a subject i with T_i > t.

With competing risks : In this setting, subjects may undergo different type of events, denoted by \delta_i in the following. Let suppose that we are interested in the event \delta_i=1. Then, a case is defined as a subject i with T_i \leq t and \delta_i=1. With competing risks, two definitions of controls were suggested: (i) a control is defined as a subject i that is free of any event, i.e with T_i > t, and (ii) a control is defined as a subject i that is not a case, i.e with T_i > t or with T_i \leq t and \delta_i \neq 1. For all outputs of this package, objects named with .1 refer to definition (i). For instance AUC.hat.1 or se.1 refer to time-dependent area under the ROC curve and its estimated standard error according to the definition (i). Objects named with .2 refer to definition (ii) .

Value

Object of class "ipcwsurvivalROC" or "ipcwcompetingrisksROC", depending on if there is competing risk or not, that is a list. For these classes, there are print, plot and confint methods. Most objects that they contain are similar, but some are specific to each class.

Specific objects of class "ipcwsurvivalROC" :

• AUC.hat : vector of time-dependent AUC estimates at each time points.

• TP : matrix of time-dependent True Positive fraction (sensitivity) estimates.

• FP : matrix of time-dependent False Positive fraction (1-specificity) estimates.

Specific objects of class "ipcwcompetingrisksROC" :

• AUC.hat.1 : vector of time-dependent AUC estimates at each time points with definition (i) of controls (see Details).

• AUC.hat.2 : vector of time-dependent AUC estimates at each time points with definition (ii) of controls (see Details).

• TP : matrix of time-dependent True Positive fraction (sensitivity) estimates.

• FP.1 : matrix of time-dependent False Positive fraction (1-specificity) estimates with definition (i) of controls (see Details).

• FP.2 : matrix of time-dependent False Positive fraction (1-specificity) estimates with definition (ii) of controls (see Details).

Objects common to both classes :

• t : the time points for which the time-dependent ROC curves were computed.

• G.hat : a object of class "IPCW", containing all informations about the weights. See ipcw function of pec package.

• computation.time : the total computation time.

• cum.incidence : the vector of estimated probabilities of being a case at each time points.

• S.t : the vector of estimated probabilities of being event-free at each time points.

• counts.at.t : a matrix containing descriptive statistics at each time points (like numbers of observed cases or censored observations before each time points).

• iid : the logical value of parameter iid used in argument.

• n : the sample size, after having omitted missing vaues.

• inference : a list that contains, among other things, iid-representations and estimated standard errors of the estimators, and that is used for computation of comparison tests and confidence intervals.

• computation.time : the computation time, in seconds.

Author(s)

Paul Blanche pabl@sund.ku.dk

References

Hung, H. and Chiang, C. (2010). Estimation methods for time-dependent AUC with survival data. Canadian Journal of Statistics, 38(1):8-26

Uno, H., Cai, T., Tian, L. and Wei, L. (2007). Evaluating prediction rules for t-years survivors with censored regression models. Journal of the American Statistical Association, 102(478):527-537.

Blanche, P., Dartigues, J. F., & Jacqmin-Gadda, H. (2013). Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Statistics in medicine, 32(30), 5381-5397.

P. Blanche, A. Latouche, V. Viallon (2013). Time-dependent AUC with right-censored data: A Survey. Risk Assessment and Evaluation of Predictions, 239-251, Springer, https://arxiv.org/abs/1210.6805.

Examples

library(survival)
Landmark <- 2
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]
n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))
 
timesS2 <- seq(Landmark,14,by=0.5)
b=0.9
h.hat.x<- get.h.hat.x(pbcT1, 'albumin', s.grid = s.grid,  T.name = 'years',  
                s.name = 'year', delta.name = 'status2', 2, 0.9) 
s.grid.2  = seq(Landmark, 14,  length=n.grid-1)
S.hat.alb<- make.S.hat(    (s.grid.2[2]-s.grid.2[1])/4 , h.hat.x)
t <- 1.5
AUC.hat.s.t<-AUC.hat.HQM(pbcT1, S.hat.alb, Landmark,t,  
              T.name = 'years', delta.name = 'status2')
AUC.hat.s.t

AUC for the High Quality Marker estimator

Description

Calculates the AUC for the HQM estimator.

Usage

AUC.hat.HQM(data, S.hat, s, t, T.name, delta.name)

Arguments

Details

The function AUC.hat.HQM implements the AUC calculation for the HQM estimator estimator.

Value

A vector of two values: the landmark time of the calculation and the AUC value.

See Also

BS.hat.HQM

Examples

library(survival)
Landmark <- 2
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]
n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))
 
timesS2 <- seq(Landmark,14,by=0.5)
b=0.9
h.hat.x<- get.h.hat.x(pbcT1, 'albumin', s.grid = s.grid,  T.name = 'years',  
                s.name = 'year', delta.name = 'status2', 2, 0.9) 
s.grid.2  = seq(Landmark, 14,  length=n.grid-1)
S.hat.alb<- make.S.hat(    (s.grid.2[2]-s.grid.2[1])/4 , h.hat.x)
t <- 1.5
AUC.hat.s.t<-AUC.hat.HQM(pbcT1, S.hat.alb, Landmark,t,  
              T.name = 'years', delta.name = 'status2')
AUC.hat.s.t

Brier score for the High Quality Marker estimator

Description

Calculates the Brier score for the HQM estimator.

Usage

BS.hat.HQM(data, S.hat, s, t, T.name, delta.name)

Arguments

Details

The function BS.hat.HQM implements the Brier score calculation for the HQM estimator estimator.

Value

Scalar: the Brier score of the HQM estimator.

See Also

AUC.hat.HQM

Examples

library(pec)
library(survival)
Landmark <- 2

pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]
n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))
 
timesS2 <- seq(Landmark,14,by=0.5)


b=0.9
h.hat.x<- get.h.hat.x(pbcT1, 'albumin', s.grid = s.grid, T.name = 'years',  
               s.name = 'year', delta.name = 'status2', 2, 0.9) 
s.grid.2  = seq(Landmark, 14,  length=n.grid-1)
S.hat.alb<- make.S.hat(    (s.grid.2[2]-s.grid.2[1])/4 , h.hat.x)


t <- 1.5
BS.hat.s.t<-BS.hat.HQM(pbcT1, S.hat.alb, Landmark,t, 
                T.name = 'years', delta.name = 'status2')
BS.hat.s.t

Cross validation bandwidth selection

Description

Implements the bandwidth selection for the future conditional hazard rate \hat h_x(t) based on K-fold cross validation.

Usage

CV.b(data,   X.name, T.name = 'years',
            s.name = 'year', delta.name = 'status2', size.I.j, b.candidates)

Arguments

Details

The function CV.b implements the cross validation bandwidth selection for the future conditional hazard rate \hat h_x(t) given by

b_{CV} = arg min_b \sum_{i = 1}^N \int_0^T \int_s^T Z_i(t)Z_i(s)(\hat{h}_{X_i(s)}(t-s)- h_{X_i(s)}(t-s))^2 dt ds,

where \hat h_x(t) is a smoothed kernel density estimator of h_x(t) and Z_i the exposure process of individual i. Note that \hat h_x(t) is dependent on b.

Value

A list with the tested bandwidths and its cross validation scores.

References

Bagkavos, I., Isakson, R., Mammen, E., Nielsen, J., and Proust–Lima, C. (2025). Biometrika, 112(2), asaf008. doi:10.1093/biomet/asaf008

See Also

g.hat.i, Q1, R.hat.K, h.hat.minus.I.j, data.minus.I.j

Examples

size.I.j = 26
# For Albumin marker:
b.candidates = seq(0.9, 1.7, 0.1)


CV.b.alb = CV.b(pbc2,  'albumin', 'years', 'year', 'status2', size.I.j, b.candidates)
CV.b.alb[[2]][which.min(CV.b.alb[[1]])]

# For Bilirubin marker:
b.candidates = seq(3, 4, 0.1)
CV.b.bil = CV.b(pbc2, 'serBilir', 'years', 'year', 'status2', size.I.j, b.candidates)
CV.b.bil[[2]][which.min(CV.b.bil[[1]])]
CV.b.bil

Cross validation index parameter selection

Description

Implements the index parameter selection for two markers based on K-fold cross validation.

Usage

CV.theta(theta, data, X.1.name, X.2.name,  
          T.name = 'years', s.name = 'year', delta.name = 'status2', size.I.j, b)

Arguments

Details

The function CV.theta implements the cross validation index parameter selection for the indexing of two markers and given by

\hat \theta = arg min_{\theta_1, \theta_2} \sum_{i = 1}^N \int_0^T \int_s^T Z_i(t)Z_i(s)(\hat{h}_{\theta_0^T X_i(s)}(t-s)- h_{\theta_0^T X_i(s)}(t-s))^2 dt ds,

where \hat h_x(t) is the HQM estimator of h_x(t), see Bagkavos et al. (2025), doi:10.1093/biomet/asaf008, and Z_i the exposure process of individual i. Note that \hat h_x(t) is dependent on b.

Value

A list with the tested parameters and its cross validation scores.

References

Bagkavos, I., Isakson, R., Mammen, E., Nielsen, J., and Proust–Lima, C. (2025). Biometrika, 112(2), asaf008. doi:10.1093/biomet/asaf008

See Also

g.hat.i, Q1, R.hat.K, h.hat.minus.I.j, data.minus.I.j

Examples

# Obtain the indexing parameters for the combination of albumin and bilirubin markers
# These are the values used in the example of function h.hat.x.t.vec
# and yielded the values:  (theta.alb, theta.bil) = ( 0.0702, 0.0856 ) that were used there.  

size.I.j = 26
b.candidates = seq(1.1, 1.2, by=0.1)

for(i in 1:length(b.candidates))
{
  theta.fit<- optim(c(0.5, 0.5), CV.theta,  data=pbc2, X.1.name='albumin', 
              X.2.name= 'serBilir', T.name = 'years', s.name = 'year', 
              delta.name = 'status2', size.I.j=26, b=b.candidates[i], method="Nelder-Mead")
  cat("i= ", i, " ", theta.fit$par, " ", theta.fit$value, "count calls to fn = ", 
      theta.fit$counts, " converge? ", theta.fit$convergence,   "\n")
  theta.fit
}

Epanechnikov kernel

Description

Implements the Epanechnikov kernel function

Usage

K.Epanechnikov(u)

Arguments

Details

Implements the Epanechnikov kernel function

K(x) = \frac{3}{4}(1-x^2)*(|x|<1)),

Value

Scalar, the value of the Epanechnikov kernel at u.

Classical (unmodified) kernel and related functionals

Description

Implements the classical kernel function and related functionals

Usage

K.b(b,u,v, K)
c1.u.integrand(b,u,v, K)
K.b.matrix(b,u,v, K)

Arguments

Details

The function K.b implements the classical kernel function calculation

h^{-1} K \left ( \frac{x-y}{h} \right )

for scalars x and y while c1.u.integrand implements the functional

h^{-1} K \left ( \frac{x-y}{h} \right )(x-y)

again for for scalars x and y. The function K.b.matrix is the vectorized version of K.b. It uses as inputs the vectors (X_1, \dots, X_n) and (Y_1, \dots, Y_n) and returns a n \times n matrix with entries

h^{-1} K \left ( \frac{X_i-Y_j}{h} \right )

Value

Scalar values for K.b and c1.u.integrand and matrix outputs for K.b.matrix.

Local linear kernel

Description

Implements the local linear kernel function.

Usage

K.b.u.v(b,u,v, K)

Arguments

Details

Implements the local linear kernel

K_{x,b}(u)= \frac{K_b(u)-K_b(u)u^T D^{-1}c_1}{c_0 - c_1^T D^{-1} c_1},

where c_1 = (c_{11}, \dots, c_{1d})^T, D = (d_{ij})_{(d+1) \times (d+1)} with

c_0 = \sum_{i=1}^n \int_0^T K_b(x-X_i(s)) Z_i(s)ds,

c_{ij} = \sum_{i=1}^n \int_0^T K_b(x-X_i(s))\{x-X_{ij}(s)\} Z_i(s)ds,

d_{jk} = \sum_{i=1}^n \int_0^T K_b(x-X_i(s))\{x-X_{ij}(s)\}\{x-X_{ik}(s)\} Z_i(s)ds,

see also Nielsen (1998), doi:10.1080/03461238.1998.10413997.

Value

Matrix output with entries the values of the kernel function at each point.

References

Nielsen (1998), Marker dependent kernel hazard estimation from local linear estimation, Scandinavian Actuarial Journal, pp. 113-124. doi:10.1080/03461238.1998.10413997

Indexing parameter objective function

Description

Objective used for optimizing indexing parameters via optim.

Usage

Q.theta(theta, data, data.id, s.grid, X.obs.list, X.i.T.i.list,  T.name = 'years', 
            s.name = 'year', delta.name = 'status2', b, x, h.true)

Arguments

Value

Scalar cross-validation score.

Examples

 

X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'


theta.1  <- 0.0702 #0.149
theta.2 <- 0.0856 #0.10
x.1 = 0 # refers to zero mean variables - slightly high
x.2 = 1.9 # refers to zero mean variable - high
b = 0.42 
x = theta.1 * x.1 + theta.2 *x.2


# first simulate true HR function:
data <- pbc2[,c(id.name, X.1.name, X.2.name, T.name, s.name, delta.name)]
n.s.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- data[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))

n.obs <- length(data$id)
n<-max(  as.double(data[,'id']) )
###################  Create bootstrap samples by group ####################
set.seed(1)  
B<- 10 # 400 #50
bootstrap.samples<-list()
for(j in 1:B)
{
  rows.star<-c()
  ids.star<-c()
  ids.drawn <- sample (n, replace = TRUE)  
  for(l in 1:n)
  {
    rows.l<-which(data[,id.name]==ids.drawn[l])
    rows.star<-c(rows.star, rows.l)
    ids.l<-rep(ids.drawn[l], times=length(rows.l))
    ids.star<-c(ids.star, ids.l)
  }
  data.star.unsorted<-data[rows.star,]
  data.star.unsorted<-data[rows.star,]
  bootstrap.samples[[j]]<- data.star.unsorted[order(data.star.unsorted$id),] #data[rows.star,]
}
h.true<- h.true.simulated(bootstrap.samples, id.name='id', n.s.grid=n.s.grid, 
              X.1.name=X.1.name, X.2.name= X.2.name, 
              T.name = T.name, 
              s.name = s.name,  delta.name = delta.name, 
              theta = c(theta.1,  theta.2), b)
##########################################################################

# Then run the optimization for the indexing parameters:

data.use<-bootstrap.samples[[1]]
data.use.id<-to.id(data.use)
data.use.id<-data.use.id[complete.cases(data.use.id), ]
X.1.obs=data.use[,X.1.name] -mean(data.use[, X.1.name])
X.1.i.T.i=data.use.id[,X.1.name] -mean(data.use.id[, X.1.name])
X.2.obs=data.use[,X.2.name]  -mean(data.use[, X.2.name])
X.2.i.T.i=data.use.id[,X.2.name] -mean(data.use.id[, X.2.name])

X.obs.list=list(X.1.obs, X.2.obs)
X.i.T.i.list=list(X.1.i.T.i, X.2.i.T.i)

s.obs<- data.use[, s.name]
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
theta.fit<- optim(par=c(theta.1, theta.2), fn=Q.theta,  data=data.use, data.id=data.use.id, 
            s.grid=s.grid, X.obs.list=X.obs.list, X.i.T.i.list=X.i.T.i.list,  T.name = T.name, 
            s.name = s.name, delta.name = delta.name, b=b, x=x, 
            h.true=h.true, method="Nelder-Mead") 

 

Bandwidth selection score Q1

Description

Calculates a part for the K-fold cross validation score.

Usage

Q1(h.hat.matrix, X.pos, n.s.grid, n, Z.i.ds)

Arguments

Details

The function implements

Q_1 = \sum_{i = 1}^N \int_0^T \int_s^T Z_i(t)Z_i(s)\hat{h}_{X_i(s)}^2(t-s) dt ds,

where \hat{h} is the hqm estimator, Z the exposure and X the marker.

Value

A value of the score Q1.

See Also

CV.b

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN <- make.dN(pbc2, data.id, X.grid, s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
            n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)

b = 1.7
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
              n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)
dN.i  <- make.dN.i(s.grid, n.s.grid, T.i, delta.i, n)

t = 2

h.hat.matrix = t(sapply(s.grid[1:99],
  function(t.point){h.hat.x.t.vec(X.grid, s.grid, n.s.grid, alpha.hat, t, b, Z.i.ds, X.pos, n)}))

Q1.b = Q1(h.hat.matrix, X.pos, n.X.grid, n, Z.i.ds)

Bandwidth selection score R

Description

Calculates a part for the K-fold cross validation score.

Usage

R.hat.K(h.hat.minus.I.j.list, X.pos, n.s.grid, Z.i.ds, dN.i, n)

Arguments

Details

The function implements the estimator

\hat{R}_K = \sum_{j = 1}^K\sum_{i \in I_j} \int_0^T g^{-I_j}_i(t) dN_i(t),

where \hat{g}^{-I_j}_i(t) = \int_0^t Z_i(s) \hat{h}^{-I_j}_{X_i(s)}(t-s) ds, and \hat{h}^{-I_j} is estimated without information from all counting processes i with i \in I_j. This function estimates

R = \sum_{i = 1}^N \int_0^T \int_s^T Z_i(t)Z_i(s)\hat{h}_{X_i(s)}(t-s) h_{X_i(s)}(t-s) dt ds .

where \hat{h} is the hqm estimator, Z the exposure and X the marker.

Value

A matrix with \hat{g}^{-I_j}_i(t) for all individuals i and time grid points t.

See Also

CV.b

Examples

data.id = to.id(pbc2)
n = max(as.numeric(pbc2$id))
b = 1.5
size.I.j = 104
h.hat.minus.I.j.list = h.hat.minus.I.j(pbc2, data.id, 'serBilir', 'years', 
                                          'year', 'status2', n, size.I.j, b)


n.s.grid <- n.X.grid <- 100
s.obs = pbc2$year
X.obs = pbc2$serBilir
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

T.i <- data.id$years
delta.i <- data.id$status2

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)
X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
              n.s.grid, n.X.grid, s.pos, X.pos, 'years', n)
dN.i  <- make.dN.i(s.grid, n.s.grid, T.i, delta.i, n)

two.R.hat.K = R.hat.K(h.hat.minus.I.j.list, X.pos, n.X.grid, Z.i.ds, dN.i, n)
two.R.hat.K

Linear interpolation

Description

Implements a linear interpolation between observered marker values.

Usage

X.i.s.interpolate(s.grid, ids, id.obs, X.obs, s.obs)

Arguments

Details

Given time points t_1,...,t_K and marker values m_1,...,m_J at different time points t^m_1,...,t^m_J, the function calculates a linear interpolation f with f(t^m_i) = m_i at the time points t_1,...,t_K for all indicated individuals. Returned are then (f(t_1),...,f(t_K)). Note that the first value is always observed at time point 0 and the function f is extrapolated constantly after the last observed marker value.

Value

A matrix with columns (f(t_1),...,f(t_K)) as described above for every individual in the vector ids.

Examples

n.s.grid <- 100
X.obs = pbc2$serBilir
s.obs = pbc2$year
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
data.id = to.id(pbc2)

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

Compute Pivot Pointwise Confidence Intervals for the Indexed Hazard Rate Estimate

Description

Computes bootstrap pivot and symmetric bootstrap-pivot confidence intervals for the hazard, log-hazard, and back-transformed (log-scale) hazard rate functions, using the indexed hazard estimator.

Usage

  bootstrap.CI.pivotal(B,  n.grid, h.hat.star.matrix, t.grid, h.hat.x, alpha.level)

Arguments

Details

This function computes several forms of pivot confidence intervals for indexed hazard rate estimates. Let \hat \sigma be the sample standard deviation of the bootstrap estimators \hat{h}_{x}^{(j)}(t), j=1,\dots, B and let k_{\alpha/2}^p, k_{1-\alpha/2}^p and \bar k_{1-\alpha}^p be the \alpha/2, 1-\alpha/2 and 1-\alpha quantiles respectively of

\frac{\hat{h}_{x}^{(j)}(t) - \hat{h}_{x}(t)}{\hat \sigma},\; j=1,\dots, B.

Then, the pivot bootstrap confidence interval (CI) for \hat{h}_{x}(t) is given by

\Bigg ( \hat{h}_{x}(t) - \hat \sigma k_{1-\alpha/2}^p, \hat{h}_{x}(t) - \hat \sigma k_{\alpha/2}^p \Bigg ).

The symmetric pivot bootstrap CI for \hat{h}_{x}(t) is

\Bigg ( \hat{h}_{x}(t) - \hat \sigma \bar k_{1-\alpha}^p, \hat{h}_{x}(t) + \hat \sigma\bar k_{1-\alpha}^p \Bigg ).

For the confidence intervals for the logarithm of the hazard rate function, first let

L_x(t) = \log \left \{ h_x(t)\right \}, \hat L_{x}(t) = \log \left \{ \hat h_{x}(t)\right \}, \hat L^{(j)}_{x}(t) = \log \left \{ \hat h_{x}^{(j)}(t)\right \},

and k_{\alpha/2}^{L,p}, k_{1-\alpha/2}^{L,p} and \bar k_{1-\alpha}^{L,p} be respectively the \alpha/2, 1-\alpha/2 and 1-\alpha quantile of

\frac{|\hat{L}_{x}^{(j)}(t) - \hat{L}_{x}(t)|}{\hat \sigma}, j=1,\dots, B.

For the log hazard function L_x(t) a pivotal bootstrap confidence interval is

\Bigg ( \hat{L}_{x}(t) - \hat \sigma k_{1-\alpha/2}^{L,p}, \hat{L}_{x}(t) - \hat \sigma k_{\alpha/2}^{L,p} \Bigg ).

A symmetric pivotal bootstrap CI for the log hazard is

\Bigg ( \hat{L}_{x}(t) - \hat \sigma \bar k_{1-\alpha}^{L,p}, \hat{L}_{x}(t) + \hat \sigma \bar k_{1-\alpha}^{L,p} \Bigg ).

These last two confidence intervals can be transformed back to yield confidence intervals for the hazard rate function h_x(t). These are:

\Bigg ( \hat{h}_{x}(t) e^{- \hat \sigma k_{1-\alpha/2}^{L,p}}, \hat{h}_{x}(t) e^{- \hat \sigma k_{\alpha/2}^{L,p}} \Bigg ),

and

\Bigg ( \hat{h}_{x}(t) e^{- \hat \sigma \bar k_{1-\alpha}^{L,p}}, \hat{h}_{x}(t) e^{ \hat \sigma \bar k_{1-\alpha}^{L,p}} \Bigg )

respectively.

Note: The bootstrap matrix h.hat.star.matrix is assumed to contain estimates produced using the same time grid as t.grid and the same estimator used to generate h.hat.x.

Value

A data frame with the following columns:

See Also

bootstrap.h.hat.star, h.hat.star.x.theta

Examples

 
X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'



theta.1  <- 0.0702 #0.149
theta.2 <- 0.0856 #0.10
x.1 = 0 # refers to zero mean variables - slightly high
x.2 = 1.9 # refers to zero mean variable - high
b = 0.42#theta.alb * b.alb + theta.bil *b.bil # 7
x = theta.1 * x.1 + theta.2 *x.2
n.s.grid<-50

#Store original sample values:
data <- pbc2[,c(id.name, X.1.name, X.2.name, T.name, s.name, delta.name)]
n.obs <- length(data$id)
n<-max(  as.double(data[,'id']) )

data.id <- to.id(data)

X.1.obs=data[,X.1.name] -mean(data[, X.1.name])
X.1.i.T.i=data.id[,X.1.name] -mean(data.id[, X.1.name])
X.2.obs=data[,X.2.name]  -mean(data[, X.2.name])
X.2.i.T.i=data.id[,X.2.name] -mean(data.id[, X.2.name])

X.obs.list=list(X.1.obs, X.2.obs)
X.i.T.i.list=list(X.1.i.T.i, X.2.i.T.i)

# Calculate the indexed HQM estimator on the original sample:
h.hat.x.theta.output<- get.h.hat.x.theta(pbc2, id.name, n.s.grid,  X.obs.list, X.i.T.i.list, 
                  T.name = 'years', s.name = 'year', delta.name = 'status2', 
                  theta= c(theta.1,  theta.2), b, x)
                          
h.hat.x<-h.hat.x.theta.output[,2] # Indexed HQM estimator on original sample
t.grid<-h.hat.x.theta.output[,1] # evaluation grid points
n.grid<- n.s.grid # length(h.hat.x)

#  Create bootstrap samples by group 
set.seed(1)  
B<- 50  #for display purposes only; for sensible results use B=1000 (slower) 
bootstrap.samples<-list()
for(j in 1:B)
{
  rows.star<-c()
  ids.star<-c()
  ids.drawn <- sample (n, replace = TRUE)  
  for(l in 1:n)
  {
    rows.l<-which(data[,id.name]==ids.drawn[l])
    rows.star<-c(rows.star, rows.l)
    ids.l<-rep(ids.drawn[l], times=length(rows.l))
    ids.star<-c(ids.star, ids.l)
  }
  data.star.unsorted<-data[rows.star,]
  data.star.unsorted<-data[rows.star,]
  bootstrap.samples[[j]]<- data.star.unsorted[order(data.star.unsorted$id),]  
}     

# Simulate true hazard rate function:   
h.true<- h.true.simulated(bootstrap.samples, id.name='id', n.grid, 
              X.1.name=X.1.name, X.2.name= X.2.name, 
              T.name = T.name, s.name = s.name,  
              delta.name = delta.name, theta = c(theta.1,  theta.2), b)

# Bootstrap the original indexed HQM estimator: 
h.hat.star.matrix <- bootstrap.h.hat.star(B, bootstrap.samples, X.1.name, X.2.name, 
                  T.name, s.name,  delta.name, b = 0.4, x = x, 
                  h.true = h.true, theta.start = c(0.07, 0.08), n.grid)

# Construct Ci's: 
alpha.level<-0.05   
CIs.pivotal<- bootstrap.CI.pivotal(B, n.grid, h.hat.star.matrix, t.grid, h.hat.x, alpha.level)

# extract Plain   + symmetric CIs 
CI.upper<-CIs.pivotal[,"CI.upper"]
CI.lower<-CIs.pivotal[,"CI.lower"]
CI.sym.upper<-CIs.pivotal[,"CI.sym.upper"]
CI.sym.lower<-CIs.pivotal[,"CI.sym.lower"]

J <- 80 #select the first 80 grid points (for display purposes only)
#Plot the selected CIs
plot(t.grid[1:J],   h.hat.x[1:J], type="l", ylim=c(0,2), ylab="Hazard rate", 
      xlab="time", lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.upper[1:J], rev(CI.lower[1:J])), 
        col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.sym.lower[1:J], lty=2, lwd=2)   

# extract transformed from Log HR + symmetric CIs 
CI.exp.upper<-CIs.pivotal[,"CI.exp.upper"]
CI.exp.lower<-CIs.pivotal[,"CI.exp.lower"]
CI.exp.sym.upper<-CIs.pivotal[,"CI.exp.sym.upper"]
CI.exp.sym.lower<-CIs.pivotal[,"CI.exp.sym.lower"]

#Plot the selected CIs
plot(t.grid[1:J],   h.hat.x[1:J] , type="l", ylim=c(0,2), ylab="Log Hazard rate", 
      xlab="time", lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.exp.upper[1:J], rev(CI.exp.lower[1:J])), 
        col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.exp.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.exp.sym.lower[1:J], lty=2, lwd=2)  

# extract Log HR + symmetric CIs  
CI.L.upper<-CIs.pivotal[,"CI.L.upper"]
CI.L.lower<-CIs.pivotal[,"CI.L.lower"]
CI.L.sym.upper<-CIs.pivotal[,"CI.L.sym.upper"]
CI.L.sym.lower<-CIs.pivotal[,"CI.L.sym.lower"]

#Plot the selected CIs
plot(t.grid[1:J], log(h.hat.x[1:J] ), type="l", ylim=c(-5,4), ylab="Hazard rate", 
      xlab="time",  lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.L.upper[1:J], rev(CI.L.lower[1:J])), 
        col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.L.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.L.sym.lower[1:J], lty=2, lwd=2)  

Compute Quantile Pointwise Confidence Intervals for for the Indexed Hazard Rate Estimate

Description

Computes quantile-bootstrap confidence intervals and its symmetric counterparts for the hazard rate, log-hazard rate, and back-transformed (from log scale) hazard rate functions, based on the indexed hazard estimator.

Usage

  bootstrap.CI.quantile(B, n.grid, h.hat.star.matrix, t.grid, h.hat.x, alpha.level)

Arguments

Details

This function computes several forms of pivot confidence intervals for indexed hazard rate estimates. Set k_{\alpha/2} = \hat{h}_{x}^{[\alpha/2]}(t) - \hat{h}_{x}(t) and k_{1-\alpha/2} = \hat{h}_{x}^{[1-\alpha/2]}(t) - \hat{h}_{x}(t) where \hat{h}_{x}^{[\alpha/2]}(t) is the \alpha/2 quantile of \hat{h}_{x}^{(j)}(t), j=1,\dots,B, obtained by ordering the bootstrap estimators in ascending order and selecting the \alpha/2-th ordered value. For example, for B=1000 bootstrap iterations and \alpha=0.05, \hat{h}_{x}^{[\alpha/2]}(t) will be the 25th smallest out of the 1000 values \hat{h}_{x}^{(j)}(t), j=1,\dots,1000. Also denote with \bar k_{1-\alpha} the 1-\alpha quantile of

| \hat{h}_{x}^{(j)}(t) - \hat{h}_{x}(t)|, j=1,\dots, B.

Then, the quantile bootstrap CI for \hat{h}_{x}(t) is given by

\Bigg ( \hat{h}_{x}(t) - k_{1-\alpha/2}, \hat{h}_{x}(t) - k_{\alpha/2} \Bigg ).

The symmetric quantile CI (basic CI) is

\Bigg ( \hat{h}_{x}(t) - \bar k_{1-\alpha}, \hat{h}_{x}(t) + \bar k_{1-\alpha} \Bigg ).

The confidence intervals for the logarithm of the hazard rate function are defined as follows. First set k_{\alpha/2}^L = \hat{L}_{x}^{[\alpha/2]}(t) - \hat{L}_{x}(t) and k_{1-\alpha/2}^L = \hat{L}_{x}^{[1-\alpha/2]}(t) - \hat{L}_{x}(t).

Also denote with \bar k_{1-\alpha}^L the 1-\alpha quantile of | \hat{L}_{x}^{(j)}(t) - \hat{L}_{x}(t)|, j=1,\dots, B. For the log hazard function L_x(t) we have the quantile confidence interval is

\Bigg ( \hat{L}_{x}(t) - k_{1-\alpha/2}^L, \hat{L}_{x}(t) - k_{\alpha/2}^L \Bigg ).

The corresponding symmetric quantile CI for the log hazard is

\Bigg ( \hat{L}_{x}(t) - \bar k_{1-\alpha}^L, \hat{L}_{x}(t) + \bar k_{1-\alpha}^L \Bigg ).

These confidence intervals are transformed back to confidence intervals for the hazard rate function h_x(t):

\Bigg ( \hat{h}_{x}(t) e^{- k_{1-\alpha/2}^L}, \hat{h}_{x}(t) e^{- k_{\alpha/2}^L} \Bigg ).

The corresponding symmetric confidence interval is

\Bigg ( \hat{h}_{x}(t) e^{- \bar k_{1-\alpha}^L}, \hat{h}_{x}(t) e^{ \bar k_{1-\alpha}^L} \Bigg ).

Note: The bootstrap matrix h.hat.star.matrix is assumed to contain estimates produced using the same time grid as t.grid and the same estimator used to generate h.hat.x.

Value

A data frame with the following columns:

See Also

bootstrap.h.hat.star, h.hat.star.x.theta

Examples

  
X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'


theta.1  <- 0.0702 #0.149
theta.2 <- 0.0856 #0.10
x.1 = 0 # refers to zero mean variables - slightly high
x.2 = 1.9 # refers to zero mean variable - high
b = 0.42# The result, on the indexed marker 'indmar' of 
         #\code{CV.b(pbc2, 'indmar', 'years', 'year', 'status2', size.I.j=26, seq(0.2,0.4,by=0.01))} 
x = theta.1 * x.1 + theta.2 *x.2
n.s.grid<-50

#Store original sample values:
data <- pbc2[,c(id.name, X.1.name, X.2.name, T.name, s.name, delta.name)]
n.obs <- length(data$id)
n<-max(  as.double(data[,'id']) )

data.id <- to.id(data)

X.1.obs=data[,X.1.name] -mean(data[, X.1.name])
X.1.i.T.i=data.id[,X.1.name] -mean(data.id[, X.1.name])
X.2.obs=data[,X.2.name]  -mean(data[, X.2.name])
X.2.i.T.i=data.id[,X.2.name] -mean(data.id[, X.2.name])

X.obs.list=list(X.1.obs, X.2.obs)
X.i.T.i.list=list(X.1.i.T.i, X.2.i.T.i)

# Calculate the indexed HQM estimator on the original sample:
h.hat.x.theta.output<- get.h.hat.x.theta(pbc2, id.name, n.s.grid,  X.obs.list, 
                          X.i.T.i.list, T.name = 'years', s.name = 'year',  
                          delta.name = 'status2', 
                         theta= c(theta.1,  theta.2), b, x)
                          
h.hat.x<-h.hat.x.theta.output[,2] # Indexed HQM estimator on original sample
t.grid<-h.hat.x.theta.output[,1] # evaluation grid points
n.grid<- n.s.grid # length(h.hat.x)

#  Create bootstrap samples by group: 
set.seed(1)  
B<- 50   #for display purposes only; for sensible results use B=1000 (slower) 
bootstrap.samples<-list()
for(j in 1:B)
{
  rows.star<-c()
  ids.star<-c()
  ids.drawn <- sample (n, replace = TRUE)  
  for(l in 1:n)
  {
    rows.l<-which(data[,id.name]==ids.drawn[l])
    rows.star<-c(rows.star, rows.l)
    ids.l<-rep(ids.drawn[l], times=length(rows.l))
    ids.star<-c(ids.star, ids.l)
  }
  data.star.unsorted<-data[rows.star,]
  data.star.unsorted<-data[rows.star,]
  bootstrap.samples[[j]]<- data.star.unsorted[order(data.star.unsorted$id),]  
}     
 
# Simulate true hazard rate function:    
h.true<- h.true.simulated(bootstrap.samples, id.name='id', n.grid, 
              X.1.name=X.1.name, X.2.name= X.2.name, 
              T.name = T.name, s.name = s.name, 
              delta.name = delta.name, theta = c(theta.1, theta.2), b)

# Bootstrap the original indexed HQM estimator: 
h.hat.star.matrix <- bootstrap.h.hat.star(B, bootstrap.samples, X.1.name, X.2.name, T.name,
                            s.name,  delta.name, b = 0.4, x = x, h.true = h.true, 
                            theta.start = c(0.07, 0.08), n.grid   )
J <- 80
alpha.level<-0.05   

# Construct Ci's: 
CIs.quantile<- bootstrap.CI.quantile(B, n.grid, h.hat.star.matrix, t.grid, h.hat.x, alpha.level)

# extract Plain   + symmetric CIs and plot them:
CI.upper<-CIs.quantile[,"CI.upper"]
CI.lower<-CIs.quantile[,"CI.lower"]
CI.sym.upper<-CIs.quantile[,"CI.sym.upper"]
CI.sym.lower<-CIs.quantile[,"CI.sym.lower"]

#Plot the selected CIs
plot(t.grid[1:J],   h.hat.x[1:J], type="l", ylim=c(0,2), ylab="Hazard rate", 
      xlab="time", lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.upper[1:J], rev(CI.lower[1:J])), 
        col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.sym.lower[1:J], lty=2, lwd=2)   

# extract transformed from Log HR + symmetric CIs 
CI.exp.upper<-CIs.quantile[,"CI.exp.upper"]
CI.exp.lower<-CIs.quantile[,"CI.exp.lower"]
CI.exp.sym.upper<-CIs.quantile[,"CI.exp.sym.upper"]
CI.exp.sym.lower<-CIs.quantile[,"CI.exp.sym.lower"]

#Plot the selected CIs
plot(t.grid[1:J],   h.hat.x[1:J], type="l", ylim=c(0,2), ylab="Hazard rate", 
      xlab="time", lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.exp.upper[1:J], rev(CI.exp.lower[1:J])), 
        col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.exp.sym.upper[1:J], lty=2, lwd=2 )#, lwd=3
lines(t.grid[1:J],  CI.exp.sym.lower[1:J], lty=2, lwd=2)  

# extract Log HR + symmetric CIs 
CI.L.upper<-CIs.quantile[,"CI.L.upper"]
CI.L.lower<-CIs.quantile[,"CI.L.lower"]
CI.L.sym.upper<-CIs.quantile[,"CI.L.sym.upper"]
CI.L.sym.lower<-CIs.quantile[,"CI.L.sym.lower"]

#Plot the selected CIs
plot(t.grid[1:J],  log(h.hat.x[1:J]), type="l", ylim=c(-5,4), ylab="Log Hazard rate", 
     xlab="time", lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.L.upper[1:J], rev(CI.L.lower[1:J])), 
      col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.L.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.L.sym.lower[1:J], lty=2, lwd=2)  
   

Compute Studentized Double Bootstrap Pointwise Confidence Intervals for the Indexed Hazard Rate Estimate

Description

Computes bootstrap confidence intervals—studentized (double bootstrap) and its symmetric versions for the hazard rate, log-hazard rate, and back-transformed (from the log scale) hazard rate functions, based on the indexed hazard estimator.

Usage

  bootstrap.CI.studentized(n.grid, double.bootstrap.output, t.grid, h.hat.x, alpha.level)

Arguments

Details

This function computes several forms of studentized confidence intervals for the indexed hazard rate function. First, for each bootstrap iteration j=1,\dots, B we construct estimators of the standard deviation, denoted by \hat \sigma_j^2 as follows: each bootstrap sample is itself bootstrapped, say B_1 times, we estimate the corresponding B_1 index parameters \hat \theta_{j,k} = (\hat \theta_{1,j,k}, \hat \theta_{2,j,k}), k=1,\dots, B_1, and use them to calculate the corresponding hazard estimators \hat{h}_{x,}^{(j,k)}(t), k=1,\dots, B_1. Finally we calculate \hat \sigma_j^2 as the sample variance of \hat{h}_{x}^{(j,1)}(t), \dots, \hat{h}_{x}^{(j,B_1)}(t). Let k_{\alpha/2}^s, k_{1-\alpha/2}^s and \bar k_{1-\alpha}^s be respectively the \alpha/2, 1-\alpha/2 and 1-\alpha quantiles of

\frac{|\hat{h}_{x}^{(j)}(t) - \hat{h}_{x}(t)|}{\hat \sigma_j},\; j=1,\dots, B.

Given the bootstrap estimators \hat{h}_{x}^{(j)}(t), j=1,\dots,B, the estimator of the original data set \hat{h}_{x}(t) and the quantiles k_{\alpha/2}^s, k_{1-\alpha/2}^s and \bar k_{1-\alpha}^s, the studentized (double bootstrap) confidence interval (CI) for \hat{h}_{x}(t) is given by

\Bigg ( \hat{h}_{x}(t) - \hat \sigma k_{1-\alpha/2}^s, \hat{h}_{x}(t) - \hat \sigma k_{\alpha/2}^s \Bigg ).

The symmetric studentized confidence interval (CI) for \hat{h}_{x}(t) defined by

\Bigg ( \hat{h}_{x}(t) - \hat \sigma \bar k_{1-\alpha}^s, \hat{h}_{x}(t) + \hat \sigma\bar k_{1-\alpha}^s \Bigg ).

For the confidence intervals for the logarithm of the hazard rate function first set k_{\alpha/2}^{L,s}, k_{1-\alpha/2}^{L,s} and \bar k_{1-\alpha}^{L,s} be the \alpha/2, 1-\alpha/2 and 1-\alpha quantile of |\hat{L}_{x}^{(j)}(t) - \hat{L}_{x}(t)|/\hat \sigma_j, j=1,\dots, B.

The studentized (double bootstrap) CI confidence interval for the logarithm of the hazard rate function is

\Bigg ( \hat{L}_{x}(t) - \hat \sigma k_{1-\alpha/2}^{L,s}, \hat{L}_{x}(t) - \hat \sigma k_{\alpha/2}^{L,s} \Bigg ).

Its symmetric studentized (double bootstrap) CI for the log hazard is

\Bigg ( \hat{L}_{x}(t) - \hat \sigma \bar k_{1-\alpha}^{L,s}, \hat{L}_{x}(t) + \hat \sigma \bar k_{1-\alpha}^{L,s} \Bigg ).

These confidence intervals are transformed back to yield the following confidence intervals for the hazard rate function h_x(t):

\Bigg ( \hat{h}_{x}(t) e^{- \hat \sigma k_{1-\alpha/2}^{L,s}}, \hat{h}_{x}(t) e^{- \hat \sigma k_{\alpha/2}^{L,s}} \Bigg ),

and the symmetric version

\Bigg ( \hat{h}_{x}(t) e^{- \hat \sigma \bar k_{1-\alpha}^{L,s}}, \hat{h}_{x}(t) e^{ \hat \sigma \bar k_{1-\alpha}^{L,s}} \Bigg ).

Note: The bootstrap matrix studentized.h.hat.star is assumed to contain estimates produced using the same time grid as t.grid and the same estimator used to generate h.hat.x.

Value

A data frame with the following columns:

See Also

bootstrap.h.hat.star, h.hat.star.x.theta

Examples

 
X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'

data <- pbc2[,c(id.name, X.1.name, X.2.name, T.name, s.name, delta.name)]
n.obs <- length(data$id)
n<-max(  as.double(data[,'id']) )

data.id <- to.id(data)

theta.1  <- 0.0702 #0.149
theta.2 <- 0.0856 #0.10
x.1 = 0 # refers to zero mean variables - slightly high
x.2 = 1.9 # refers to zero mean variable - high
b = 0.42#theta.alb * b.alb + theta.bil *b.bil # 7
x = theta.1 * x.1 + theta.2 *x.2
n.s.grid<-50

X.1.obs=data[,X.1.name] -mean(data[, X.1.name])
X.1.i.T.i=data.id[,X.1.name] -mean(data.id[, X.1.name])
X.2.obs=data[,X.2.name]  -mean(data[, X.2.name])
X.2.i.T.i=data.id[,X.2.name] -mean(data.id[, X.2.name])

X.obs.list=list(X.1.obs, X.2.obs)
X.i.T.i.list=list(X.1.i.T.i, X.2.i.T.i)

# Calculate the indexed HQM estimator on the original sample:
h.hat.x.theta.output<- get.h.hat.x.theta(pbc2, id.name, n.s.grid,  X.obs.list, X.i.T.i.list, 
                      T.name = 'years', s.name = 'year',  delta.name = 'status2', 
                      theta= c(theta.1,  theta.2), b, x)
                          
h.hat.x<-h.hat.x.theta.output[,2] # Indexed HQM estimator on original sample
t.grid<-h.hat.x.theta.output[,1] # evaluation grid points
n.grid<- n.s.grid # length(h.hat.x)



#  Create bootstrap samples by group 
set.seed(1)  
B<-10 # 20 # 5 for display purposes only; for sensible results use B=200 (slower) 
B1<- 10 # 20  # 5 for display purposes only; use B1=20 (slower)
bootstrap.samples<-list()
for(j in 1:B)
{
  rows.star<-c()
  ids.star<-c()
  ids.drawn <- sample (n, replace = TRUE)  
  for(l in 1:n)
  {
    rows.l<-which(data[,id.name]==ids.drawn[l])
    rows.star<-c(rows.star, rows.l)
    ids.l<-rep(ids.drawn[l], times=length(rows.l))
    ids.star<-c(ids.star, ids.l)
  }
  data.star.unsorted<-data[rows.star,]
  data.star.unsorted<-data[rows.star,]
  bootstrap.samples[[j]]<- data.star.unsorted[order(data.star.unsorted$id),]  
}
 


# Simulate true hazard rate function:
h.true<- h.true.simulated(bootstrap.samples, id.name='id', n.grid, 
              X.1.name=X.1.name, X.2.name= X.2.name, 
              T.name = T.name, s.name = s.name, 
              delta.name = delta.name, theta = c(theta.1,  theta.2), b)

# Bootstrap the original indexed HQM estimator:                              
double.bootstrap.output<-double.bootstrap.h.hat.star(B, B1, bootstrap.samples, X.1.name, X.2.name, 
             T.name, s.name, delta.name, b, x, h.true=h.true, 
             theta.start=c(theta.1,  theta.2), h.hat.x=h.hat.x, id.name= 'id', data=data)
             
# Construct Ci's:  
alpha.level<-0.05 
CIs.studentized<-bootstrap.CI.studentized(n.grid, double.bootstrap.output, t.grid, h.hat.x, 
                                          alpha.level)

# extract Plain + symmetric CIs
CI.upper<-CIs.studentized[,"CI.upper"]
CI.lower<-CIs.studentized[,"CI.lower"]
CI.sym.upper<-CIs.studentized[,"CI.sym.upper"]
CI.sym.lower<-CIs.studentized[,"CI.sym.lower"]

#Plot the selected CIs
J<-80 #select the first 80 grid points (for display purposes only)
plot(t.grid[1:J], h.hat.x[1:J], type="l", ylim=c(0,2), ylab="Hazard rate", xlab="time", 
            lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.upper[1:J],  rev(CI.lower[1:J])), 
      col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.sym.lower[1:J], lty=2, lwd=2)   

# extract transformed from Log HR + symmetric CIs
CI.L.upper<-CIs.studentized[,"CI.L.upper"]
CI.L.lower<-CIs.studentized[,"CI.L.lower"]
CI.L.sym.upper<-CIs.studentized[,"CI.L.sym.upper"]
CI.L.sym.lower<-CIs.studentized[,"CI.L.sym.lower"]

#Plot the selected CI's
plot(t.grid[1:J], log(h.hat.x[1:J]), type="l", ylim=c(-5,4), ylab="Log Hazard rate", 
      xlab="time",  lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.L.upper[1:J], rev(CI.L.lower[1:J])), 
      col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.L.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.L.sym.lower[1:J], lty=2, lwd=2)  

# extract Log HR + symmetric CIs  
CI.exp.upper<-CIs.studentized[,"CI.exp.upper"]
CI.exp.lower<-CIs.studentized[,"CI.exp.lower"]
CI.exp.sym.upper<-CIs.studentized[,"CI.exp.sym.upper"]
CI.exp.sym.lower<-CIs.studentized[,"CI.exp.sym.lower"]

#Plot the selected CIs
plot(t.grid[1:J],   h.hat.x[1:J] , type="l", ylim=c(0,2), ylab="Hazard rate", xlab="time", 
      lwd=2)
polygon(x = c(t.grid[1:J], rev(t.grid[1:J])), y = c(CI.exp.upper[1:J], rev(CI.exp.lower[1:J])), 
    col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(t.grid[1:J], CI.exp.sym.upper[1:J], lty=2, lwd=2 ) 
lines(t.grid[1:J],  CI.exp.sym.lower[1:J], lty=2, lwd=2)  

Bootstrap Estimation of Hazard Function and Index Parameters

Description

Performs bootstrap estimation of hazard rates and corresponding index parameters for a given set of bootstrap samples. For each bootstrap replicate, the function re-estimates the index parameters via optimisation and computes hazard estimates using h.hat.star.x.theta. The output is used as input in functions bootstrap.CI.quantile and bootstrap.CI.pivotal

Usage

bootstrap.h.hat.star(B, bootstrap.samples, X.1.name, X.2.name,T.name, 
                             s.name, delta.name, b, x, h.true, theta.start, n.grid)

Arguments

Details

For each bootstrap iteration k = 1, \dots, B, the function:

1. Extracts the bootstrap sample data.star.

2. Computes centred marker values at the subject and observation level.

3. Estimates index parameters by minimising Q.theta using optim.

4. Computes the bootstrap hazard estimate via h.hat.star.x.theta.

The outputs are matrices collecting the hazard estimates and estimated index parameter vectors across bootstrap replicates.

Value

A matrix of dimension n.grid × B containing the bootstrap hazard estimates.

See Also

h.hat.star.x.theta, Q.theta, to.id

Examples

X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'

theta.1  <- 0.0702  
theta.2 <- 0.0856  
x.1 = 0 # refers to zero mean variables - slightly high
x.2 = 1.9 # refers to zero mean variable - high
b = 0.42 
x = theta.1 * x.1 + theta.2 *x.2

# first simulate true HR function:
data <- pbc2[,c(id.name, X.1.name, X.2.name, T.name, s.name, delta.name)]
n.obs <- length(data$id)
n<-max(  as.double(data[,'id']) )
data.id <- to.id(data)


#  Create bootstrap samples by group:
set.seed(1)  
B<- 10 # 400 #50
bootstrap.samples<-list()
for(j in 1:B)
{
  rows.star<-c()
  ids.star<-c()
  ids.drawn <- sample (n, replace = TRUE)  
  for(l in 1:n)
  {
    rows.l<-which(data[,id.name]==ids.drawn[l])
    rows.star<-c(rows.star, rows.l)
    ids.l<-rep(ids.drawn[l], times=length(rows.l))
    ids.star<-c(ids.star, ids.l)
  }
  data.star.unsorted<-data[rows.star,]
  data.star.unsorted<-data[rows.star,]
  bootstrap.samples[[j]]<- data.star.unsorted[order(data.star.unsorted$id),] #data[rows.star,]
}
h.true<- h.true.simulated(bootstrap.samples, id.name='id', 100,  X.1.name=X.1.name,  
                    X.2.name= X.2.name, T.name = T.name, 
                    s.name = s.name,  delta.name = delta.name, 
                    theta = c(theta.1,  theta.2), b)
                              
h.hat.star.matrix <- bootstrap.h.hat.star( B,  bootstrap.samples, X.1.name, X.2.name, 
          T.name,  s.name,  delta.name ,  b = 0.4, x = 1.0,  
          h.true = h.true, theta.start = c(0.07, 0.08), n.grid = 100)
          
#return bootstrap hazard rate estimators in marix format:
h.hat.star.matrix

D matrix entries, used for the implementation of the local linear kernel

Description

Calculates the entries of the D matrix in the definition of the local linear kernel

Usage

c2.u.integrand(b,u,v, K)

Arguments

Details

Implements the caclulation of all d \times d entries of matrix D, which is part of the definition of the local linear kernel. The actual calculation is performed by

d_{jk} = \sum_{i=1}^n \int_0^T K_b(x-X_i(s))\{x-X_{ij}(s)\}\{x-X_{ik}(s)\} Z_i(s)ds,

Value

scalar value, the result of d_{jk}.

Local linear weight functions

Description

Implements the weights to be used in the local linear HQM estimator.

Usage

c0.u(w, u, b, K)
c1.u(w, u, b, K)
c2.u(w, u, b, K)

Arguments

Details

The function implements the local linear weights in the definition of the estimator \hat{h}_x(t), see also h.hat.x.t

Value

A vector of s_n(x) for all values x on the marker grid.

Confidence bands

Description

Implements the uniform and pointwise confidence bands for the future conditional hazard rate based on the last observed marker measure.

Usage

confidence.bands(data,   X.name, T.name = 'years',
            s.name = 'year', delta.name = 'status2', x, b)

Arguments

Details

The function confidence.bands implements the pointwise and uniform confidence bands for the estimator of the future conditional hazard rate \hat h_x(t). The confidence bands are based on a wild bootstrap approach {h^*}_{{x_*},B}(t).

Pointwise: For a given t\in (0,T) generate {h^*}_{{x_*},B}^{(1)}(t),...,{h^*}_{{x_*},B}^{(N)}(t) for N = 1000 and order it {h^*}_{{x_*},B}^{[1]}(t)\leq ...\leq {h^*}_{{x_*},B}^{[N]}(t). Then

\hat{I}^1_{n,N} = \Bigg[\hat{h}_{x_*}(t) - \hat{\sigma}_{{G}_{x_*}}(t)\frac{{h^*}_{{x_*},B}^{[ N(1-\frac{\alpha}{2})]}(t)}{\sqrt{n}}, \hat{h}_{x_*}(t) - \hat{\sigma}_{ {G}_x}(t)\frac{{h^*}_{{x_*},B}^{[ N\frac{\alpha}{2}]}(t)}{\sqrt{n}}\Bigg]

is a 1-\alpha pointwise confidence band for h_{x_*}(t), where \hat{\sigma}_{{G}_{x_*}}(t) is a bootrap estimate of the variance. For more details on the wild bootstrap approach, please see wild.bootstrap.summands and g.hat.t.x.

Uniform: Generate \bar{h}_{{x_*},B}^{(1)}(t),...,\bar{h}_{{x_*},B}^{(N)}(t) for N = 1000 for all t\in [\delta_T,T-\delta_T] and define W^{(i)} = \sup_{t\in[0,T]}\big|\bar{h}_{{x_*},B}^{(i)}(t)| for i = 1,...,N. Order W^{[1]} \leq ... \leq W^{[N]}. Then

\hat{I}^2_{n,N} = \Bigg[\hat{h}_{x_*}(t) \pm \hat{\sigma}_{{G}_{x_*}}(t) \frac{W^{[ N(1 - \alpha)]}}{\sqrt{n}} \Bigg]

is a 1-\alpha uniform confidence band for h_{x_*}(t).

Value

A list with pointwise, uniform confidence bands and the estimator \hat h_x(t) for all possible time points t:

See Also

g.hat.t.x, wild.bootstrap.summands

Examples

b = 10
x = 3
n.s.grid<-100
s.obs = pbc2$year
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))


bands = confidence.bands(pbc2,   'serBilir', T.name = 'years',
                    s.name = 'year', delta.name = 'status2', x, b)

J = 80
plot(s.grid[1:J], bands$h.hat.x.t[1:J], type = "l", ylim = c(0,1), ylab = 'Hazard', xlab = 'Years')

lines(s.grid[1:J], bands$I.hat.1.upper[1:J], col = "red")
lines(s.grid[1:J], bands$I.hat.1.lower[1:J], col = "red")
lines(s.grid[1:J], bands$I.hat.2.upper[1:J], col = "blue")
lines(s.grid[1:J], bands$I.hat.2.lower[1:J], col = "blue")

Split dataset for K-fold cross validation

Description

Creates multiple splits of a dataset which is then used in the bandwidth selection with K-fold cross validation.

Usage

data.minus.I.j(size.I.j, data)

Arguments

Details

The function data.minus.I.j takes a data frame and transforms it into K = n/I data frames with I individuals missing from each data frame. Let I_j be sets of indices with \cup_{j=1}^K I_j = \{1,...,n\}, I_k\cap I_j = \emptyset and |I_j| = |I_k| = I for all j, k \in \{1,...,K\}. Then data frames with \{1,...,n \}/I_j individuals are created.

Value

A list of data frames with size.I.j individuals missing in the above way.

See Also

CV.b

Examples

data.minus.I.j.list = data.minus.I.j(26, pbc2)

Bootstrap Estimation of Hazard Function and Index Parameters

Description

Performs bootstrap estimation of hazard rates and their standard deviation at each grid point (double boostrap) together with the corresponding index parameters for a given set of bootstrap samples. The output of the function is used as input in bootstrap.CI.studentized.

Usage

double.bootstrap.h.hat.star(B, B1, bootstrap.samples, X.1.name, X.2.name, 
                            T.name, s.name, delta.name, b, x, h.true, 
                            theta.start,  h.hat.x, id.name, data)

Arguments

Details

For each bootstrap iteration k = 1, \dots, B, the function:

1. Extracts the bootstrap sample data.star.

2. Computes centred marker values at the subject and observation level.

3. Estimates index parameters by minimising Q.theta using optim.

4. Computes the bootstrap hazard estimate via h.hat.star.x.theta.

The outputs are matrices collecting the hazard estimates and estimated index parameter vectors across bootstrap replicates.

Value

A list of matrices of dimension n.grid × B containing the bootstrap hazard estimates, the logarithm of the hazard rate estimates and and two vectors of the estimate's standard deviations at each grid point.

See Also

h.hat.star.x.theta, Q.theta, to.id

Examples

  
# See the example for  function: bootstrap.CI.studentized 
 
 

Preparations for bandwidth selection

Description

Calculates an intermediate part for the K-fold cross validation.

Usage

g.hat.i(h.hat.matrix, X.pos, n.s.grid, n, Z.i.ds)

Arguments

Details

The function g.hat.i calculates a key component for the bandwidth selection

\hat{g}^{-I_j}_i(t) = \int_0^t Z_i(s) \hat{h}^{-I_j}_{X_i(s)}(t-s) ds,

where \hat{h}^{-I_j} is estimated without information from all counting processes i with i \in I_j and Z is the exposure.

Value

A matrix with \hat{g}^{-I_j}_i(t) for all individuals i and time grid points t.

See Also

CV.b

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN <- make.dN(pbc2, data.id, X.grid=X.grid, s.grid=s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid, n.X.grid,
            s.pos, X.pos, 'years', n)

b = 1.7
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos,
              X.pos,'years', n)
dN.i  <- make.dN.i(s.grid=s.grid, n.s.grid, T.i, delta.i, n)

t = 2

h.hat.matrix = t(sapply(s.grid[1:99], function(t.point){
          h.hat.x.t.vec(X.grid, s.grid, n.s.grid, alpha.hat, t, b, Z.i.ds, X.pos, n)}))

g.hat.i(h.hat.matrix, X.pos, n.X.grid, n, Z.i.ds)

Computation of a key component for wild bootstrap

Description

Implements a key part for the wild bootstrap of the hqm estimator.

Usage

g.hat.t.x(X.grid, s.grid, n.s.grid, X.pos, x, t, b, Z.i.ds, Z.ds, n)

Arguments

Details

The function implements

\hat{g}_{t,x}(z) = \frac{1}{n} \sum_{ j = 1}^n \int^{T-t}_0 \hat{E}(X_j(t+s))^{-1} K_b(z,X_j(t+s)) Z_j(t+s)Z_j(s)K_b(x,X_j(s))ds,

for every value z on the marker grid, where \hat{E}(x) = \frac{1}{n} \sum_{j=1}^n \int_0^T K_b(x,X_j(s))Z_j(s)ds, Z the exposure and X the marker.

Value

A vector of \hat{g}_{t,x}(z) for all values z on the marker grid.

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
X.obs = pbc2$serBilir
s.obs = pbc2$year
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

Z.i.ds<-make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos, X.pos,
            'years', n)
Z.ds<-make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid,n.s.grid,n.X.grid, s.pos,X.pos,
            'years', n)

t = 2
x = 2
b = 10

g.hat.t.x(X.grid, s.grid, n.s.grid, X.pos, x, t, b, Z.i.ds, Z.ds, n)

Marker-only hazard rate

Description

Calculates the marker-only hazard rate for time dependent data.

Usage

get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Arguments

Details

The function get.alpha.hat implements the marker-only hazard estimator

\hat{\alpha}_i(z) = \frac{\sum_{k\neq i}\int_0^T K_{b_1}(z-X_k(s))\mathrm{d}N_k(s)}{\sum_{k\neq i}\int_0^T K_{b_1}(z-X_k(s))Z_k(s)\mathrm{d}s},

where X is the marker and Z is the exposure. The marker-only hazard is defined as the underlying hazard which is not dependent on time

\alpha(X(t),t) = \alpha(X(t)).

Value

A vector of marker-only values for X.grid.

See Also

h.hat.x.t

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN  <- make.dN(pbc2, data.id, X.grid=X.grid, s.grid=s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid,
            n.X.grid, s.pos, X.pos, T.name = 'years', n)

b = 1.7
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Local constant future conditional hazard rate estimator

Description

Calculates the (indexed) local constant future hazard rate function, conditional on a marker value x, across across a set of time values t.

Usage

get.h.hat.x(data, X.name, s.grid, T.name, s.name, delta.name, x, b)

Arguments

Details

The function get.h.hat.x implements the indexed local linear future conditional hazard estimator

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(\hat \theta^T X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-\hat \theta^T X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- \hat \theta^T X_i(s))\mathrm {d}s},

across a grid of possible time values t, where, for a positive integer p, \hat \theta^T = (\hat \theta_1, \dots, \hat \theta_p ) is the vector of the estimated indexing parameters, X_i = (X_{1, i}, \dots, X_{i,p}) is a vector of markers for indexing, Z_i is the exposure and \alpha(z) is the marker-only hazard, see get.alpha.hat for more details and K_b = b^{-1}K(./b) is an ordinary kernel function, e.g. the Epanechnikov kernel. For p=1 and \hat \theta = 1 the above estimator becomes the HQM hazard rate estimator conditional on one covariate,

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i( X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s},

defined in equation (2) of Bagkavos et al. (2025), doi:10.1093/biomet/asaf008.

Value

A vector of \hat h_x(t) for a grid of possible time values t.

References

Bagkavos, I., Isakson, R., Mammen, E., Nielsen, J., and Proust–Lima, C. (2025). Biometrika, 112(2), asaf008. doi:10.1093/biomet/asaf008

See Also

get.alpha.hat, h.hat.x.t, get.h.tilde.x

Examples

library(survival)
b = 10
x = 3
Landmark <- 2
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]

n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))

b=0.9

h.tilde.x<-get.h.tilde.x(pbcT1,'albumin', s.grid = s.grid, T.name='years',
                  s.name='year',delta.name='status2',2,0.9)
h.hat.x<-get.h.hat.x(pbcT1,'albumin', s.grid = s.grid, T.name='years',
                  s.name='year',delta.name='status2',2,0.9)

#Caclulate the local contant and local linear survival functions
s.grid.2  = seq(Landmark, 14,  length= n.grid-1 )
S.tilde.alb<- make.S.hat(    (s.grid.2[2]-s.grid.2[1])/4 , h.tilde.x)
S.hat.alb<- make.S.hat(    (s.grid.2[2]-s.grid.2[1])/4 , h.hat.x)

#For comparison, also calculate the Kaplan-Meier
kma2<- survfit(Surv(years , status2) ~ 1, data = pbcT1)

#Plot the survival functions:
plot(s.grid.2, S.tilde.alb,  type="l", col=1, lwd=2, ylab="Survival probability",
     xlab="Marker level")
lines(s.grid.2, S.hat.alb,  lty=2, lwd=2, col=2)
lines(kma2$time, kma2$surv, type="s",  lty=2, lwd=2, col=3)

legend("topright", c(  "Local linear HQM", "Local constant HQM",
        "Kaplan-Meier"), lty=c(1, 2, 2), col=1:3, lwd=2, cex=1.7)


## Not run: 
#Example of get.h.hat.x with two indexed covariates:
#First, estimate the joint model for Albumin and Bilirubin combined:
library(JM)
lmeFit <- lme(albumin + serBilir~ year, random = ~ year | id, data = pbc2)
coxFit <- coxph(Surv(years, status2) ~ albumin + serBilir, data = pbc2.id, x = TRUE)
jointFit <- jointModel(lmeFit, coxFit, timeVar = "year", method = "piecewise-PH-aGH")

Landmark <- 1
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]


n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))


# Index Albumin and Bilirubin:
x.alb = 3  # slightly low albumin value
x.bil = 1.9  # slightly  high bilirubin value

theta.alb  <- 0.0702
theta.bil <- 0.0856
theta.T.X.obs = theta.alb * pbcT1$albumin + theta.bil *pbcT1$serBilir 
 # theta.T.X.obs is now the indexed marker
x = theta.alb * x.alb + theta.bil *x.bil #conditioning value
## store the combine variable in place of 'drug' column which is redundant
pbcT1$drug<- theta.T.X.obs 
## i.e. 'drug' corresponds to indexed bilirubin and albumin

timesS2 <- seq(Landmark,14,by=0.5)
predT1 <- survfitJM(jointFit, newdata = pbcT1, survTimes = timesS2)
nm<-length(predT1$summaries)

mat.out1<-matrix(nrow=length(timesS2), ncol=nm)
for(r in 1:nm)
{
  SurvLand <- predT1$summaries[[r]][,"Mean"][1] #obtain mean predictions
  mat.out1[,r] <- predT1$summaries[[r]][,"Mean"]/SurvLand
}
JM.surv.est<-rowMeans(mat.out1, na.rm=TRUE) #average the resulting JM estimates

# calculate indexed local linear HQM estimator for bilirubin and albumin
b.alb = 1.5
b.bil = 4
b.hqm   =  theta.alb * b.alb + theta.bil *b.bil # bandwidth for HQM estimator
h.hat.x<- get.h.hat.x(pbcT1, 'drug', s.grid  =s.grid, T.name = 'years',  s.name = 'year',
                 delta.name = 'status2', x, b.hqm)
s.grid.2  = seq(Landmark, 14,  length=n.grid-1) #grid points for HMQ estimator
S.hat.x<- make.S.hat((s.grid.2[2]-s.grid.2[1])/5, h.hat.x) # transform HR to Survival func.

km.land<- survfit(Surv(years , status2) ~ 1, data = pbcT1) #KM estimate

#Plot the survival functions:
plot(s.grid.2, S.hat.x, type="l", ylim=c(0,1), xlim=c(Landmark,14),
            ylab="Survival probability", xlab="years",lwd=2)
lines(timesS2, JM.surv.est, col=2,  lwd=2, lty=2)
lines(km.land$time, km.land$surv, type="s",lty=2, lwd=2, col=3)
legend("bottomleft", c("HQM est.", "Joint Model est.", "Kaplan-Meier"),
      lty=c(1,2,2),  col=1:3, lwd=2, cex=1.7)


## End(Not run)

Local linear future conditional hazard estimator (wrapper)

Description

Compute the indexed local linear future conditional hazard rate estimator on a grid.

Usage

get.h.hat.x.theta(data, id.name, n.grid,  X.obs.list, X.i.T.i.list, T.name = 'years', 
                            s.name = 'year', delta.name = 'status2', theta, b, x)

Arguments

Value

A data frame with columns t and h.hat.x.t.

Examples

b.alb = 0.9   
b.bil = 4

x.alb = 1 # refers to zero mean variables - slightly high
x.bil = 1.9 # refers to zero mean variable - high

theta.alb  <- 0.0702 #0.149
theta.bil <- 0.0856 #0.10


b = 0.42  
x = x.alb * theta.alb + x.bil *theta.bil

X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'
n.s.grid<-50

data.use<-pbc2
data.use.id<-to.id(data.use)
#data.use.id<-data.use.id[complete.cases(data.use.id), ]

# mean adjust the data:
X.1.obs=data.use[,X.1.name] -mean(data.use[, X.1.name])
X.1.i.T.i=data.use.id[,X.1.name] -mean(data.use.id[, X.1.name])
X.2.obs=data.use[,X.2.name]  -mean(data.use[, X.2.name])
X.2.i.T.i=data.use.id[,X.2.name] -mean(data.use.id[, X.2.name])

X.obs.list=list(X.1.obs, X.2.obs)
X.i.T.i.list=list(X.1.i.T.i, X.2.i.T.i)

h.hat.x.theta.output<-get.h.hat.x.theta(pbc2, id.name, n.s.grid,  X.obs.list, 
                                        X.i.T.i.list, T.name = 'years', 
    s.name = 'year',  delta.name = 'status2', theta= c(theta.alb,  theta.bil), b, x)
h.hat.x<-h.hat.x.theta.output[,2]
t.grid<-h.hat.x.theta.output[,1]
n.grid<-length(h.hat.x)

Local linear future conditional hazard rate estimator

Description

Calculates the (indexed) local linear future hazard rate function, conditional on a marker value x, across a set of time values t.

Usage

 get.h.tilde.x(data, X.name, s.grid, T.name, s.name, delta.name, x, b)

Arguments

Details

The function get.h.tilde.x uses a local linear kernel to implement the indexed local linear future conditional hazard estimator

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(\hat \theta^T X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-\hat \theta^T X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- \hat \theta^T X_i(s))\mathrm {d}s},

across a grid of possible time values t, where, for a positive integer p, \hat \theta^T = (\hat \theta_1, \dots, \hat \theta_p ) is the vector of the estimated indexing parameters, X_i = (X_{1, i}, \dots, X_{i,p}) is a vector of markers for indexing, Z_i is the exposure and \alpha(z) is the marker-only hazard, see get.alpha.hat for more details. For p=1 and \hat \theta = 1 the above estimator becomes the HQM hazard rate estimator conditional on one covariate,

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i( X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s},

defined in equation (2) of Bagkavos et al. (2025) doi:10.1093/biomet/asaf008. In the place of K_b(), get.h.tilde.x uses the kernel

K_{x,b}(u)= \frac{K_b(u)-K_b(u)u^T D^{-1}c_1}{c_0 - c_1^T D^{-1} c_1},

where K_b() = b^{-1}K(./b) with K being an ordinary kernel, e.g. the Epanechnikov kernel, c_1 = (c_{11}, \dots, c_{1d})^T, D = (d_{ij})_{(d+1) \times (d+1)} with

c_0 = \sum_{i=1}^n \int_0^T K_b(x-\hat \theta^T X_i(s)) Z_i(s)ds,

c_{ij} = \sum_{i=1}^n \int_0^T K_b(x-\hat \theta^T X_i(s))\{x-\hat \theta^T X_{ij}(s)\} Z_i(s)ds,

d_{jk} = \sum_{i=1}^n \int_0^T K_b(x-\hat \theta^T X_i(s))\{x-\hat \theta^T X_{ij}(s)\}\{x-\hat \theta^T X_{ik}(s)\} Z_i(s)ds,

see also Nielsen (1998), doi:10.1080/03461238.1998.10413997.

Value

A vector of \hat h_x(t) for a grid of possible time values t.

References

Bagkavos, I., Isakson, R., Mammen, E., Nielsen, J., and Proust–Lima, C. (2025). Biometrika, 112(2), asaf008. doi:10.1093/biomet/asaf008

Nielsen (1998), Marker dependent kernel hazard estimation from local linear estimation, Scandinavian Actuarial Journal, pp. 113-124. doi:10.1080/03461238.1998.10413997

See Also

get.alpha.hat, h.hat.x.t, get.h.hat.x

Examples

library(survival)
library(JM)

# Compare Local constant and local linear estimator for a single covariate,
# use KM for reference.
# Albumin marker, use landmarking
Landmark <- 2
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]
b=0.9

n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))


h.tilde.x<-get.h.tilde.x(pbcT1, 'albumin', s.grid = s.grid, T.name = 'years', s.name = 'year',
                  delta.name = 'status2', 2, 0.9)
h.hat.x<-get.h.hat.x(pbcT1, 'albumin', s.grid = s.grid, T.name = 'years', s.name = 'year',
                delta.name = 'status2', 2, 0.9)

#Calculate the local contant and local linear survival functions
s.grid.2  = seq(Landmark, 14,  length=n.grid-1)
S.tilde.alb<- make.S.hat((s.grid.2[2]-s.grid.2[1])/4 , h.tilde.x)
S.hat.alb<- make.S.hat((s.grid.2[2]-s.grid.2[1])/4 , h.hat.x)

#For comparison, also calculate the Kaplan-Meier
kma2<- survfit(Surv(years , status2) ~ 1, data = pbcT1)

#Plot the survival functions:
plot(s.grid.2, S.tilde.alb,  type="l", col=1, lwd=2, ylab="Survival probability",
                                                        xlab="Marker level")
lines(s.grid.2, S.hat.alb,  lty=2, lwd=2, col=2)
lines(kma2$time, kma2$surv, type="s",  lty=2, lwd=2, col=3)
legend("topright", c(  "Local linear HQM", "Local constant HQM", "Kaplan-Meier"),
        lty=c(1, 2, 2), col=1:3, lwd=2, cex=1.7)


## Not run: 
#Example of get.h.tilde.x with a single covariate (no indexing):
#Compare JM, HQM and KM for Bilirubin
b = 10
Landmark <- 1
lmeFit <- lme(serBilir ~ year, random = ~ year | id, data = pbc2)
coxFit <- coxph(Surv(years, status2) ~ serBilir, data = pbc2.id, x = TRUE)

jointFit0 <- jointModel(lmeFit, coxFit, timeVar = "year",
                                method = "piecewise-PH-aGH")
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]
n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))

timesS1 <- seq(1,14,by=0.5)
predT1 <- survfitJM(jointFit0, newdata = pbcT1,survTimes = timesS1)
nm<-length(predT1$summaries)

mat.out1<-matrix(nrow=length(timesS1), ncol=nm)
for(r in 1:nm)
{
  SurvLand <- predT1$summaries[[r]][,"Mean"][1]
  mat.out1[,r] <- predT1$summaries[[r]][,"Mean"]/SurvLand
}
sfit1y<-rowMeans(mat.out1, na.rm=TRUE)

h.tilde.x<- get.h.tilde.x(pbcT1, 'serBilir', s.grid= s.grid,  T.name = 'years',
        s.name = 'year', delta.name = 'status2', 1, 10)
s.grid.1  = seq(Landmark, 14,  length=n.grid-1)
S.tilde.bil<- make.S.hat((s.grid.1[2]-s.grid.1[1])/5.4 , h.tilde.x)
kma1<- survfit(Surv(years , status2) ~ 1, data = pbcT1)

plot(s.grid.1, S.tilde.bil, type="l", ylim=c(0,1), xlim=c(Landmark,14),
                    ylab="Survival probability", xlab="years",lwd=2)
lines(timesS1, sfit1y, col=2,  lwd=2, lty=2)
lines(kma1$time, kma1$surv, type="s", lty=2, lwd=2, col=3 )
legend("bottomleft", c("HQM est.", "Joint Model est.", "Kaplan-Meier"),
                                    lty=c(1,2,2), col=1:3, lwd=2, cex=1.7)

#Example of get.h.tilde.x with two indexed covariates:

#First, estimate the joint model for Albumin and Bilirubin combined:
lmeFit <- lme(albumin + serBilir~ year, random = ~ year | id, data = pbc2)
coxFit <- coxph(Surv(years, status2) ~ albumin + serBilir, data = pbc2.id,
                                                                  x = TRUE)
jointFit <- jointModel(lmeFit, coxFit, timeVar = "year",
                                               method = "piecewise-PH-aGH")

Landmark <- 1
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]
n.grid<-50 # 50 grid points to evaluate the estimates
s.obs<- pbcT1[, 'year']
s.grid <-  seq(0, max(s.obs), max(s.obs)/( n.grid-1))


# Index Albumin and Bilirubin:
x.alb = 3  # slightly low albumin value
x.bil = 1.9  # slightly  high bilirubin value

theta.alb  <- 0.0702
theta.bil <- 0.0856
# theta.T.X.obs is now the indexed marker
theta.T.X.obs = theta.alb * pbcT1$albumin + theta.bil *pbcT1$serBilir 
x = theta.alb * x.alb + theta.bil *x.bil #conditioning value

pbcT1$drug<- theta.T.X.obs ## store theta.T.X.obs in place of 'drug' column which is redundant here
## i.e. 'drug' corresponds to indexed bilirubin and albumin

timesS2 <- seq(Landmark,14,by=0.5)
predT1 <- survfitJM(jointFit, newdata = pbcT1,survTimes = timesS2)
nm<-length(predT1$summaries)

mat.out1<-matrix(nrow=length(timesS2), ncol=nm)
for(r in 1:nm)
{
  SurvLand <- predT1$summaries[[r]][,"Mean"][1] #obtain mean predictions
  mat.out1[,r] <- predT1$summaries[[r]][,"Mean"]/SurvLand
}
JM.surv.est<-rowMeans(mat.out1, na.rm=TRUE) #average the resulting JM estimates

# calculate indexed local linear HQM estimator for bilirubin and albumin
b.alb = 1.5
b.bil = 4
b.hqm   =  theta.alb * b.alb + theta.bil *b.bil # bandwidth for HQM estimator
h.tilde.x<- get.h.tilde.x(pbcT1, 'drug', s.grid = s.grid, T.name = 'years',  s.name = 'year',
                                               delta.name = 'status2', x, b.hqm)
s.grid.2  = seq(Landmark, 14,  length=n.grid-1) #grid points for HMQ estimator
S.tilde.x<- make.S.hat((s.grid.2[2]-s.grid.2[1])/5  ,h.tilde.x) # transform HR to Survival func.

km.land<- survfit(Surv(years , status2) ~ 1, data = pbcT1) #KM estimate

#Plot the survival functions:
plot(s.grid.2, S.tilde.x, type="l", ylim=c(0,1), xlim=c(Landmark,14),
    ylab="Survival probability", xlab="years",lwd=2)
lines(timesS2, JM.surv.est, col=2,  lwd=2, lty=2)
lines(km.land$time, km.land$surv, type="s",lty=2, lwd=2, col=3)
legend("bottomleft", c("HQM est.", "Joint Model est.", "Kaplan-Meier"),
        lty=c(1,2,2),  col=1:3, lwd=2, cex=1.7)

## End(Not run)

Prepare for Cross validation bandwidth selection

Description

Implements the calculation of the hqm estimator on cross validation data sets. This is a preparation for the cross validation bandwidth selection technique for future conditional hazard rate estimation based on marker information data.

Usage

h.hat.minus.I.j(data, data.id, X.name, T.name = 'years',
        s.name = 'year',delta.name = 'status2', n, size.I.j, b)

Arguments

Details

The function splits the data set via data.minus.I.j and calculates for every splitted data set the hqm estimator

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x-X_i(s))\mathrm {d}s},

for all x on the marker grid and t on the time grid, where X is the marker, Z is the exposure and \alpha(z) is the marker-only hazard, see get.alpha.hat for more details.

Value

A list of matrices for every cross validation data set with \hat{h}_x(t) for all x on the marker grid and t on the time grid.

See Also

CV.b

Examples

data.id = to.id(pbc2)
n = max(as.numeric(pbc2$id))
b = 1.5
size.I.j = 26
h.hat.minus.I.j.list = h.hat.minus.I.j(pbc2, data.id, 'serBilir', 'years', 'year', 'status2', 
                                        n, size.I.j, b)

Prepare for Cross validation index parameter selection

Description

Implements the calculation of the hqm estimator on cross validation data sets. This is a preparation for the cross validation index selection technique for future conditional hazard rate estimation based on marker information data.

Usage

h.hat.minus.I.j.theta(theta, data, data.id, X.1.name, X.2.name, T.name = 'years',
        s.name = 'year',delta.name = 'status2', n, size.I.j, b)

Arguments

Details

The function splits the data set via data.minus.I.j and calculates for every splitted data set the hqm estimator

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(\theta_0^T X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-\theta_0^T X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x-\theta_0^T X_i(s))\mathrm {d}s},

for all x on the marker grid and t on the time grid, where X is the marker, Z is the exposure and \alpha(z) is the marker-only hazard, see get.alpha.hat for more details.

Value

A list of matrices for every cross validation data set with \hat{h}_x(t) for all x on the marker grid and t on the time grid.

See Also

CV.theta

Indexed HQM hazard estimator for one bootstrap sample

Description

Compute the bootstrap estimate of the indexed HQM hazard estimate for a single (bootstrap) sample. The function is meant to be used internally for the calculation of confidence intervals

Usage

h.hat.star.x.theta(theta, data, data.id, n.grid, X.obs.list, X.i.T.i.list,  
                     T.name = 'years', s.name = 'year', delta.name = 'status2', b, x)

Arguments

Value

Numeric vector with estimated hazard on the grid.

Examples

#Single instance of the bootstrap version of the bootstrap version
#of the indexed HQM estimator

b.alb = 0.9   
b.bil = 4

x.alb = 1 # refers to zero mean variables - slightly high
x.bil = 1.9 # refers to zero mean variable - high

theta.alb  <- 0.0702 #0.149
theta.bil <- 0.0856 #0.10

 
b = 0.42 # The result, on the indexed marker 'indmar' of 
         #\code{CV.b(pbc2, 'indmar', 'years', 'year', 'status2', size.I.j=26, seq(0.2,0.4,by=0.01))} 
x = x.alb * theta.alb + x.bil *theta.bil

X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'
n.s.grid<-50

data.use<-pbc2
data.use.id<-to.id(data.use)
data.use.id<-data.use.id[complete.cases(data.use.id), ]

# mean adjust the data:
X.1.obs=data.use[,X.1.name] -mean(data.use[, X.1.name])
X.1.i.T.i=data.use.id[,X.1.name] -mean(data.use.id[, X.1.name])
X.2.obs=data.use[,X.2.name]  -mean(data.use[, X.2.name])
X.2.i.T.i=data.use.id[,X.2.name] -mean(data.use.id[, X.2.name])

X.obs.list=list(X.1.obs, X.2.obs)
X.i.T.i.list=list(X.1.i.T.i, X.2.i.T.i)
h.hat.star<-h.hat.star.x.theta (c(theta.alb,theta.bil), data.use, data.use.id, n.grid=n.s.grid, 
                                X.obs.list, X.i.T.i.list, T.name = 'years', s.name = 'year', 
                                delta.name = 'status2', b, x)

Local constant future conditional hazard rate estimation at a single time point

Description

Calculates the (indexed) future conditional hazard rate for a marker value x and a time value t.

Usage

h.hat.x.t(X.grid, s.grid, X.pos, n.s.grid, alpha.hat, x,t, b, Z.i.ds,n)

Arguments

Details

Function h.hat.x.t implements the future conditional hazard estimator

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(\hat \theta^T X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-\hat \theta^T X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- \hat \theta^T X_i(s))\mathrm {d}s},

where, for a positive integer p, \hat \theta^T = (\hat \theta_1, \dots, \hat \theta_p ) is the vector of the estimated indexing parameters, X_i = (X_{1, i}, \dots, X_{i,p}) is a vector of markers for indexing, Z_i is the exposure and \alpha(z) is the marker-only hazard, see get.alpha.hat for more details and K_b = b^{-1}K(./b) is an ordinary kernel function, e.g. the Epanechnikov kernel. For p=1 and \hat \theta = 1 the above estimator becomes the HQM hazard rate estimator conditional on one covariate,

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x-X_i(s))\mathrm {d}s},

defined in equation (2) of Bagkavos et al. (2025) doi:10.1093/biomet/asaf008.

The future conditional hazard is defined as

h_{x,T}(t) = P\left(T_i\in (t+T, t+T+dt)| \theta_0^T X_i(T)=x, T_i > t+T\right),

where T_i is the survival time and X_i the marker of individual i observed in the time frame [0,T].

Value

A single numeric value of \hat h_x(t).

References

Bagkavos, I., Isakson, R., Mammen, E., Nielsen, J., and Proust–Lima, C. (2025). Biometrika, 112(2), asaf008. doi:10.1093/biomet/asaf008

See Also

get.alpha.hat

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN <- make.dN(pbc2, data.id, X.grid=X.grid, s.grid=s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos, X.pos,
            'years', n)

b = 1.7
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos, X.pos,
              'years', n)

x = 2
t = 2
h.hat = h.hat.x.t(X.grid, s.grid, X.pos, n.s.grid, alpha.hat, x, t, b, Z.i.ds, n)

Hqm estimator on the marker grid

Description

Computes the (indexed) hqm estimator on the marker grid.

Usage

h.hat.x.t.vec(X.grid, s.grid, n.s.grid, alpha.hat, t, b, Z.i.ds, X.pos, n)

Arguments

Details

The function implements the future conditional hazard estimator

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(\hat \theta^T X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-\hat \theta^T X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- \hat \theta^T X_i(s))\mathrm {d}s},

for every x on the marker grid, where, for a positive integer p, \hat \theta^T = (\hat \theta_1, \dots, \hat \theta_p ) is the vector of the estimated indexing parameters, X_i = (X_{1, i}, \dots, X_{i,p}) is a vector of markers for indexing, Z_i is the exposure and \alpha(z) is the marker-only hazard, see get.alpha.hat for more details and K_b = b^{-1}K(./b) is an ordinary kernel function, e.g. the Epanechnikov kernel. For p=1 and \hat \theta = 1 the above estimator becomes the HQM hazard rate estimator conditional on one covariate,

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i( X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s},

defined in equation (2) of Bagkavos et al. (2025), doi:10.1093/biomet/asaf008.

Value

A vector of \hat{h}_{x}(t) for all values x on the marker grid.

References

Bagkavos, I., Isakson, R., Mammen, E., Nielsen, J., and Proust–Lima, C. (2025). Biometrika, 112(2), asaf008. doi:10.1093/biomet/asaf008

Examples

# Longitudinal data example

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN <- make.dN(pbc2, data.id, X.grid, s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
            n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)

b = 1.7
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
              n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)

t = 2

h.hat.x.t.vec(X.grid, s.grid, n.s.grid, alpha.hat, t, b, Z.i.ds, X.pos, n)

# Time-invariant data example:
# pbc2 dataset, single event per individual version:
# 312 observations, most recent event per individual.
# Use landmarking to produce comparable curve with KM.
library(survival)
Landmark <- 3 #set the landmark to 3 years
pbc2.use<- to.id(pbc2) # keep only the most recent row per patient
pbcT1 <- pbc2.use[which(pbc2.use$year< Landmark  & pbc2.use$years> Landmark),]

timesS2 <- seq(Landmark,14,by=0.5)
b=0.9
h.hat.x<- get.h.hat.x(pbcT1, 'albumin', s.grid=s.grid,  T.name = 'years',  s.name = 'year',
                                                    delta.name = 'status2', 2, b)
s.grid.2  = seq(Landmark, 14,  length=99)
S.hat.alb<- make.S.hat(    (s.grid.2[2]-s.grid.2[1])/1.35 , h.hat.x)
kma2<- survfit(Surv(years , status2) ~ 1, data = pbcT1)

#Plot the survival functions:
plot(s.grid.2, S.hat.alb, type="l", ylim=c(0,1), xlim=c(Landmark,14), ylab="Survival probability",
            xlab="years",lwd=2, main="HQM and KM survival functions, conditional on albumin=2,
            for the time-invariant pbc dataset")
lines(kma2$time, kma2$surv, type="s",lty=2, lwd=2, col=2)
legend("bottomleft", c("HQM est.", "Kaplan-Meier"), lty=c(1,2), col=1:2, lwd=2, cex=1.7)


############# Indexed HR estimator example 1: ##############
### The example is based on indexing two covariates: albumin and bilirubin

index.use<-84
t.grid<-s.grid[1:index.use] #non need to plot the last year

b.alb = 1.5 # results from CV rule of Bagkavos et al. (2025), i.e.:
 # b.candidates = seq(1.5, 3, 0.05)
 # CV.b.alb = CV.b(pbc2, 'albumin', 'years', 'year', 'status2', size.I.j, b.candidates)
 # b.alb =  CV.b.alb[[2]][which.min(CV.b.alb[[1]])]

b.bil = 4 # results from CV rule of Bagkavos et al. (2025) :
  # b.candidates.bil = seq(3, 4, 0.05)
  # CV.b.bil = CV.b(pbc2, 'serBilir', 'years', 'year', 'status2', size.I.j, b.candidates.bil)
  # b.bil =  CV.b.alb[[2]][which.min(CV.b.alb[[1]])] - result = 4

x.alb = 1 # refers to zero mean variables
          #corresponds to approximately normal albumin level
x.bil = 1.9 # refers to zero mean variable - corresponds to
          #elevated bilirubin levelm approximately 7mg/dl

theta.alb  <- 0.0702 # results from the indexing process
theta.bil <- 0.0856 # results from the indexing process

# First, mean adjust the albumin covariate
X.alb.obs = pbc2$albumin - mean(pbc2$albumin)
X.alb.i.T.i = data.id$albumin - mean(data.id$albumin)

X.alb.grid = seq(min(X.alb.obs), max(X.alb.obs), (max(X.alb.obs)-min(X.alb.obs))/( n.X.grid-1))
X.alb.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.alb.obs, s.obs)
X.alb.pos <- findInterval(X.alb.i.s, X.alb.grid)
dN.alb <- make.dN(pbc2, data.id, X.alb.grid, s.grid, T.i, X.alb.i.T.i, delta.i)
Z.ds.alb <- make.Z.ds(pbc2, data.id, X.alb.i.s, X.alb.grid, s.grid,
            n.s.grid, n.X.grid, s.pos, X.alb.pos, T.name = 'years', n)

alpha.hat.alb<-get.alpha.hat(dN.alb, Z.ds.alb, b.alb, X.alb.grid, K=K.Epanechnikov )
Z.i.ds.alb <- make.Z.i.ds(pbc2, data.id, X.alb.i.s, X.alb.grid, s.grid, n.s.grid, n.X.grid,
                                          s.pos, X.alb.pos, T.name = 'years', n)

#calculate HQM estimator for original albumin covariate:
h.hat.alb<- h.hat.x.t.vec(X.alb.grid, s.grid, n.s.grid, alpha.hat.alb, x.alb, b.alb, Z.i.ds.alb,
                                                                  X.alb.pos, n)

y.grid.alb<-h.hat.alb[1:index.use]

# Now, mean adjust the bilirubon covariate:
X.bil.obs = pbc2$serBilir - mean(pbc2$serBilir)
X.bil.i.T.i = data.id$serBilir - mean(data.id$serBilir)

X.bil.grid = seq(min(X.bil.obs), max(X.bil.obs), (max(X.bil.obs)-min(X.bil.obs))/( n.X.grid-1))
X.bil.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.bil.obs, s.obs)
X.bil.pos <- findInterval(X.bil.i.s, X.bil.grid)
dN.bil <- make.dN(pbc2, data.id, X.bil.grid, s.grid, T.i, X.bil.i.T.i, delta.i)
Z.ds.bil <- make.Z.ds(pbc2, data.id, X.bil.i.s, X.bil.grid, s.grid,
            n.s.grid, n.X.grid, s.pos, X.bil.pos, T.name = 'years', n)

alpha.hat.bil<-get.alpha.hat(dN.bil, Z.ds.bil, b.bil, X.bil.grid, K=K.Epanechnikov )
Z.i.ds.bil <- make.Z.i.ds(pbc2, data.id, X.bil.i.s, X.bil.grid, s.grid, n.s.grid, n.X.grid,
                                        s.pos, X.bil.pos, T.name = 'years', n)

#calculate HQM estimator for original bilirubin covariate:
h.hat.bil<- h.hat.x.t.vec(X.bil.grid, s.grid, n.s.grid, alpha.hat.bil, x.bil, b.bil, Z.i.ds.bil,
                                                                    X.bil.pos, n)
y.grid.bil<-h.hat.bil[1:index.use]

# calculate the indexed variables:
theta.T.X.obs = theta.alb * X.alb.obs + 	theta.bil *X.bil.obs
theta.T.X.i.T.i = theta.alb * X.alb.i.T.i + theta.bil *X.bil.i.T.i
b = 0.4 # results from CV rule in Bagkavos et al (2025)
x = theta.alb * x.alb + theta.bil *x.bil

X.grid = seq(min(theta.T.X.obs), max(theta.T.X.obs), 
         (max(theta.T.X.obs)-min(theta.T.X.obs))/( n.X.grid-1))

theta.T.X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, theta.T.X.obs, s.obs)
X.pos <- findInterval(theta.T.X.i.s, X.grid)
dN <- make.dN(pbc2, data.id, X.grid=X.grid, s.grid=s.grid, T.i, theta.T.X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, theta.T.X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos,
                                                             X.pos, 'years', n)
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )
Z.i.ds <- make.Z.i.ds(pbc2, data.id, theta.T.X.i.s, X.grid, s.grid, n.s.grid, n.X.grid,
                                                      s.pos, X.pos,'years', n)

# Now calculate indexed HR estimator:
h.hat.x.theta<- h.hat.x.t.vec(X.grid, s.grid, n.s.grid, alpha.hat, x, b, Z.i.ds, X.pos, n)
y.grid2<-h.hat.x.theta[1:index.use]


############## Plot the results:
plot(t.grid, y.grid2, type="l", ylim=c(0,2), xlab="Years", ylab="Hazard rate", lwd=2 )
lines(t.grid, y.grid.bil , lty=2, col=4, lwd=2)
lines(t.grid, y.grid.alb, lty=2, col=2, lwd=2)

legend("topleft", lty=c(1, 2, 2), lwd=c(2,2,2),
       c("Indexed (Bilirubin and Albumin) HQM hazard rate est.",
         "HQM hazard rate est. conditional on Albumin=2mg/dl",
         "HQM hazard rate est. conditional on Bilirubin=7mg/dl" ),
       col=c(1,2,4), cex=2, bty="n")

############# Indexed HR estimator example 2: ##############
### The example is based on indexing two covariates:  bilirubin and prothrombin time

index.use<-84
t.grid<-s.grid[1:index.use]

b.pro = 3
b.bil = 4

x.pro = 1.5 # refers to zero mean variables - slightly high
x.bil = 1.9 # refers to zero mean variable - high

theta.pro  <- 0.349 #0.5 #0.149
theta.bil <- 0.0776 #0.10

############ prothrombin time #############
ind1 <- which( is.na( pbc2$prothrombin ) )
X.pro.obs =   pbc2$prothrombin
X.pro.obs[ind1]<-mean(pbc2$prothrombin, na.rm=TRUE)
X.pro.obs <- X.pro.obs - mean(X.pro.obs)

ind2 <- which( is.na( data.id$prothrombin ) )
X.pro.i.T.i =   data.id$prothrombin
X.pro.i.T.i[ind2]<-mean(data.id$prothrombin, na.rm=TRUE)
X.pro.i.T.i <- X.pro.i.T.i - mean(X.pro.i.T.i)

X.pro.grid = seq(min(X.pro.obs), max(X.pro.obs), (max(X.pro.obs)-min(X.pro.obs))/( n.X.grid-1))
X.pro.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.pro.obs, s.obs)
X.pro.pos <- findInterval(X.pro.i.s, X.pro.grid)
dN.pro <- make.dN(pbc2, data.id, X.pro.grid, s.grid, T.i, X.pro.i.T.i, delta.i)
Z.ds.pro <- make.Z.ds(pbc2, data.id, X.pro.i.s, X.pro.grid, s.grid,
                n.s.grid, n.X.grid, s.pos, X.pro.pos, T.name = 'years', n)

alpha.hat.pro<-get.alpha.hat(dN.pro, Z.ds.pro, b.pro, X.pro.grid, K=K.Epanechnikov )
Z.i.ds.pro <- make.Z.i.ds(pbc2, data.id, X.pro.i.s, X.pro.grid, s.grid, n.s.grid, n.X.grid, s.pos,
                                                          X.pro.pos, T.name = 'years', n)
h.hat.pro<- h.hat.x.t.vec(X.pro.grid, s.grid, n.s.grid, alpha.hat.pro, x.pro, b.pro, Z.i.ds.pro, 
                            X.pro.pos, n)
y.grid.pro<-h.hat.pro[1:index.use]


############ Bilirubin #############
X.bil.obs = pbc2$serBilir - mean(pbc2$serBilir)
X.bil.i.T.i = data.id$serBilir - mean(data.id$serBilir)

X.bil.grid = seq(min(X.bil.obs), max(X.bil.obs), (max(X.bil.obs)-min(X.bil.obs))/( n.X.grid-1))
X.bil.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.bil.obs, s.obs)
X.bil.pos <- findInterval(X.bil.i.s, X.bil.grid)
dN.bil <- make.dN(pbc2, data.id, X.bil.grid, s.grid, T.i, X.bil.i.T.i, delta.i)
Z.ds.bil <- make.Z.ds(pbc2, data.id, X.bil.i.s, X.bil.grid, s.grid,
            n.s.grid, n.X.grid, s.pos, X.bil.pos, T.name = 'years', n)

alpha.hat.bil<-get.alpha.hat(dN.bil, Z.ds.bil, b.bil, X.bil.grid, K=K.Epanechnikov )
Z.i.ds.bil <- make.Z.i.ds(pbc2, data.id, X.bil.i.s, X.bil.grid, s.grid, n.s.grid, n.X.grid, s.pos,
                                                          X.bil.pos, T.name = 'years', n)
h.hat.bil<- h.hat.x.t.vec(X.bil.grid, s.grid, n.s.grid, alpha.hat.bil, x.bil, b.bil, Z.i.ds.bil, 
                          X.bil.pos, n)
y.grid.bil<-h.hat.bil[1:index.use]

######### Index prothrombin time and bilirubin
theta.T.X.obs = theta.pro * X.pro.obs + 	theta.bil *X.bil.obs
theta.T.X.i.T.i = theta.pro * X.pro.i.T.i + theta.bil *X.bil.i.T.i

b = 1.5
x = theta.pro * x.pro + theta.bil *x.bil

X.grid = seq(min(theta.T.X.obs), max(theta.T.X.obs), 
          (max(theta.T.X.obs)-min(theta.T.X.obs))/( n.X.grid-1))

theta.T.X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, theta.T.X.obs, s.obs)
X.pos <- findInterval(theta.T.X.i.s, X.grid)
dN <- make.dN(pbc2, data.id, X.grid=X.grid, s.grid=s.grid, T.i, theta.T.X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, theta.T.X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos, X.pos,
                                                                          'years', n)
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )
Z.i.ds <- make.Z.i.ds(pbc2, data.id, theta.T.X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos, 
                        X.pos,'years', n)
h.hat.x.theta<- h.hat.x.t.vec(X.grid, s.grid, n.s.grid, alpha.hat, x, b, Z.i.ds, X.pos, n)
y.grid2<-h.hat.x.theta[1:index.use]



##############Plot the results
plot(t.grid, y.grid2, type="l", ylim=c(0,2), xlab="Years", ylab="Hazard rate", lwd=2 )
lines(t.grid, y.grid.bil , lty=2, col=4, lwd=2)

lines(t.grid, y.grid.pro, lty=2, col=2, lwd=2)
legend("topleft", lty=c(1, 2, 2), lwd=c(2,2,2),
       c("Indexed (Bilirubin and PT) HQM hazard rate est.",
        "HQM hazard rate est. conditional on PT=15.5s",
        "HQM hazard rate est. conditional on Bilirubin=7mg/dl" ),
        col=c(1,2,4), cex=2, bty="n")

Local linear future conditional hazard rate estimation at a single time point

Description

Calculates the (indexed) local linear future conditional hazard rate for a marker value x and a time value t.

Usage

h.tilde.x.t(X.grid, s.grid, X.pos, n.s.grid, alpha.hat, x,t, b, Z.i.ds,n, Z.ds)

Arguments

Details

Function h.tilde.x.t implements the future local linear conditional hazard estimator

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i(\hat \theta^T X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x-\hat \theta^T X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- \hat \theta^T X_i(s))\mathrm {d}s},

for every x on the marker grid, where, for a positive integer p, \hat \theta^T = (\hat \theta_1, \dots, \hat \theta_p ) is the vector of the estimated indexing parameters, X_i = (X_{1, i}, \dots, X_{i,p}) is a vector of markers for indexing, Z_i is the exposure and \alpha(z) is the marker-only hazard, see get.alpha.hat for more details. For p=1 and \hat \theta = 1 the above estimator becomes the HQM hazard rate estimator conditional on a single covariate,

\hat{h}_x(t) = \frac{\sum_{i=1}^n \int_0^T\hat{\alpha}_i( X_i(t+s))Z_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s}{\sum_{i=1}^n\int_0^TZ_i(t+s)Z_i(s)K_{b}(x- X_i(s))\mathrm {d}s},

defined in equation (2) of Bagkavos et al. (2025) doi:10.1093/biomet/asaf008. In the place of K_b(), h.tilde.x.t uses the kernel

K_{x,b}(u)= \frac{K_b(u)-K_b(u)u^T D^{-1}c_1}{c_0 - c_1^T D^{-1} c_1},

where K_b() = b^{-1}K(./b) with K being an ordinary kernel, e.g. the Epanechnikov kernel, c_1 = (c_{11}, \dots, c_{1d})^T, D = (d_{ij})_{(d+1) \times (d+1)} with

c_0 = \sum_{i=1}^n \int_0^T K_b(x-\hat \theta^T X_i(s)) Z_i(s)ds,

c_{ij} = \sum_{i=1}^n \int_0^T K_b(x-\hat \theta^T X_i(s))\{x-\hat \theta^T X_{ij}(s)\} Z_i(s)ds,

d_{jk} = \sum_{i=1}^n \int_0^T K_b(x-\hat \theta^T X_i(s))\{x-\hat \theta^T X_{ij}(s)\}\{x-\hat \theta^T X_{ik}(s)\} Z_i(s)ds,

see also Nielsen (1998) doi:10.1080/03461238.1998.10413997.

The future conditional hazard is defined as

h_{x,T}(t) = P\left(T_i\in (t+T, t+T+dt)| \theta_0^T X_i(T)=x, T_i > t+T\right),

where T_i is the survival time and X_i the marker of individual i observed in the time frame [0,T].

Value

A single numeric value of \hat h_x(t).

References

Bagkavos, I., Isakson, R., Mammen, E., Nielsen, J., and Proust–Lima, C. (2025). Biometrika, 112(2), asaf008. doi:10.1093/biomet/asaf008

Nielsen (1998), Marker dependent kernel hazard estimation from local linear estimation, Scandinavian Actuarial Journal, pp. 113-124. doi:10.1080/03461238.1998.10413997

See Also

get.alpha.hat, c2.u.integrand

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN <- make.dN(pbc2, data.id, X.grid=X.grid, s.grid=s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos, X.pos,
            'years', n)

b = 1.7
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid, n.s.grid, n.X.grid, s.pos, X.pos,
              'years', n)

x = 2
t = 2
h.tilde = h.tilde.x.t(X.grid, s.grid, X.pos, n.s.grid, alpha.hat, x, t, b, Z.i.ds, n, Z.ds)

Simulated true hazard (bootstrap average)

Description

Compute the Monte Carlo / bootstrap averaged alpha (true hazard) across a list of datasets.

Usage

h.true.simulated(bootstrap.samples, id.name, n.s.grid, X.1.name = X.1.name, 
       X.2.name = X.2.name, T.name = T.name, 
       s.name = s.name, delta.name = delta.name, theta, b)

Arguments

Value

Numeric vector (row means of alpha estimates).

Examples

 
X.1.name <- 'albumin'
X.2.name <-  'serBilir'
T.name <- 'years' 
s.name <- 'year' 
delta.name <- 'status2'
id.name<-'id'

n.s.grid<-50 #user supplied grid length on the s.name argument

theta.1  <- 0.0702  #indexing parameter for marker 1
theta.2 <- 0.0856  #indexing parameter for marker 2

x.1 = 0 # conditioning level for marker 1, assumes zero mean variable
x.2 = 1.9 # conditioning level for marker 2, assumes zero mean variable
b = 0.42 # The result, on the indexed marker 'indmar' of 
         #\code{CV.b(pbc2,'indmar','years','year','status2',size.I.j=26,seq(0.2,0.4,by=0.01))}
x = theta.1 * x.1 + theta.2 *x.2

data <- pbc2[,c(id.name, X.1.name, X.2.name, T.name, s.name, delta.name)]
n.obs <- length(data$id)
n<-max(  as.double(data[,'id']) )

#  Create bootstrap samples by group 
set.seed(1)  
B<-100 #set bootstrap iterations here
bootstrap.samples<-list()
for(j in 1:B)
{
  rows.star<-c()
  ids.star<-c()
  ids.drawn <- sample (n, replace = TRUE)  
  for(l in 1:n)
  {
    rows.l<-which(data[,id.name]==ids.drawn[l])
    rows.star<-c(rows.star, rows.l)
    ids.l<-rep(ids.drawn[l], times=length(rows.l))
    ids.star<-c(ids.star, ids.l)
  }
  data.star.unsorted<-data[rows.star,]
  data.star.unsorted<-data[rows.star,]
  bootstrap.samples[[j]]<- data.star.unsorted[order(data.star.unsorted$id),]  
}

#compute the bootstrap simulated true HR function:
h.true<- h.true.simulated(bootstrap.samples, id.name='id', n.s.grid=n.s.grid, 
              X.1.name=X.1.name, X.2.name= X.2.name, 
              T.name = T.name, s.name = s.name,  
              delta.name = delta.name, theta = c(theta.1,  theta.2), b)

 

IID decomposition for time-dependent AUC with IPCW

Description

Internal function implementing the influence function (iid) decomposition of the time-dependent Area Under the Curve (AUC) in the presence of right-censoring, using inverse probability of censoring weighting (IPCW).

This function reproduces the internal behavior of the corresponding routines from the timeROC package. It supports both standard survival settings and competing risks settings through automatic dispatch.

Usage

  iid.AUC.hat(t, n, cause, F.1.t, S.t,
                            G.hat, T.i, delta.i, pi.i, int.dM.C.over.S)

Arguments

Details

The function computes the iid (influence function) representation of the time-dependent AUC estimator based on IPCW methodology.

Two cases are handled internally:

• Survival setting: only one event type is present (no competing risks).

• Competing risks setting: multiple event types are present, and the AUC is defined using different control definitions.

The implementation relies on pairwise comparisons between cases and controls, weighted by inverse probability of censoring weights. The iid representation is derived from a U-statistic formulation and incorporates the influence of the censoring distribution via martingale-based corrections.

The computational complexity is quadratic in the sample size.

Value

A list with the following components:

Note

This is an internal utility function and is not intended for direct use by end users. It is provided for reproducibility and to avoid dependency on archived packages.

Author(s)

Adapted from the timeROC package.

References

Blanche, P., Dartigues, J. F., and Jacqmin-Gadda, H. (2013). Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Statistics in Medicine, 32(30), 5381–5397.

Uno, H., Cai, T., Tian, L., and Wei, L. J. (2007). Evaluating prediction rules for t-year survivors with censored regression models. Journal of the American Statistical Association, 102(478), 527–537.

See Also

timeROC (in the original timeROC package), pec::ipcw, iid.G.hat

iid.G.hat function implementation

Description

Calculates the AUC for the HQM estimator.

Usage

iid.G.hat(T.i,delta.i) 

Arguments

Details

test

Value

A vector of two values: the landmark time of the calculation and the AUC value.

See Also

BS.hat.HQM

Survival function from a hazard

Description

Creates a survival function from a hazard rate which was calculated on a grid.

Usage

make.S.hat(du, h.hat.x)

Arguments

Details

The function make.S.hat calculates the survival function

S(t) = \exp (-\int_0^t h(t) dt),

where h is the hazard rate. Here, a discritisation via an equidistant grid \{ t_i \} on [0,t] is used to calculate the integral and it is assumed that h has been calculated for exactly these time points t_i .

Value

A vector of values S(t_i).

Examples

make.S.hat(0.1, rep(0.1,10))

Occurance and Exposure on grids

Description

Auxiliary functions that help automate the process of calculating integrals with occurances or exposure processes.

Usage

make.dN(data, data.id, X.grid, s.grid, T.i, X.i.T.i, delta.i)
make.dN.i(s.grid, n.s.grid, T.i, delta.i, n)
make.Z.ds(data, data.id, X.i.s, X.grid, s.grid,
        n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)
make.Z.i.ds(data, data.id, X.i.s, X.grid, s.grid,
        n.s.grid, n.X.grid, s.pos,X.pos, T.name = 'years', n)

Arguments

Details

Implements matrices for the computation of integrals with occurences and exposures of the form

\int f(s)Z(s)Z(s+t)ds, \int f(s) Z(s)ds, \int f(s) dN(s).

where N is a 0-1 counting process, Z the exposure and f an arbitrary function.

Value

The functions make.dN and make.Z.ds return a matrix on the time grid and marker grid for occurence and exposure, respectively, while make.dN.i and make.Z.i.ds return a matrix on the time grid for evey individual again for occurence and exposure, respectively.

See Also

h.hat.x.t, g.hat.t.x, get.alpha.hat

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN <- make.dN(pbc2, data.id, X.grid, s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
            n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)
Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
              n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)
dN.i  <- make.dN.i(s.grid, n.s.grid, T.i, delta.i, n)

Mayo Clinic Primary Biliary Cirrhosis Data

Description

Followup of 312 randomised patients with primary biliary cirrhosis, a rare autoimmune liver disease, at Mayo Clinic.

Usage

  pbc2 

Format

A data frame with 1945 observations on the following 20 variables.

id

patients identifier; in total there are 312 patients.

years

number of years between registration and the earlier of death, transplantion, or study analysis time.

status

a factor with levels alive, transplanted and dead.

drug

a factor with levels placebo and D-penicil.

age

at registration in years.

sex

a factor with levels male and female.

year

number of years between enrollment and this visit date, remaining values on the line of data refer to this visit.

ascites

a factor with levels No and Yes.

hepatomegaly

a factor with levels No and Yes.

spiders

a factor with levels No and Yes.

edema

a factor with levels No edema (i.e., no edema and no diuretic therapy for edema), edema no diuretics (i.e., edema present without diuretics, or edema resolved by diuretics), and edema despite diuretics (i.e., edema despite diuretic therapy).

serBilir

serum bilirubin in mg/dl.

serChol

serum cholesterol in mg/dl.

albumin

albumin in gm/dl.

alkaline

alkaline phosphatase in U/liter.

SGOT

SGOT in U/ml.

platelets

platelets per cubic ml / 1000.

prothrombin

prothrombin time in seconds.

histologic

histologic stage of disease.

status2

a numeric vector with the value 1 denoting if the patient was dead, and 0 if the patient was alive or transplanted.

References

Fleming, T. and Harrington, D. (1991) Counting Processes and Survival Analysis. Wiley, New York.

Therneau, T. and Grambsch, P. (2000) Modeling Survival Data: Extending the Cox Model. Springer-Verlag, New York.

Examples

  summary(pbc2)
  

Event data frame

Description

Creates a data frame with only one entry per individual from a data frame with time dependent data. The resulting data frame focusses on the event time and the last observed marker value.

Usage

to.id(data)

Arguments

Details

The function to.id uses a data frame of time dependent marker data to create a smaller data frame with only one entry per individual, the last observed marker value and the event time. Note that the column indicating the individuals must have the name id. Note also that this data frame is similar to pbc2.id from the JM package with the difference that the last observed marker value instead of the first one is captured.

Value

A data frame with only one entry per individual.

Examples

data.id = to.id(pbc2)

Precomputation for wild bootstrap

Description

Implements key components for the wild bootstrap of the hqm estimator in preparation for obtaining confidence bands.

Usage

wild.bootstrap.summands(g.hat, alpha.hat, dN.i, Z.i.ds, n.s.grid, X.grid, s.grid, 
                        t, b, X.pos, x, n)

Arguments

Details

The function implements

A_B(t) = \frac{1}{\sqrt{n}} \sum_{i=1}^n \int^{T}_0 \hat{g}_{i,t,x_*}(X_i(s)) V_i\{dN_i(s) - \hat{\alpha}_i(X_i(s))Z_i(s)ds\},

and

B_B(t) = \frac{1}{\sqrt{n}}\sum_{i = 1}^n V_i\{\hat{\Gamma}(t,x_*)^{-1}W_i(t,x_*) - \hat{h}_{x_*}(t)\},

where V \sim N(0,1),

W_i(t) =\int_0^T\hat{\alpha}_i(X_i(t+s))Z_i(t+s)Z_i(s)K_b(x_*,X_i(s))\mathrm {d}s,

and

\hat{\Gamma}(t,x) = \frac{1}{n} \sum_{i = 1}^n \int_{0}^{T-t} Z_i(t+s)Z_i(s) K_b(x,X_i(s))ds,

with Z being the exposure and X the marker.

Value

A list of 5 items. The first two are vectors for calculating A_B and the third one a vector for B_B. The 4th one is the value of the hqm estimator that can also be obtained by h.hat.x.t and the last one is the value of \Gamma.

See Also

confidence.bands

Examples

data.id = to.id(pbc2)
n.s.grid <- n.X.grid <- 100
n = max(as.numeric(pbc2$id))
s.obs = pbc2$year
X.obs = pbc2$serBilir
X.i.T.i = data.id$serBilir
T.i <- data.id$years
delta.i <- data.id$status2
s.grid = seq(0, max(s.obs), max(s.obs)/( n.s.grid-1))
X.grid = seq(min(X.obs), max(X.obs), (max(X.obs)-min(X.obs))/( n.X.grid-1))

X.i.s = X.i.s.interpolate(s.grid, data.id$id, pbc2$id, X.obs, s.obs)

X.pos <- findInterval(X.i.s, X.grid)
s.pos = rep(1:length(s.grid), n)

dN <- make.dN(pbc2, data.id, X.grid, s.grid, T.i, X.i.T.i, delta.i)
Z.ds <- make.Z.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
            n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)

b = 1.7
alpha.hat<-get.alpha.hat(dN, Z.ds, b, X.grid, K=K.Epanechnikov )

Z.i.ds <- make.Z.i.ds(pbc2, data.id, X.i.s, X.grid, s.grid,
              n.s.grid, n.X.grid, s.pos, X.pos, T.name = 'years', n)
dN.i  <- make.dN.i(s.grid, n.s.grid, T.i, delta.i, n)

t = 2
x = 2

g.hat = g.hat.t.x(X.grid, s.grid, n.s.grid, X.pos, x, t, b, Z.i.ds, Z.ds, n)

summands = wild.bootstrap.summands(g.hat, alpha.hat, dN.i, Z.i.ds, n.s.grid, X.grid, 
                  s.grid, t, b, X.pos, x, n)
summands