Bayesian survival analysis

Question

I am going through R's function indeptCoxph() in the spBayesSurv package which fits a bayesian Cox model. I am confused by some of the input parameters to this functions. (I also had some questions about the R code which I have posted separately on Stack Overflow: Stuck with package example code in R - simulating data to fit a model).

What is the role of the "prediction" input parameter? Should it not only contain the predictor covariates? In the R example, the authors have included a vector "s" which was used to initially simulate the survival times data in their example as well as the predictors. I'm not sure what this "s" is.

Given that my data is just a set of survival times between 0 and 100, along with censored (yes/no) information, how would I use this function and how should I handle the input "s"?

(I have also posted on SO, but posting here too since I would like to understand the theory behind this model ).

Answer 1

The function example is conducted under the framework of spatial copula models (i.e. the function spCopulaCoxph). Briefly speaking, you just need to ignore the spred=s0 in the prediction settings, that is, prediction=list(xpred=xpred) is sufficient. To be more clear, a new example is attached at the end.

Alternatively, the newly developed function survregbayes (https://rdrr.io/cran/spBayesSurv/man/survregbayes.html) is more user-friendly to use, which fits three popular semiparametric survival models (either non-, iid-, CAR-, or GRF-frailties): proportional hazards, accelerated failure time, and proportional odds.

    ###############################################################
    # A simulated data: Cox PH
    ###############################################################
    rm(list=ls())
    library(survival)
    library(spBayesSurv)
    library(coda)
    library(MASS)
    ## True parameters 
    betaT = c(1,1); 
    n=500; npred=30; ntot=n+npred;
    ## Baseline Survival
    f0oft = function(t) 0.5*dlnorm(t, -1, 0.5)+0.5*dlnorm(t,1,0.5);
    S0oft = function(t) (0.5*plnorm(t, -1, 0.5, lower.tail=FALSE)+
                           0.5*plnorm(t, 1, 0.5, lower.tail=FALSE))
    ## The Survival function:
    Sioft = function(t,x)  exp( log(S0oft(t))*exp(sum(x*betaT)) ) ;
    fioft = function(t,x) exp(sum(x*betaT))*f0oft(t)/S0oft(t)*Sioft(t,x);
    Fioft = function(t,x) 1-Sioft(t,x);
    ## The inverse for Fioft
    Finv = function(u, x) uniroot(function (t) Fioft(t,x)-u, lower=1e-100, 
                                  upper=1e100, extendInt ="yes", tol=1e-6)$root

    ## generate x 
    x1 = rbinom(ntot, 1, 0.5); x2 = rnorm(ntot, 0, 1); X = cbind(x1, x2);
    ## generate survival times
    u = runif(ntot);
    tT = rep(0, ntot);
    for (i in 1:ntot){
      tT[i] = Finv(u[i], X[i,]);
    }

    ## right censoring
    t_obs=tT 
    Centime = runif(ntot, 2, 6);
    delta = (tT<=Centime) +0 ; 
    length(which(delta==0))/ntot; # censoring rate
    rcen = which(delta==0);
    t_obs[rcen] = Centime[rcen]; ## observed time 
    ## make a data frame
    dtotal = data.frame(t_obs=t_obs, x1=x1, x2=x2, delta=delta, 
                        tT=tT);
    ## Hold out npred=30 for prediction purpose
    predindex = sample(1:ntot, npred);
    dpred = dtotal[predindex,];
    dtrain = dtotal[-predindex,];

    # Prediction settings 
    xpred = cbind(dpred$x1,dpred$x2);
    prediction = list(xpred=xpred);

    ###############################################################
    # Independent Cox PH
    ###############################################################
    # MCMC parameters
    nburn=1000; nsave=1000; nskip=0;
    # Note larger nburn, nsave and nskip should be used in practice.
    mcmc=list(nburn=nburn, nsave=nsave, nskip=nskip, ndisplay=1000);
    prior = list(M=10);
    state <- NULL;
    # Fit the Cox PH model
    res1 = indeptCoxph( y = dtrain$t_obs, delta =dtrain$delta, 
                        x = cbind(dtrain$x1, dtrain$x2),RandomIntervals=FALSE, 
                        prediction=prediction,  prior=prior, mcmc=mcmc, state=state);
    save.beta = res1$beta; row.names(save.beta)=c("x1","x2")
    apply(save.beta, 1, mean); # coefficient estimates
    apply(save.beta, 1, sd); # standard errors
    apply(save.beta, 1, function(x) quantile(x, probs=c(0.025, 0.975))) # 95% CI
    ## traceplot
    par(mfrow = c(2,1))
    traceplot(mcmc(save.beta[1,]), main="beta1")
    traceplot(mcmc(save.beta[2,]), main="beta2")
    res1$ratebeta; # adaptive MH acceptance rate
    ## LPML
    LPML1 = sum(log(res1$cpo)); LPML1;
    ## MSPE
    mean((dpred$tT-apply(res1$Tpred, 1, median))^2); 

    ## plots
    par(mfrow = c(2,1))
    x1new = c(0, 0);
    x2new = c(0, 1)
    xpred = cbind(x1new, x2new); 
    nxpred = nrow(xpred);
    tgrid = seq(1e-10, 4, 0.03);
    ngrid = length(tgrid);
    estimates = GetCurves(res1, xpred, log(tgrid), CI=c(0.05, 0.95));
    fhat = estimates$fhat; 
    Shat = estimates$Shat;
    ## density in t
    plot(tgrid, fioft(tgrid, xpred[1,]), "l", lwd=2,  ylim=c(0,3), main="density")
    for(i in 1:nxpred){
      lines(tgrid, fioft(tgrid, xpred[i,]), lwd=2)
      lines(tgrid, fhat[,i], lty=2, lwd=2, col=4);
    }
    ## survival in t
    plot(tgrid, Sioft(tgrid, xpred[1,]), "l", lwd=2, ylim=c(0,1), main="survival")
    for(i in 1:nxpred){
      lines(tgrid, Sioft(tgrid, xpred[i,]), lwd=2)
      lines(tgrid, Shat[,i], lty=2, lwd=2, col=4);
      lines(tgrid, estimates$Shatup[,i], lty=2, lwd=1, col=4);
      lines(tgrid, estimates$Shatlow[,i], lty=2, lwd=1, col=4);
    }