hdslc {hds}

Hazard discrimination summary estimator

Description

Returns local constant HDS estimates at all specified evaluation times. See Liang and Heagerty (2016) for details on HDS.

Usage

hdslc(times, status, m, evaltimes = times[order(times)], h = 1.06 *
  sd(times) * (length(times)^(-0.2)), se = TRUE)

Arguments

times

A vector of observed follow up times.

status

A vector of status indicators, usually 0=alive, 1=dead.

m

A matrix or data frame of covariate values, with a column for each covariate and each observation is on a separate row. Non-numeric values are acceptable, as the values will be transformed into a numeric model matrix through survival::coxph.

evaltimes

A vector of times at which to estimate HDS. Defaults to all the times specified by the times vector. If there are a lot of observations, then you may want to enter in a sparser vector of times for faster computation.

h

A single numeric value representing the bandwdith to use, on the time scale. The default bandwidth is a very ad hoc estimate using Silverman's rule of thumb

se

TRUE or FALSE. TRUE: calculate and return standard error estimates. FALSE: do not calculate standard errors estimates and return NAs. Defaults to TRUE. May want to set to FALSE to save computation time if using this function to compute bootstrap standard errors.

Details

A local constant version of hds. While hds estimates HDS(t) assuming the Cox proportional hazards model, hdslc estimates HDS(t) using a relaxed, local-in-time Cox model. Specifically, the hazard ratios are allowed to vary with time. See Cai and Sun (2003) and Tian Zucker Wei (2005) for details on the local-in-time Cox model.

Point estimates use hdslc.fast and standard errors use hdslcse.fast. hdslc.fast requires an estimate of beta(t) (time-varying hazard ratio), which is estimated using finda(); and subject specific survival, which is estimated using sssf.fast(). hdslcse.fast requires the same and in addition standard error estimates of beta(t), which are estimated using betahatse.fast().

The covariate values m are centered for numerical stability. This is particularly relevant for the standard error calculations.

Value

A data frame with three columns: 1) the evaluation times, 2) the HDS estimates at each evaluation time, and 3) the standard error estimates at each evaluation time

References

Liang CJ and Heagerty PJ (2016). A risk-based measure of time-varying prognostic discrimination for survival models. Biometrics. doi:10.1111/biom.12628

Cai Z and Sun Y (2003). Local linear estimation for time-dependent coefficients in Cox's regression models. Scandinavian Journal of Statistics, 30: 93-111. doi:10.1111/1467-9469.00320

Tian L, Zucker D, and Wei LJ (2005). On the Cox model with time-varying regression coefficients. Journal of the American Statistical Association, 100(469):172-83. doi:10.1198/016214504000000845

See Also

hds, finda

Examples

## Not run: 
head(hdslc(times = survival::pbc[1:312, 2],
           status = (survival::pbc[1:312, 3]==2)*1,
           m = survival::pbc[1:312, 5]))

hdsres   <- hds(times=pbc5[,1], status=pbc5[,2], m=pbc5[,3:7])
hdslcres <- hdslc(times = pbc5[,1], status=pbc5[,2], m = pbc5[,3:7], h = 730)
Survt    <- summary(survival::survfit(survival::Surv(pbc5[,1], pbc5[,2])~1))
Survtd   <- cbind(Survt$time, c(0,diff(1-Survt$surv)))
tden     <- density(x=Survtd[,1], weights=Survtd[,2], bw=100, kernel="epanechnikov")

par(mar=c(2.25,2.25,0,0)+0.1, mgp=c(1.25,0.5,0))
plot(c(hdslcres[,1], hdslcres[,1]), c(hdslcres[,2] - 1.96*hdslcres[,3],
                                      hdslcres[,2] + 1.96*hdslcres[,3]),
     type="n", xlab="days", ylab="HDS(t)", cex.lab=.75, cex.axis=.75,
     ylim=c(-3,15), xlim=c(0,3650))
polygon(x=c(hdsres[,1], hdsres[312:1,1]), col=rgb(1,0,0,.25), border=NA,
        fillOddEven=TRUE, y=c(hdsres[,2]+1.96*hdsres[,3],
                              (hdsres[,2]-1.96*hdsres[,3])[312:1]))
polygon(x=c(hdslcres[,1], hdslcres[312:1, 1]), col=rgb(0,0,1,.25), border=NA,
        fillOddEven=TRUE, y=c(hdslcres[,2] + 1.96*hdslcres[,3],
                              (hdslcres[,2] - 1.96*hdslcres[,3])[312:1]))
lines(hdsres[,1], hdsres[,2], lwd=2, col="red")
lines(hdslcres[,1], hdslcres[,2], lwd=2, col="blue")
abline(h=1, lty=3)
legend(x=1200, y=14, legend=c("Proportional hazards",
                              "Local-in-time proportional hazards",
                              "Time density"), col=c("red", "blue", "gray"),
       lwd=2, bty="n", cex=0.75)
with(tden, polygon(c(x, x[length(x):1]),
                   c(y*3/max(y)-3.5, rep(-3.5, length(x))),
                   col="gray", border=NA, fillOddEven=TRUE))

## End(Not run)