Skip to content

Bayesian cox model: spBayesSurv package

4 messages · David Winsemius, radhika sundar

Original
radhika sundar
# 
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 function.

What is the role of the "prediction" input parameter? Should it not only
contain the predictor covariates? In the R code 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"?


Thanks for any help!
David Winsemius
# 
Cross posting is deprecated on rhelp but if you do so, please at least post a link to the stackoverflow address for the duplicate question.

Sent from my iPhone
radhika sundar
# 
The example code for the indeptCoxph function  in the spBayesSurv package
has been updated(this is not cross-posted anywhere else). See code below.
The author simulates data to illustrate the Cox model -  I am stuck trying
to understand the role of the functions f0oft, S0oft, fioft and Sioft as
also Finv.

Am I right in thinking that he is only using the above mentioned functions
to simulate his data? If I wanted to run the indeptCoxph function with
different data, do I need to define those functions again?
Roughly, in the author's example, I can understand he fits a Cox model with
2 predictors x1 and x2. He simulates survival time data (but this is where
I am confused).
As for the bayesian model itself, the only prior he uses is for M, the
number of cutpoints in the baseline hazard function. (There is no function
listed as a prior for Survival times in the ideptCoxph call).

Sorry --- this is a novice question relating to understanding both the
statistical set-up and the R-code.
Thanks for any help!
Author's(updated)  code:
###############################################################
# 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);
}


On Wed, Dec 28, 2016 at 2:11 AM, David Winsemius <dwinsemius at comcast.net>
wrote:

  
  


  


      

    

    

      
      

        


  

        

      David Winsemius
    

    

      
      #
    

  


  

          

      
      
      
    

        
Your preamble very much resembles the 2 questions from yesterday of the anonymous questioner: user2450223

http://stackoverflow.com/questions/41344300/bayesian-survival-analysis
http://stackoverflow.com/questions/41342836/stuck-with-package-example-code-in-r-simulating-data-to-fit-a-model

One of those questions was apparently answered by the package's author yesterday, although you (at least I think it must have been you) have not acknowledged it yet. If you are asking question on Rhelp about little-used packages you are advised in the Posting Guide to first contact the package author and if unsuccessful, then post to Rhelp (preferably with visible CC: to the author as well as to the list.)

Learn the `maintainer` function:

            

      
      
      
    

        
[1] "Haiming Zhou <zhouh at niu.edu>"

            

      
      
      
    

        
Rhelp is not set up to handle questions that are fundamentally statistical. When the issue is the underlying statistical theory for running R code the best place start would be the package author, and only when unsuccessful post to whatever alternate location he suggests in the packageDescription (although I don't see one) or then post to:

http://stats.stackexchange.com/

Given the package's "spatial" capacities, it might also have been appropriate on:

https://stat.ethz.ch/mailman/listinfo/r-sig-geo


Neither R-help nor StackOverflow are good fits to this question in my opinion.