View source: R/NPICbayesSurv.R
NPICbayesSurv | R Documentation |
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 |
upp |
upper limits of observed intervals with |
choice |
a character string indicating the initial guess (S^*) of the survival distribution |
cc |
parameter of the Dirichlet process prior |
parm |
a numeric vector of parameters for the initial guess:
|
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
a 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)
grid of time points (excluding an “infinity” value)
a matrix with sampled weights
a matrix with sampled values of n
a matrix with sampled values of the survival function
parameter of the Dirichlet process prior which was used
character indicating the initial guess
parameters of the initial guess
length of the burn-in
number of sampled values
Emmanuel Lesaffre emmanuel.lesaffre@kuleuven.be, Arnošt Komárek arnost.komarek@mff.cuni.cz
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.
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.