> On Dec 28, 2016, at 4:30 AM, radhika sundar <radhikasundar...@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 <zh...@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 <dwinsem...@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 <radhikasundar...@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@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 ______________________________________________ R-help@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.
