NPICbayesSurv {icensBKL}

Bayesian non-parametric estimation of a survival curve with interval-censored data

Description

Bayesian non-parametric estimation of a survival curve for right-censored data as proposed by Susarla and Van Ryzin (1976, 1978)

Usage

NPICbayesSurv(low, upp, 
  choice = c("exp", "weibull", "lnorm"), cc, parm,
  n.sample = 5000, n.burn = 5000, cred.level = 0.95)

Arguments

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 (S^*) 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

Value

A list with the following components

S

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)

t

grid of time points (excluding an “infinity” value)

w

a matrix with sampled weights

n

a matrix with sampled values of n

Ssample

a matrix with sampled values of the survival function

c

parameter of the Dirichlet process prior which was used

choice

character indicating the initial guess

parm

parameters of the initial guess

n.burn

length of the burn-in

n.sample

number of sampled values

Author(s)

Emmanuel Lesaffre [e-mail], Arnošt Komárek [e-mail]

References

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.

Examples

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