Package 'HQM'

AUC for the High Quality Marker estimator

Description

Calculates the AUC for the HQM estimator.

Usage

auc.hqm(xin, est, landm, th, event_time_name, status_name)

Arguments

xin	A data frame containing event times and the patient status.
est	The HQM estimator values, typically the output of get_h_x.
landm	Landmark time.
th	Time horizon.
event_time_name	The column name of the event times in the xin data frame.
status_name	The column name of the status variable in the xin frame.

Details

The function auc.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.hqm

Examples

library(survival)
Landmark <- 2
pbcT1 <- pbc2[which(pbc2$year< Landmark  & pbc2$years> Landmark),]
ls<-50 # 50 grid points to evaluate the estimates
s.out<- pbcT1[, 'year']
s.out.use <-  seq(0, max(s.out), max(s.out)/( ls-1))
 
timesS2 <- seq(Landmark,14,by=0.5)
b=0.9
arg1<- get_h_x(pbcT1, 'albumin', br_s = s.out.use,  event_time_name = 'years',  
                time_name = 'year', event_name = 'status2', 2, 0.9) 
br_s2  = seq(Landmark, 14,  length=ls-1)
sfalb2<- make_sf(    (br_s2[2]-br_s2[1])/4 , arg1)
tHor <- 1.5
auc.hq.use<-auc.hqm(pbcT1, sfalb2, Landmark,tHor,  
              event_time_name = 'years', status_name = 'status2')
auc.hq.use

---

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

b_selection(data,   marker_name, event_time_name = 'years',
            time_name = 'year', event_name = 'status2', I, b_list)

Arguments

data	A data frame of time dependent data points. Missing values are allowed.
marker_name	The column name of the marker values in the data frame data.
event_time_name	The column name of the event times in the data frame data.
time_name	The column name of the times the marker values were observed in the data frame data.
event_name	The column name of the events in the data frame data.
I	Number of observations leave out for a K cross validation.
b_list	Vector of bandwidths that need to be tested.

Details

The function b_selection 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

b_selection_prep_g, Q1, R_K, prep_cv, dataset_split

Examples

I = 26
# For Albumin marker:
b_list = seq(0.9, 1.7, 0.1)
 

b_scores_alb = b_selection(pbc2,  'albumin', 'years', 'year', 'status2', I, b_list)
b_scores_alb[[2]][which.min(b_scores_alb[[1]])]

# For Bilirubin marker:
b_list = seq(3, 4, 0.1)
b_scores_bil = b_selection(pbc2, 'serBilir', 'years', 'year', 'status2', I, b_list)
b_scores_bil[[2]][which.min(b_scores_bil[[1]])]
b_scores_bil

---

Cross validation index parameter selection

Description

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

Usage

b_selection_index_optim(in.par, data, marker_name1, marker_name2,  
          event_time_name = 'years', time_name = 'year', event_name = 'status2', I, b)

Arguments

in.par	Vector of candidate values for the index parameters.
data	A data frame of time dependent data points. Missing values are allowed.
marker_name1	The column name of the first marker values in the data frame data.
marker_name2	The column name of the second marker values in the data frame data.
event_time_name	The column name of the event times in the data frame data.
time_name	The column name of the times the marker values were observed in the data frame data.
event_name	The column name of the events in the data frame data.
I	Number of observations leave out for a K cross validation.
b	scalar bandwidth for the HQM estimator.

Details

The function b_selection_index_optim 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

b_selection_prep_g, Q1, R_K, prep_cv, dataset_split

Examples

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

I = 26
b_list = seq(1.1, 1.2, by=0.1)

for(i in 1:length(b_list))
{
  res<- optim(c(0.5, 0.5), b_selection_index_optim,  data=pbc2, marker_name1='albumin', 
              marker_name2= 'serBilir', event_time_name = 'years', time_name = 'year', 
              event_name = 'status2', I=26, b=b_list[i], method="Nelder-Mead")
  cat("i= ", i, " ", res$par, " ", res$value, "count calls to fn = ", res$counts, " converge? ", 
        res$convergence,   "\n")
  res
}

---

Preparations for bandwidth selection

Description

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

Usage

b_selection_prep_g(h_mat, int_X, size_X_grid, n, Yi)

Arguments

h_mat	A matrix of the estimator for the future conditional hazard rate for all values x and t.
int_X	Vector of the position of the observed marker values in the grid for marker values.
size_X_grid	Numeric value indicating the number of grid points for marker values.
n	Number of individuals.
Yi	A matrix made by make_Yi indicating the exposure.

Details

The function b_selection_prep_g 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

b_selection

Examples

pbc2_id = to_id(pbc2)
size_s_grid <- size_X_grid <- 100
n = max(as.numeric(pbc2$id))
s = pbc2$year
X = pbc2$serBilir
XX = pbc2_id$serBilir
ss <- pbc2_id$years
delta <- pbc2_id$status2
br_s = seq(0, max(s), max(s)/( size_s_grid-1))
br_X = seq(min(X), max(X), (max(X)-min(X))/( size_X_grid-1))

X_lin = lin_interpolate(br_s, pbc2_id$id, pbc2$id, X, s)

int_X <- findInterval(X_lin, br_X)
int_s = rep(1:length(br_s), n)

N <- make_N(pbc2, pbc2_id, breaks_X=br_X, breaks_s=br_s, ss, XX, delta)
Y <- make_Y(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, 
            int_s, int_X, 'years', n)

b = 1.7
alpha<-get_alpha(N, Y, b, br_X, K=Epan )

Yi <- make_Yi(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, int_s, 
              int_X,'years', n)
Ni  <- make_Ni(breaks_s=br_s, size_s_grid, ss, delta, n)

t = 2

h_xt_mat = t(sapply(br_s[1:99], function(si){
          h_xt_vec(br_X, br_s, size_s_grid, alpha, t, b, Yi, int_X, n)}))

b_selection_prep_g(h_xt_mat, int_X, size_X_grid, n, Yi)

---

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

Boot.hqm(in.par, data, data.id, ls, X1, XX1,  event_time_name = 'years', 
                 time_name = 'year', event_name = 'status2', b, t)

Arguments

in.par	Numeric vector, the values of the indexing parameters.
data	Bootstrap data.frame.
data.id	Id-level data.frame (result of to_id).
ls	user supplied grid length on the time_name argument.
X1	List of vectors for indexing: each vector corresponds to a biomarker and contains one summary measurement per individual.
XX1	List of vectors for indexing: each vector corresponds to a single biomarker and contains its longitudinal measurements across all individuals/time points.
event_time_name	Name of event time column.
time_name	Name of observation time column.
event_name	Name of event indicator column.
b	Bandwidth parameter.
t	Conditioning level.

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

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

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

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

marker_name1 <- 'albumin'
marker_name2 <-  'serBilir'
event_time_name <- 'years' 
time_name <- 'year' 
event_name <- 'status2'
id<-'id'
ls<-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:
X1t=data.use[,marker_name1] -mean(data.use[, marker_name1])
XX1t=data.use.id[,marker_name1] -mean(data.use.id[, marker_name1])
X2t=data.use[,marker_name2]  -mean(data.use[, marker_name2])
XX2t=data.use.id[,marker_name2] -mean(data.use.id[, marker_name2])

X1=list(X1t, X2t)
XX1=list(XX1t, XX2t)
boot.haz<-Boot.hqm (c(par.alb,par.bil), data.use, data.use.id, ls=ls, X1, XX1,  
                    event_time_name = 'years', time_name = 'year', event_name = 'status2',   b, t)

---

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 Boot.hqm. The output is used as input in functions Quantile.Index.CIs and Pivot.Index.CIs

Usage

Boot.hrandindex.param(B, Boot.samples, marker_name1, marker_name2,event_time_name, 
                             time_name, event_name, b, t, true.haz, v.param, n.est.points)

Arguments

B	Integer. Number of bootstrap iterations.
Boot.samples	A list of bootstrap datasets. Each element corresponds to one replicate.
marker_name1	Character string. Name of the first longitudinal marker.
marker_name2	Character string. Name of the second longitudinal marker.
event_time_name	Name of the event time variable in the data.
time_name	Name of the time variable for the longitudinal marker measurements.
event_name	Name of the event indicator variable.
b	Numeric. Bandwidth parameter used in hazard estimation.
t	Numeric. Evaluation point for the conditional hazard.
true.haz	Numeric vector. The true or reference hazard used in the optimisation criterion.
v.param	Numeric vector. Starting values of the indexing parameters for the optimisation of the index coefficients.
n.est.points	Integer. Number of estimation grid points for the hazard curve.

Details

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

1. Extracts the bootstrap sample data.use.

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

3. Estimates index parameters by minimising index_optim using optim.

4. Computes the bootstrap hazard estimate via Boot.hqm.

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

Value

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

See Also

Boot.hqm, index_optim, to_id

Examples

marker_name1 <- 'albumin'
marker_name2 <-  'serBilir'
event_time_name <- 'years' 
time_name <- 'year' 
event_name <- 'status2'
id<-'id'

par.x1  <- 0.0702  
par.x2 <- 0.0856  
t.x1 = 0 # refers to zero mean variables - slightly high
t.x2 = 1.9 # refers to zero mean variable - high
b = 0.42 
t = par.x1 * t.x1 + par.x2 *t.x2

# first simulate true HR function:
xin <- pbc2[,c(id, marker_name1, marker_name2, event_time_name, time_name, event_name)]
n <- length(xin$id)
nn<-max(  as.double(xin[,'id']) )
xin.id <- to_id(xin)


#  Create bootstrap samples by group:
set.seed(1)  
B<- 10 # 400 #50
Boot.samples<-list()
for(j in 1:B)
{
  i.use<-c()
  id.use<-c()
  index.nn <- sample (nn, replace = TRUE)  
  for(l in 1:nn)
  {
    i.use2<-which(xin[,id]==index.nn[l])
    i.use<-c(i.use, i.use2)
    id.use2<-rep(index.nn[l], times=length(i.use2))
    id.use<-c(id.use, id.use2)
  }
  xin.i<-xin[i.use,]
  xin.i<-xin[i.use,]
  Boot.samples[[j]]<- xin.i[order(xin.i$id),] #xin[i.use,]
}
true.hazard<- Sim.True.Hazard(Boot.samples, id='id', 100,  marker_name1=marker_name1,  
                    marker_name2= marker_name2, event_time_name = event_time_name, 
                    time_name = time_name,  event_name = event_name, 
                    in.par = c(par.x1,  par.x2), b)
                              
res <- Boot.hrandindex.param( B,  Boot.samples, marker_name1, marker_name2, 
          event_time_name,  time_name,  event_name ,  b = 0.4, t = 1.0,  
          true.haz = true.hazard, v.param = c(0.07, 0.08), n.est.points = 100)
          
#return bootstrap hazard rate estimators in marix format:
res

---

Brier score for the High Quality Marker estimator

Description

Calculates the Brier score for the HQM estimator.

Usage

bs.hqm(xin, est, landm, th, event_time_name, status_name)

Arguments

xin	A data frame containing event times and the patient status.
est	The HQM estimator values, typically the output of get_h_x.
landm	Landmark time.
th	Time horizon.
event_time_name	The column name of the event times in the data frame xin.
status_name	The column name of the status variable in the data frame xin.

Details

The function bs.hqm implements the Brier score calculation for the HQM estimator estimator.

Value

Scalar: the Brier score of the HQM estimator.

See Also

auc.hqm

Examples

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

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


b=0.9
arg1<- get_h_x(pbcT1, 'albumin', br_s = s.out.use, event_time_name = 'years',  
               time_name = 'year', event_name = 'status2', 2, 0.9) 
br_s2  = seq(Landmark, 14,  length=ls-1)
sfalb2<- make_sf(    (br_s2[2]-br_s2[1])/4 , arg1)


tHor <- 1.5
bs.use<-bs.hqm(pbcT1, sfalb2, Landmark,tHor, 
                event_time_name = 'years', status_name = 'status2')
bs.use

---

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 StudentizedBwB.Index.CIs.

Usage

BwB.HRandIndex.param(B, B1, Boot.samples, marker_name1, marker_name2, 
                            event_time_name, time_name, event_name, b, t, true.haz, 
                            v.param,  hqm.est, id, xin)

Arguments

B	Integer. Number of bootstrap samples.
B1	Integer. Number of bootstrap re-samples.
Boot.samples	A list of bootstrap datasets. Each element corresponds to one replicate.
marker_name1	Character string. Name of the first longitudinal marker.
marker_name2	Character string. Name of the second longitudinal marker.
event_time_name	Name of the event time variable in the data.
time_name	Name of the time variable for the longitudinal marker measurements.
event_name	Name of the event indicator variable.
b	Numeric. Bandwidth parameter used in hazard estimation.
t	Numeric. Evaluation point for the conditional hazard.
true.haz	Numeric vector. The true or reference hazard used in the optimisation criterion.
v.param	Numeric vector. Starting values of the indexing parameters for the optimisation of the index coefficients.
hqm.est	HQM estimator on the original sample.
id	label of id variable of dataset.
xin	original sample.

Details

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

1. Extracts the bootstrap sample data.use.

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

3. Estimates index parameters by minimising index_optim using optim.

4. Computes the bootstrap hazard estimate via Boot.hqm.

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

Value

A list of matrices of dimension n.est.points × 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

Boot.hqm, index_optim, to_id

Examples

# See the example for  function: StudentizedBwB.Index

---

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

compute_iid_decomposition(t, n, cause, F01t, St,
                            weights, T, delta, marker, MatInt0TcidhatMCksurEff)

Arguments

t	Numeric scalar. Time point at which the AUC is evaluated.
n	Integer. Sample size.
cause	Integer or numeric. Event type of interest.
F01t	Numeric. Estimated cumulative incidence function at time t for the event of interest.
St	Numeric. Estimated survival probability at time t.
weights	List. IPCW weights object returned by pec::ipcw. Must contain components IPCW.subjectTimes, IPCW.times, and times.
T	Numeric vector. Observed event or censoring times.
delta	Integer or numeric vector. Event indicator: 0 for censoring, cause for the event of interest, and other values for competing events.
marker	Numeric vector. Prognostic marker or risk score.
MatInt0TcidhatMCksurEff	Numeric matrix. Influence function representation of the censoring distribution, typically obtained from Compute.iid.KM().

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:

iid_representation_AUC	Numeric vector of length n. Influence function values for the AUC estimator.
iid_representation_AUCstar	Numeric vector of length n. Alternative influence function corresponding to an alternative control definition (mainly relevant in competing risks).
seAUC	Numeric scalar. Estimated standard error of the AUC.
seAUCstar	Numeric scalar. Estimated standard error corresponding to AUCstar.
computation_times	Placeholder (numeric or NA). Included for compatibility with timeROC.

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, Compute.iid.KM

---

Compute.iid.KM function implementation

Description

Calculates the AUC for the HQM estimator.

Usage

Compute.iid.KM(times,status)

Arguments

times	The vector of (censored) event-times.
status	The vector of event indicators at the corresponding value of the vector T. Censored observations must be denoted by the value 0.

Details

test

Value

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

See Also

bs.hqm

---

Confidence bands

Description

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

Usage

Conf_bands(data,   marker_name, event_time_name = 'years',
            time_name = 'year', event_name = 'status2', x, b)

Arguments

data	A data frame of time dependent data points. Missing values are allowed.
marker_name	The column name of the marker values in the data frame data.
event_time_name	The column name of the event times in the data frame data.
time_name	The column name of the times the marker values were observed in the data frame data.
event_name	The column name of the events in the data frame data.
x	Numeric value of the last observed marker value.
b	Bandwidth.

Details

The function Conf_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 prep_boot and g_xt.

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_xt, prep_boot

Examples

b = 10
x = 3
size_s_grid<-100
s = pbc2$year
br_s = seq(0, max(s), max(s)/( size_s_grid-1))


c_bands = Conf_bands(pbc2,   'serBilir', event_time_name = 'years',
                    time_name = 'year', event_name = 'status2', x, b)

J = 80
plot(br_s[1:J], c_bands$h_hat[1:J], type = "l", ylim = c(0,1), ylab = 'Hazard', xlab = 'Years')

lines(br_s[1:J], c_bands$I_p_up[1:J], col = "red")
lines(br_s[1:J], c_bands$I_p_do[1:J], col = "red")
lines(br_s[1:J], c_bands$I_nu[1:J], col = "blue")
lines(br_s[1:J], c_bands$I_nd[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

dataset_split(I, data)

Arguments

data	A data frame of time dependent data points. Missing values are allowed.
I	The number of individuals that should be left out. Optimally, $K = n/I$ should be an integer, where $n$ is the number of individuals.

Details

The function dataset_split 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 I individuals missing in the above way.

See Also

b_selection

Examples

splitted_dataset = dataset_split(26, pbc2)

---

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

dij(b,x,y, K)

Arguments

x	A vector of design points where the kernel will be evaluated.
y	A vector of sample data points.
b	The bandwidth to use (a scalar).
K	The kernel function to use.

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}$.

---

Epanechnikov kernel

Description

Implements the Epanechnikov kernel function

Usage

Epan(x)

Arguments

x	A vector of design points where the kernel will be evaluated.

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 $x$.

---

Computation of a key component for wild bootstrap

Description

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

Usage

g_xt(br_X, br_s, size_s_grid, int_X, x, t, b, Yi, Y, n)

Arguments

br_X	Marker value grid points that will be used in the evaluatiuon.
br_s	Time value grid points that will be used in the evaluatiuon.
size_s_grid	Size of the time grid.
int_X	Position of the linear interpolated marker values on the marker grid.
x	Numeric value of the last observed marker value.
t	Numeric value of the time the function should be evaluated.
b	Bandwidth.
Yi	A matrix made by make_Yi indicating the exposure.
Y	A matrix made by make_Y indicating the exposure.
n	Number of individuals.

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

pbc2_id = to_id(pbc2)
size_s_grid <- size_X_grid <- 100
n = max(as.numeric(pbc2$id))
X = pbc2$serBilir
s = pbc2$year
br_s = seq(0, max(s), max(s)/( size_s_grid-1))
br_X = seq(min(X), max(X), (max(X)-min(X))/( size_X_grid-1))

X_lin = lin_interpolate(br_s, pbc2_id$id, pbc2$id, X, s)

int_X <- findInterval(X_lin, br_X)
int_s = rep(1:length(br_s), n)

Yi<-make_Yi(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, int_s, int_X, 
            'years', n)
Y<-make_Y(pbc2, pbc2_id, X_lin, br_X, br_s,size_s_grid,size_X_grid, int_s,int_X, 
            'years', n)

t = 2
x = 2
b = 10

g_xt(br_X, br_s, size_s_grid, int_X, x, t, b, Yi, Y, n)

---

Marker-only hazard rate

Description

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

Usage

get_alpha(N, Y, b, br_X, K=Epan )

Arguments

N	A matrix made by make_N indicating the occurences of events.
Y	A matrix made by make_Y indicating the exposure.
b	Bandwidth.
br_X	Vector of grid points for the marker values $X$.
K	Used kernel function.

Details

The function get_alpha 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 br_X.

See Also

h_xt

Examples

pbc2_id = to_id(pbc2)
size_s_grid <- size_X_grid <- 100
n = max(as.numeric(pbc2$id))
s = pbc2$year
X = pbc2$serBilir
XX = pbc2_id$serBilir
ss <- pbc2_id$years
delta <- pbc2_id$status2
br_s = seq(0, max(s), max(s)/( size_s_grid-1))
br_X = seq(min(X), max(X), (max(X)-min(X))/( size_X_grid-1))

X_lin = lin_interpolate(br_s, pbc2_id$id, pbc2$id, X, s)

int_X <- findInterval(X_lin, br_X)
int_s = rep(1:length(br_s), n)

N  <- make_N(pbc2, pbc2_id, breaks_X=br_X, breaks_s=br_s, ss, XX, delta)
Y <- make_Y(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid,
            size_X_grid, int_s, int_X, event_time = 'years', n)

b = 1.7
alpha<-get_alpha(N, Y, b, br_X, K=Epan )

---

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_x(data, marker_name, br_s, event_time_name, time_name, event_name, x, b)

Arguments

data	A data frame of time dependent data points. Missing values are allowed.
marker_name	The column name of the marker values in the data frame data.
br_s	User defined grid mesh on time_name variable
event_time_name	The column name of the event times in the data frame data.
time_name	The column name of the times the marker values were observed in the data frame data.
event_name	The column name of the events in the data frame data.
x	Numeric value of the last observed marker value.
b	Bandwidth parameter.

Details

The function get_h_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 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, h_xt, get_h_xll

Examples

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

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

b=0.9

arg1ll<-get_h_xll(pbcT1,'albumin', br_s = s.out.use, event_time_name='years',
                  time_name='year',event_name='status2',2,0.9) 
arg1lc<-get_h_x(pbcT1,'albumin', br_s = s.out.use, event_time_name='years',  
                  time_name='year',event_name='status2',2,0.9) 

#Caclulate the local contant and local linear survival functions
br_s  = seq(Landmark, 14,  length= ls-1 )
sfalb2ll<- make_sf(    (br_s[2]-br_s[1])/4 , arg1ll)
sfalb2lc<- make_sf(    (br_s[2]-br_s[1])/4 , arg1lc)

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

#Plot the survival functions:
plot(br_s, sfalb2ll,  type="l", col=1, lwd=2, ylab="Survival probability", 
     xlab="Marker level")
lines(br_s, sfalb2lc,  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_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),]

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


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

par.alb  <- 0.0702  
par.bil <- 0.0856 
X = par.alb * pbcT1$albumin + par.bil *pbcT1$serBilir # X is now the indexed marker
t = par.alb * t.alb + par.bil *t.bil #conditioning value

pbcT1$drug<- X ## store the combine variable in place of 'drug' column which is redundant
## 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   =  par.alb * b.alb + par.bil *b.bil # bandwidth for HQM estimator 
arg1<- get_h_x(pbcT1, 'drug', br_s  =s.out.use, event_time_name = 'years',  time_name = 'year', 
                 event_name = 'status2', t, b.hqm)
br_s2  = seq(Landmark, 14,  length=ls-1) #grid points for HMQ estimator
hqm.surv.est<- make_sf((br_s2[2]-br_s2[1])/5, arg1) # transform HR to Survival func.

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

#Plot the survival functions:
plot(br_s2, hqm.surv.est, 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 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_xll(data, marker_name, br_s, event_time_name, time_name, event_name, x, b)

Arguments

data	A data frame of time dependent data points. Missing values are allowed.
marker_name	The column name of the marker values in the data frame data.
br_s	User defined grid mesh on time_name variable
event_time_name	The column name of the event times in the data frame data.
time_name	The column name of the times the marker values were observed in the data frame data.
event_name	The column name of the events in the data frame data.
x	Numeric value of the last observed marker value.
b	Bandwidth parameter.

Details

The function get_h_xll 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 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_xll 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, h_xt, get_h_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

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


arg1ll<-get_h_xll(pbcT1, 'albumin', br_s = s.out.use, event_time_name = 'years', time_name = 'year',
                  event_name = 'status2', 2, 0.9) 
arg1lc<-get_h_x(pbcT1, 'albumin', br_s = s.out.use, event_time_name = 'years', time_name = 'year',
                event_name = 'status2', 2, 0.9) 

#Calculate the local contant and local linear survival functions
br_s  = seq(Landmark, 14,  length=ls-1)
sfalb2ll<- make_sf((br_s[2]-br_s[1])/4 , arg1ll)
sfalb2lc<- make_sf((br_s[2]-br_s[1])/4 , arg1lc)

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

#Plot the survival functions:
plot(br_s, sfalb2ll,  type="l", col=1, lwd=2, ylab="Survival probability", 
                                                        xlab="Marker level")
lines(br_s, sfalb2lc,  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_xll 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),]
ls<-50 # 50 grid points to evaluate the estimates
s.out<- pbcT1[, 'year']
s.out.use <-  seq(0, max(s.out), max(s.out)/( ls-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)

arg1<- get_h_xll(pbcT1, 'serBilir', br_s= s.out.use,  event_time_name = 'years',  
        time_name = 'year', event_name = 'status2', 1, 10) 
br_s1  = seq(Landmark, 14,  length=ls-1)
sfbil1<- make_sf((br_s1[2]-br_s1[1])/5.4 , arg1)
kma1<- survfit(Surv(years , status2) ~ 1, data = pbcT1)

plot(br_s1, sfbil1, 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_xll 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),]
ls<-50 # 50 grid points to evaluate the estimates
s.out<- pbcT1[, 'year']
s.out.use <-  seq(0, max(s.out), max(s.out)/( ls-1))


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

par.alb  <- 0.0702  
par.bil <- 0.0856 
X = par.alb * pbcT1$albumin + par.bil *pbcT1$serBilir # X is now the indexed marker
t = par.alb * t.alb + par.bil *t.bil #conditioning value

pbcT1$drug<- X ## store X 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   =  par.alb * b.alb + par.bil *b.bil # bandwidth for HQM estimator 
arg1<- get_h_xll(pbcT1, 'drug', br_s = s.out.use, event_time_name = 'years',  time_name = 'year', 
                                               event_name = 'status2', t, b.hqm)
br_s2  = seq(Landmark, 14,  length=ls-1) #grid points for HMQ estimator
hqm.surv.est<- make_sf((br_s2[2]-br_s2[1])/5  ,arg1) # transform HR to Survival func.

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

#Plot the survival functions:
plot(br_s2, hqm.surv.est, 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 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_xt(br_X, br_s, int_X, size_s_grid, alpha, x,t, b, Yi,n)

Arguments

br_X	Vector of grid points for the marker values $X$.
br_s	Vector of grid points for the time values $s$.
int_X	Position of the linear interpolated marker values on the marker grid.
size_s_grid	Size of the time grid.
alpha	Marker-hazard obtained from get_alpha.
x	Numeric value of the last observed marker value.
t	Numeric time value.
b	Bandwidth.
Yi	A matrix made by make_Yi indicating the exposure.
n	Number of individuals.

Details

Function h_xt 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 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

Examples

pbc2_id = to_id(pbc2)
size_s_grid <- size_X_grid <- 100
n = max(as.numeric(pbc2$id))
s = pbc2$year
X = pbc2$serBilir
XX = pbc2_id$serBilir
ss <- pbc2_id$years
delta <- pbc2_id$status2
br_s = seq(0, max(s), max(s)/( size_s_grid-1))
br_X = seq(min(X), max(X), (max(X)-min(X))/( size_X_grid-1))

X_lin = lin_interpolate(br_s, pbc2_id$id, pbc2$id, X, s)

int_X <- findInterval(X_lin, br_X)
int_s = rep(1:length(br_s), n)

N <- make_N(pbc2, pbc2_id, breaks_X=br_X, breaks_s=br_s, ss, XX, delta)
Y <- make_Y(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, int_s, int_X, 
            'years', n)

b = 1.7
alpha<-get_alpha(N, Y, b, br_X, K=Epan )

Yi <- make_Yi(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, int_s, int_X,
              'years', n)

x = 2
t = 2
h_hat = h_xt(br_X, br_s, int_X, size_s_grid, alpha, x, t, b, Yi, n)

---

Hqm estimator on the marker grid

Description

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

Usage

h_xt_vec(br_X, br_s, size_s_grid, alpha, t, b, Yi, int_X, n)

Arguments

br_X	Marker value grid points that will be used in the evaluatiuon.
br_s	Time value grid points that will be used in the evaluatiuon.
size_s_grid	Size of the time grid.
alpha	Marker-hazard obtained from get_alpha.
t	Numeric value of the time the function should be evaluated.
b	Bandwidth.
Yi	A matrix made by make_Yi indicating the exposure.
int_X	Position of the linear interpolated marker values on the marker grid.
n	Number of individuals.

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 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

pbc2_id = to_id(pbc2)
size_s_grid <- size_X_grid <- 100
n = max(as.numeric(pbc2$id))
s = pbc2$year
X = pbc2$serBilir
XX = pbc2_id$serBilir
ss <- pbc2_id$years
delta <- pbc2_id$status2
br_s = seq(0, max(s), max(s)/( size_s_grid-1))
br_X = seq(min(X), max(X), (max(X)-min(X))/( size_X_grid-1))

X_lin = lin_interpolate(br_s, pbc2_id$id, pbc2$id, X, s)

int_X <- findInterval(X_lin, br_X)
int_s = rep(1:length(br_s), n)

N <- make_N(pbc2, pbc2_id, br_X, br_s, ss, XX, delta)
Y <- make_Y(pbc2, pbc2_id, X_lin, br_X, br_s,
            size_s_grid, size_X_grid, int_s, int_X, event_time = 'years', n)

b = 1.7
alpha<-get_alpha(N, Y, b, br_X, K=Epan )

Yi <- make_Yi(pbc2, pbc2_id, X_lin, br_X, br_s,
              size_s_grid, size_X_grid, int_s, int_X, event_time = 'years', n)

t = 2

h_xt_vec(br_X, br_s, size_s_grid, alpha, t, b, Yi, int_X, 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
arg1<- get_h_x(pbcT1, 'albumin', br_s=br_s,  event_time_name = 'years',  time_name = 'year', 
                                                    event_name = 'status2', 2, b) 
br_s2  = seq(Landmark, 14,  length=99)
sfalb2<- make_sf(    (br_s2[2]-br_s2[1])/1.35 , arg1)
kma2<- survfit(Surv(years , status2) ~ 1, data = pbcT1)

#Plot the survival functions:
plot(br_s2, sfalb2, 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
x.grid<-br_s[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_list = seq(1.5, 3, 0.05)
 # b_scores_alb = b_selection(pbc2, 'albumin', 'years', 'year', 'status2', I, b_list)
 # b.alb =  b_scores_alb[[2]][which.min(b_scores_alb[[1]])]
 
b.bil = 4 # results from CV rule of Bagkavos et al. (2025) : 
  # b_list.bil = seq(3, 4, 0.05)
  # b_scores_bil = b_selection(pbc2, 'serBilir', 'years', 'year', 'status2', I, b_list.bil)
  # b.bil =  b_scores_alb[[2]][which.min(b_scores_alb[[1]])] - result = 4

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

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

# First, mean adjust the albumin covariate
Xalb = pbc2$albumin - mean(pbc2$albumin)
XXalb = pbc2_id$albumin - mean(pbc2_id$albumin)

br_Xalb = seq(min(Xalb), max(Xalb), (max(Xalb)-min(Xalb))/( size_X_grid-1))
X_linalb = lin_interpolate(br_s, pbc2_id$id, pbc2$id, Xalb, s)
int_Xalb <- findInterval(X_linalb, br_Xalb)
N.alb <- make_N(pbc2, pbc2_id, br_Xalb, br_s, ss, XXalb, delta)
Y.alb <- make_Y(pbc2, pbc2_id, X_linalb, br_Xalb, br_s,
            size_s_grid, size_X_grid, int_s, int_Xalb, event_time = 'years', n)

alpha.alb<-get_alpha(N.alb, Y.alb, b.alb, br_Xalb, K=Epan )
Yi.alb <- make_Yi(pbc2, pbc2_id, X_linalb, br_Xalb, br_s, size_s_grid, size_X_grid, 
                                          int_s, int_Xalb, event_time = 'years', n)

#calculate HQM estimator for original albumin covariate:
arg.alb<- h_xt_vec(br_Xalb, br_s, size_s_grid, alpha.alb, t.alb, b.alb, Yi.alb, 
                                                                  int_Xalb, n)

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

# Now, mean adjust the bilirubon covariate:
Xbil = pbc2$serBilir - mean(pbc2$serBilir)
XXbil = pbc2_id$serBilir - mean(pbc2_id$serBilir)

br_Xbil = seq(min(Xbil), max(Xbil), (max(Xbil)-min(Xbil))/( size_X_grid-1))
X_linbil = lin_interpolate(br_s, pbc2_id$id, pbc2$id, Xbil, s)
int_Xbil <- findInterval(X_linbil, br_Xbil)
N.bil <- make_N(pbc2, pbc2_id, br_Xbil, br_s, ss, XXbil, delta)
Y.bil <- make_Y(pbc2, pbc2_id, X_linbil, br_Xbil, br_s,
            size_s_grid, size_X_grid, int_s, int_Xbil, event_time = 'years', n)

alpha.bil<-get_alpha(N.bil, Y.bil, b.bil, br_Xbil, K=Epan )
Yibil <- make_Yi(pbc2, pbc2_id, X_linbil, br_Xbil, br_s, size_s_grid, size_X_grid, 
                                        int_s, int_Xbil, event_time = 'years', n)

#calculate HQM estimator for original bilirubin covariate:
arg1.bil<- h_xt_vec(br_Xbil, br_s, size_s_grid, alpha.bil, t.bil, b.bil, Yibil, 
                                                                    int_Xbil, n)
y.grid.bil<-arg1.bil[1:index.use]

# calculate the indexed variables:
X = par.alb * Xalb + 	par.bil *Xbil
XX = par.alb * XXalb + par.bil *XXbil
b = 0.4 # results from CV rule in Bagkavos et al (2025) 
t = par.alb * t.alb + par.bil *t.bil

br_X = seq(min(X), max(X), (max(X)-min(X))/( size_X_grid-1))

X_lin = lin_interpolate(br_s, pbc2_id$id, pbc2$id, X, s)
int_X <- findInterval(X_lin, br_X)
N <- make_N(pbc2, pbc2_id, breaks_X=br_X, breaks_s=br_s, ss, XX, delta)
Y <- make_Y(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, int_s, 
                                                             int_X, 'years', n)
alpha<-get_alpha(N, Y, b, br_X, K=Epan )
Yi <- make_Yi(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, 
                                                      int_s, int_X,'years', n)

# Now calculate indexed HR estimator:
arg2<- h_xt_vec(br_X, br_s, size_s_grid, alpha, t, b, Yi, int_X, n) 
y.grid2<-arg2[1:index.use] 
 

############## Plot the results:
plot(x.grid, y.grid2, type="l", ylim=c(0,2), xlab="Years", ylab="Hazard rate", lwd=2 )
lines(x.grid, y.grid.bil , lty=2, col=4, lwd=2)
lines(x.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
x.grid<-br_s[1:index.use]

b.pro = 3
b.bil = 4

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

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

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

ind2 <- which( is.na( pbc2_id$prothrombin ) )
XXpro1 =   pbc2_id$prothrombin   
XXpro1[ind2]<-mean(pbc2_id$prothrombin, na.rm=TRUE)
XXpro1 <- XXpro1 - mean(XXpro1)

br_Xbil = seq(min(Xpro1), max(Xpro1), (max(Xpro1)-min(Xpro1))/( size_X_grid-1))
X_linbil = lin_interpolate(br_s, pbc2_id$id, pbc2$id, Xpro1, s)
int_Xbil <- findInterval(X_linbil, br_Xbil)
N.bil <- make_N(pbc2, pbc2_id, br_Xbil, br_s, ss, XXpro1, delta)
Y.bil <- make_Y(pbc2, pbc2_id, X_linbil, br_Xbil, br_s,
                size_s_grid, size_X_grid, int_s, int_Xbil, event_time = 'years', n)

alpha.bil<-get_alpha(N.bil, Y.bil, b.pro, br_Xbil, K=Epan )
Yibil <- make_Yi(pbc2, pbc2_id, X_linbil, br_Xbil, br_s, size_s_grid, size_X_grid, int_s, 
                                                          int_Xbil, event_time = 'years', n)
arg1.pro<- h_xt_vec(br_Xbil, br_s, size_s_grid, alpha.bil, t.pro, b.pro, Yibil, int_Xbil, n)
y.grid.pro<-arg1.pro[1:index.use]


############ Bilirubin #############
Xbil = pbc2$serBilir - mean(pbc2$serBilir)
XXbil = pbc2_id$serBilir - mean(pbc2_id$serBilir)

br_Xbil = seq(min(Xbil), max(Xbil), (max(Xbil)-min(Xbil))/( size_X_grid-1))
X_linbil = lin_interpolate(br_s, pbc2_id$id, pbc2$id, Xbil, s)
int_Xbil <- findInterval(X_linbil, br_Xbil)
N.bil <- make_N(pbc2, pbc2_id, br_Xbil, br_s, ss, XXbil, delta)
Y.bil <- make_Y(pbc2, pbc2_id, X_linbil, br_Xbil, br_s,
            size_s_grid, size_X_grid, int_s, int_Xbil, event_time = 'years', n)

alpha.bil<-get_alpha(N.bil, Y.bil, b.bil, br_Xbil, K=Epan )
Yibil <- make_Yi(pbc2, pbc2_id, X_linbil, br_Xbil, br_s, size_s_grid, size_X_grid, int_s, 
                                                          int_Xbil, event_time = 'years', n)
arg1.bil<- h_xt_vec(br_Xbil, br_s, size_s_grid, alpha.bil, t.bil, b.bil, Yibil, int_Xbil, n)
y.grid.bil<-arg1.bil[1:index.use]

######### Index prothrombin time and bilirubin
X = par.pro * Xpro1 + 	par.bil *Xbil
XX = par.pro * XXpro1 + par.bil *XXbil

b = 1.5  
t = par.pro * t.pro + par.bil *t.bil

br_X = seq(min(X), max(X), (max(X)-min(X))/( size_X_grid-1))

X_lin = lin_interpolate(br_s, pbc2_id$id, pbc2$id, X, s)
int_X <- findInterval(X_lin, br_X)
N <- make_N(pbc2, pbc2_id, breaks_X=br_X, breaks_s=br_s, ss, XX, delta)
Y <- make_Y(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, int_s, int_X, 
                                                                          'years', n)
alpha<-get_alpha(N, Y, b, br_X, K=Epan )
Yi <- make_Yi(pbc2, pbc2_id, X_lin, br_X, br_s, size_s_grid, size_X_grid, int_s, int_X,
                                                                          'years', n)
arg2<- h_xt_vec(br_X, br_s, size_s_grid, alpha, t, b, Yi, int_X, n)
y.grid2<-arg2[1:index.use]
 


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

lines(x.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")

---