Pivot.Index.CIs: 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

Pivot.Index.CIs(B,  n.est.points, Mat.boot.haz.rate, time.grid, hqm.est, a.sig)

Value

A data frame with the following columns:

time: The evaluation grid points.
est: Indexed hazard rate estimator hqm.est.
downci, upci: Lower and upper endpoints of basic pivot CIs.
docisym, upcisym: Lower and upper endpoints of symmetric pivot CIs.
logdoci, logupci: Lower and upper endpoints of pivot CIs on the log-scale.
logdocisym, logupcisym: Symmetric pivot CIs on the log-scale.
log.est: The logarithm of the indexed hazard rate estimate, log(hqm.est).
tLogDoCI, tLogUpCI: Transformed-log CIs based on 2*log(hqm.est) - log-quantiles.
tSymLogDoCI, tSymLogUpCI: Symmetric transformed-log CIs.

Arguments

B: Integer. Number of bootsrap replications
n.est.points: Integer. Number of estimation points at which the indexed hazard estimates and confidence intervals are evaluated.
Mat.boot.haz.rate: A matrix of bootstrap estimated hazard rates with dimensions n.est.points × B, where each column corresponds to one bootstrap replicate.
time.grid: Numeric vector of length n.est.points: the grid points at which the indexed hazard estimates and confidence intervals are calculated.
hqm.est: Indexed hazard estimator, calculated at the grid points time.grid and using the original sample.
a.sig: The significance level (e.g., 0.05) which will be used in computing the confidence intervals.

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 Mat.boot.haz.rate is assumed to contain estimates produced using the same time grid as time.grid and the same estimator used to generate hqm.est.

Examples

Run this code

# \donttest{ 
marker_name1 <- 'albumin'
marker_name2 <-  'serBilir'
event_time_name <- 'years' 
time_name <- 'year' 
event_name <- 'status2'
id<-'id'



par.x1  <- 0.0702 #0.149
par.x2 <- 0.0856 #0.10
t.x1 = 0 # refers to zero mean variables - slightly high
t.x2 = 1.9 # refers to zero mean variable - high
b = 0.42#par.alb * b.alb + par.bil *b.bil # 7
t = par.x1 * t.x1 + par.x2 *t.x2
ls<-50

#Store original sample values:
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)

X1t=xin[,marker_name1] -mean(xin[, marker_name1])
XX1t=xin.id[,marker_name1] -mean(xin.id[, marker_name1])
X2t=xin[,marker_name2]  -mean(xin[, marker_name2])
XX2t=xin.id[,marker_name2] -mean(xin.id[, marker_name2])

X1=list(X1t, X2t)
XX1=list(XX1t, XX2t)

# Calculate the indexed HQM estimator on the original sample:
arg2<- SingleIndCondFutHaz(pbc2, id, ls,  X1, XX1, event_time_name = 'years', 
        time_name = 'year',  event_name = 'status2', in.par= c(par.x1,  par.x2), b, t)
                          
hqm.est<-arg2[,2] # Indexed HQM estimator on original sample
time.grid<-arg2[,1] # evaluation grid points
n.est.points<- ls # length(hqm.est)

#  Create bootstrap samples by group 
set.seed(1)  
B<- 50  #for display purposes only; for sensible results use B=1000 (slower) 
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),]  
}     

# Simulate true hazard rate function:   
true.hazard<- Sim.True.Hazard(Boot.samples, id='id', n.est.points, 
              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)

# Bootstrap the original indexed HQM estimator: 
res <- Boot.hrandindex.param(B, Boot.samples, marker_name1, marker_name2, 
                  event_time_name, time_name,  event_name, b = 0.4, t = t, 
                  true.haz = true.hazard, v.param = c(0.07, 0.08), n.est.points)

# Construct Ci's: 
a.sig<-0.05   
all.cis.pivot<- Pivot.Index.CIs(B, n.est.points, res, time.grid, hqm.est, a.sig)

# extract Plain   + symmetric CIs 
UpCI<-all.cis.pivot[,"UpCI"]
DoCI<-all.cis.pivot[,"DoCI"]
SymUpCI<-all.cis.pivot[,"SymUpCI"]
SymDoCI<-all.cis.pivot[,"SymDoCI"]

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

# extract transformed from Log HR + symmetric CIs 
LogUpCI<-all.cis.pivot[,"LogUpCI"]
LogDoCI<-all.cis.pivot[,"LogDoCI"]
SymLogUpCI<-all.cis.pivot[,"LogSymUpCI"]
SymLogDoCI<-all.cis.pivot[,"LogSymDoCI"]

#Plot the selected CIs
plot(time.grid[1:J],   hqm.est[1:J] , type="l", ylim=c(0,2), ylab="Log Hazard rate", 
      xlab="time", lwd=2)
polygon(x = c(time.grid[1:J], rev(time.grid[1:J])), y = c(LogUpCI[1:J], rev(LogDoCI[1:J])), 
        col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(time.grid[1:J], SymLogUpCI[1:J], lty=2, lwd=2 ) 
lines(time.grid[1:J],  SymLogDoCI[1:J], lty=2, lwd=2)  

# extract Log HR + symmetric CIs  
tLogUpCI<-all.cis.pivot[,"LogtUpCI"]
tLogDoCI<-all.cis.pivot[,"LogTDoCI"]
tSymLogUpCI<-all.cis.pivot[,"SymLogtUpCI"]
tSymLogDoCI<-all.cis.pivot[,"SymLogTDoCI"]

#Plot the selected CIs
plot(time.grid[1:J], log(hqm.est[1:J] ), type="l", ylim=c(-5,4), ylab="Hazard rate", 
      xlab="time",  lwd=2)
polygon(x = c(time.grid[1:J], rev(time.grid[1:J])), y = c(tLogUpCI[1:J], rev(tLogDoCI[1:J])), 
        col =  adjustcolor("red", alpha.f = 0.50), border = NA )
lines(time.grid[1:J], tSymLogUpCI[1:J], lty=2, lwd=2 ) 
lines(time.grid[1:J],  tSymLogDoCI[1:J], lty=2, lwd=2)  
# }