Description Usage Arguments Details Value Author(s) References See Also Examples
This function fits a Bayesian proportional hazards model for non-spatial right censored time-to-event data.
1 2 3 |
formula |
a formula expression with the response returned by the |
data |
a data frame in which to interpret the variables named in the |
na.action |
a missing-data filter function, applied to the |
prediction |
a list giving the information used to obtain conditional inferences. The list includes the following element: |
mcmc |
a list giving the MCMC parameters. The list must include the following elements: |
prior |
a list giving the prior information. The list includes the following elements: |
state |
a list giving the current value of the parameters. This list is used if the current analysis is the continuation of a previous analysis. |
scale.designX |
flag to indicate wheter the design matrix X will be centered by column means and scaled by column standard deviations, where |
This function fits a Bayesian proportional hazards model (Zhou, Hanson and Zhang, 2017) for non-spatial right censored time-to-event data.
The results include the MCMC chains for the parameters discussed in Zhou, Hanson and Zhang (2017).
Use names
to find out what they are.
Haiming Zhou <[email protected]> and Tim Hanson <[email protected]>
Zhou, H., Hanson, T., and Knapp, R. (2015). Marginal Bayesian nonparametric model for time to disease arrival of threatened amphibian populations. Biometrics, 71(4): 1101-1110.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 | ###############################################################
# A simulated data: Cox PH
###############################################################
rm(list=ls())
library(survival)
library(spBayesSurv)
library(coda)
## True parameters
betaT = c(1,1);
n=100;
## 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(n, 1, 0.5); x2 = rnorm(n, 0, 1); X = cbind(x1, x2);
## generate survival times
u = runif(n);
tT = rep(0, n);
for (i in 1:n){
tT[i] = Finv(u[i], X[i,]);
}
### ----------- right-censored -------------###
t_obs=tT
Centime = runif(n, 2, 6);
delta = (tT<=Centime) +0 ;
length(which(delta==0))/n; # censoring rate
rcen = which(delta==0);
t_obs[rcen] = Centime[rcen]; ## observed time
## make a data frame
d = data.frame(tobs=t_obs, x1=x1, x2=x2, delta=delta, tT=tT);
table(d$delta)/n;
###############################################################
# Independent Cox PH
###############################################################
# MCMC parameters
nburn=500; nsave=300; 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, r0=1);
# Fit the Cox PH model
res1 = indeptCoxph(formula = Surv(tobs, delta)~x1+x2, data=d,
prior=prior, mcmc=mcmc);
sfit1=summary(res1); sfit1;
## traceplot
par(mfrow = c(2,2))
traceplot(mcmc(res1$beta[1,]), main="beta1")
traceplot(mcmc(res1$beta[2,]), main="beta2")
traceplot(mcmc(res1$h.scaled[2,]), main="h")
traceplot(mcmc(res1$h.scaled[3,]), main="h")
############################################
## Curves
############################################
wide=0.1;
tgrid = seq(1e-10,4,wide);
ngrid = length(tgrid);
p = length(betaT); # number of covariates
xpred = rbind(c(0,0), c(0,1));
estimates=plot(res1, xpred=xpred, tgrid=tgrid);
Shat = estimates$Shat;
## plot
par(mfrow = c(1,1))
plot(tgrid, Sioft(tgrid, xpred[2,]), "l", lwd=3);
lines(tgrid, Sioft(tgrid, xpred[1,]), "l", lwd=3);
lines(estimates$tgrid, estimates$Shat[,1], lty=2, lwd=3)
lines(estimates$tgrid, estimates$Shatlow[,1], lty=3, lwd=1)
lines(estimates$tgrid, estimates$Shatup[,1], lty=3, lwd=1)
lines(estimates$tgrid, estimates$Shat[,2], lty=2, lwd=3)
lines(estimates$tgrid, estimates$Shatlow[,2], lty=3, lwd=1)
lines(estimates$tgrid, estimates$Shatup[,2], lty=3, lwd=1)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.