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]]