[R] Bayesian cox model: spBayesSurv package
David Winsemius
dwinsemius at comcast.net
Wed Dec 28 18:07:32 CET 2016
> On Dec 28, 2016, at 4:30 AM, radhika sundar <radhikasundar520 at gmail.com> wrote:
>
> 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.
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:
> maintainer( 'spBayesSurv')
[1] "Haiming Zhou <zhouh at niu.edu>"
>
> 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).
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.
--
David.
>
> 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.
>
>
David Winsemius
Alameda, CA, USA
More information about the R-help
mailing list