hdslc: Hazard discrimination summary estimator
In liangcj/hds: Hazard Discrimination Summary

Description Usage Arguments Details Value References See Also Examples

View source: R/hds_lc_hat.R

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

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

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

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.

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

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

hds, finda

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