[R] Bayesian cox model: spBayesSurv package

Wed Dec 28 13:30:34 CET 2016

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:

> 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
>
> > On Dec 27, 2016, at 6:52 AM, radhika sundar <radhikasundar520 at gmail.com>
> wrote:
> >
> > 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!
> >
> >    [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/
> posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
>

	[[alternative HTML version deleted]]