Bayesian non-parametric estimation of a survival curve with interval-censored data
Bayesian non-parametric estimation of a survival curve for right-censored data as proposed by Susarla and Van Ryzin (1976, 1978)
NPICbayesSurv(low, upp, choice = c("exp", "weibull", "lnorm"), cc, parm, n.sample = 5000, n.burn = 5000, cred.level = 0.95)
low: lower limits of observed intervals with NA if left-censored
upp: upper limits of observed intervals with NA if right-censored
choice: a character string indicating the initial guess () of the survival distribution
cc: parameter of the Dirichlet process prior
parm: a numeric vector of parameters for the initial guess: rate parameter for exponential (see also Exponential), a two-element vector with shape
and scale parameters for weibull (see also Weibull), a two-element vector with meanlog and sdlog parameters for log-normal (see also Lognormal). If not given, parameters for the initial guess are taken from the ML fit
n.sample: number of iterations of the Gibbs sampler after the burn-in
n.burn: length of the burn-in
cred.level: credibility level of calculated pointwise credible intervals for values of the survival function
A list with the following components
data.frame with columns: t (time points), S (posterior mean of the value of the survival function at t), Lower, Upper (lower and upper bound of the pointwise credible interval for the value of the survival function)Calle, M. L. and Gómez, G. (2001). Nonparametric Bayesian estimation from interval-censored data using Monte Carlo methods. Journal of Statistical Planning and Inference, 98 (1-2), 73-87.
Emmanuel Lesaffre emmanuel.lesaffre@kuleuven.be , Arnošt Komárek arnost.komarek@mff.cuni.cz
## Breast Cancer study (radiotherapy only group) ## Dirichlet process approach to estimate nonparametrically ## the survival distribution with interval-censored data data("breastCancer", package = "icensBKL") breastR <- subset(breastCancer, treat == "radio only", select = c("low", "upp")) ### Lower and upper interval limit to be used here low <- breastR[, "low"] upp <- breastR[, "upp"] ### Common parameters for sampling ### (quite low, only for testing) n.sample <- 100 n.burn <- 100 ### Gibbs sampler set.seed(19680821) Samp <- NPICbayesSurv(low, upp, choice = "weibull", n.sample = n.sample, n.burn = n.burn) print(ncol(Samp$w)) ## number of supporting intervals print(nrow(Samp$S)) ## number of grid points (without "infinity") print(Samp$S[, "t"]) ## grid points (without "infinity") print(Samp$t) ## grid points (without "infinity") print(Samp$S) ## posterior mean and pointwise credible intervals print(Samp$w[1:10,]) ## sampled weights (the first 10 iterations) print(Samp$n[1:10,]) ## sampled latend vectors (the first 10 iterations) print(Samp$Ssample[1:10,]) ## sampled S (the first 10 iterations) print(Samp$parm) ## parameters of the guess ### Fitted survival function including pointwise credible bands ngrid <- nrow(Samp$S) plot(Samp$S[1:(ngrid-1), "t"], Samp$S[1:(ngrid-1), "Mean"], type = "l", xlim = c(0, 50), ylim = c(0, 1), xlab = "Time", ylab = expression(hat(S)(t))) polygon(c(Samp$S[1:(ngrid-1), "t"], Samp$S[(ngrid-1):1, "t"]), c(Samp$S[1:(ngrid-1), "Lower"], Samp$S[(ngrid-1):1, "Upper"]), col = "grey95", border = NA) lines(Samp$S[1:(ngrid - 1), "t"], Samp$S[1:(ngrid - 1), "Lower"], col = "grey", lwd = 2) lines(Samp$S[1:(ngrid - 1), "t"], Samp$S[1:(ngrid - 1), "Upper"], col = "grey", lwd = 2) lines(Samp$S[1:(ngrid - 1), "t"], Samp$S[1:(ngrid - 1), "Mean"], col = "black", lwd = 3)