crisk2.bart | R Documentation |
Here we have implemented another approach to utilize BART for competing risks that is very flexible, and is akin to discrete-time survival analysis. Following the capabilities of BART, we allow for maximum flexibility in modeling the dependence of competing failure times on covariates. In particular, we do not impose proportional hazards.
Similar to crisk.bart
, we utilize two BART models, yet they are
two different BART models than previously considered. First, given an
event of either cause occurred, we employ a typical binary BART model to
discriminate between cause 1 and 2. Next, we proceed as if it were a typical
survival analysis with BART for an absorbing event from either cause.
To elaborate, consider data in the form: (s_i, \delta_i,
{x}_i)
where s_i
is the event time;
\delta_i
is an indicator distinguishing events,
\delta_i=h
due to cause h in {1, 2}
, from
right-censoring, \delta_i=0
; {x}_i
is a vector of
covariates; and i=1, ..., N
indexes subjects.
We denote the K
distinct event/censoring times by
0<t_{(1)}<...<t_{(K)}<\infty
thus
taking t_{(j)}
to be the j^{th}
order statistic
among distinct observation times and, for convenience,
t_{(0)}=0
.
First, consider event indicators for an event from either cause:
y_{1ij}
for each subject i
at each distinct time
t_{(j)}
up to and including the subject's last observation
time s_i=t_{(n_i)}
with n_i=\arg \max_j [t_{(j)}\leq
s_i]
. We denote by p_{1ij}
the
probability of an event at time t_{(j)}
conditional on no
previous event. We now write the model for y_{1ij}
as a
nonparametric probit (or logistic) regression of y_{1ij}
on
the time t_{(j)}
and the covariates {x}_{1i}
,
and then utilize BART for binary responses. Specifically,
y_{1ij}\ =\ I[\delta_i>0] I[s_i=t_{(j)}],\ j=1, ..., n_i
. Therefore, we have p_{1ij} =
F(mu_{1ij}),\ mu_{1ij} = mu_1+f_1(t_{(j)}, {x}_{1i})
where F
denotes the Normal (or
Logistic) cdf.
Next, we denote by p_{2i}
the probability of a cause 1
event at time s_i
conditional on an event having
occurred. We now write the model for y_{2i}
as a
nonparametric probit (or logistic) regression of y_{2i}
on
the time s_i
and the covariates {x}_{2i}
,
via BART for binary responses. Specifically,
y_{2i}\ =\ I[\delta_i=1]
. Therefore, we
have p_{2i} = F(mu_{2i}),\ mu_{2i} = mu_2+f_2(s_i,
{x}_{2i})
where
F
denotes the Normal (or Logistic) cdf. Although, we modeled
p_{2i}
at the time of an event, s_i
, we can
estimate this probability at any other time points on the grid via
p(t_{(j)}, x_2)=F( mu_2+f_2(t_{(j)}, {x}_2))
.
Finally, based on these probabilities,
p_{hij}
, we can construct targets of inference such as the
cumulative incidence functions.
crisk2.bart(x.train=matrix(0,0,0), y.train=NULL,
x.train2=x.train, y.train2=NULL,
times=NULL, delta=NULL, K=NULL,
x.test=matrix(0,0,0), x.test2=x.test,
sparse=FALSE, theta=0, omega=1,
a=0.5, b=1, augment=FALSE,
rho=NULL, rho2=NULL,
xinfo=matrix(0,0,0), xinfo2=matrix(0,0,0),
usequants=FALSE,
rm.const=TRUE, type='pbart',
ntype=as.integer(
factor(type, levels=c('wbart', 'pbart', 'lbart'))),
k=2, power=2, base=0.95,
offset=NULL, offset2=NULL,
tau.num=c(NA, 3, 6)[ntype],
ntree=50, numcut=100, ndpost=1000, nskip=250,
keepevery = 10L,
printevery=100L,
id=NULL, ## crisk2.bart only
seed=99, ## mc.crisk2.bart only
mc.cores=2, ## mc.crisk2.bart only
nice=19L ## mc.crisk2.bart only
)
mc.crisk2.bart(x.train=matrix(0,0,0), y.train=NULL,
x.train2=x.train, y.train2=NULL,
times=NULL, delta=NULL, K=NULL,
x.test=matrix(0,0,0), x.test2=x.test,
sparse=FALSE, theta=0, omega=1,
a=0.5, b=1, augment=FALSE,
rho=NULL, rho2=NULL,
xinfo=matrix(0,0,0), xinfo2=matrix(0,0,0),
usequants=FALSE,
rm.const=TRUE, type='pbart',
ntype=as.integer(
factor(type, levels=c('wbart', 'pbart', 'lbart'))),
k=2, power=2, base=0.95,
offset=NULL, offset2=NULL,
tau.num=c(NA, 3, 6)[ntype],
ntree=50, numcut=100, ndpost=1000, nskip=250,
keepevery = 10L,
printevery=100L,
id=NULL, ## crisk2.bart only
seed=99, ## mc.crisk2.bart only
mc.cores=2, ## mc.crisk2.bart only
nice=19L ## mc.crisk2.bart only
)
x.train |
Covariates for training (in sample) data for an event. |
y.train |
Event binary response for training (in sample)
data. |
x.train2 |
Covariates for training (in sample)
data of for a cause 1 event. Similar to |
y.train2 |
Cause 1 event binary response for training (in sample) data.
Similar to |
times |
The time of event or right-censoring, |
delta |
The event indicator: 1 for cause 1, 2 for cause 2 and 0 is censored. |
K |
If provided, then coarsen |
x.test |
Covariates for test (out of sample) data of an event. |
x.test2 |
Covariates for test (out of sample) data of a cause 1 event.
Similar to |
sparse |
Whether to perform variable selection based on a sparse Dirichlet prior; see Linero 2016. |
theta |
Set |
omega |
Set |
a |
Sparse parameter for |
b |
Sparse parameter for |
rho |
Sparse parameter: typically |
rho2 |
Sparse parameter: typically |
augment |
Whether data augmentation is to be performed in sparse variable selection. |
xinfo |
You can provide the cutpoints to BART or let BART
choose them for you. To provide them, use the |
xinfo2 |
Cause 2 cutpoints. |
usequants |
If |
rm.const |
Whether or not to remove constant variables. |
type |
Whether to employ probit BART via Albert-Chib,
|
ntype |
The integer equivalent of |
k |
k is the number of prior standard deviations
|
power |
Power parameter for tree prior. |
base |
Base parameter for tree prior. |
offset |
Cause 1 binary offset. |
offset2 |
Cause 2 binary offset. |
tau.num |
The numerator in the |
ntree |
The number of trees in the sum. |
numcut |
The number of possible values of cutpoints (see
|
ndpost |
The number of posterior draws returned. |
nskip |
Number of MCMC iterations to be treated as burn in. |
keepevery |
Every |
printevery |
As the MCMC runs, a message is printed every |
id |
|
seed |
|
mc.cores |
|
nice |
|
crisk2.bart
returns an object of type crisk2bart
which is
essentially a list. Besides the items listed
below, the list has offset
, offset2
,
times
which are the unique times, K
which is the number of unique times, tx.train
and
tx.test
, if any.
yhat.train |
A matrix with |
yhat.test |
Same as |
surv.test |
test data fits for the survival function, |
surv.test.mean |
mean of |
prob.test |
The probability of suffering an event. |
prob.test2 |
The probability of suffering a cause 1 event. |
cif.test |
The cumulative incidence function of cause 1,
|
cif.test2 |
The cumulative incidence function of cause 2,
|
cif.test.mean |
mean of |
cif.test2.mean |
mean of |
varcount |
a matrix with |
varcount2 |
For each variable the total count of the number of times this variable is used for a cause 1 event in a tree decision rule is given. |
surv.pre.bart
, predict.crisk2bart
,
mc.crisk2.pwbart
, crisk.bart
data(transplant)
pfit <- survfit(Surv(futime, event) ~ abo, transplant)
# competing risks for type O
plot(pfit[4,], xscale=7, xmax=735, col=1:3, lwd=2, ylim=c(0, 1),
xlab='t (weeks)', ylab='Aalen-Johansen (AJ) CI(t)')
legend(450, .4, c("Death", "Transplant", "Withdrawal"), col=1:3, lwd=2)
## plot(pfit[4,], xscale=30.5, xmax=735, col=1:3, lwd=2, ylim=c(0, 1),
## xlab='t (months)', ylab='Aalen-Johansen (AJ) CI(t)')
## legend(450, .4, c("Death", "Transplant", "Withdrawal"), col=1:3, lwd=2)
delta <- (as.numeric(transplant$event)-1)
## recode so that delta=1 is cause of interest; delta=2 otherwise
delta[delta==1] <- 4
delta[delta==2] <- 1
delta[delta>1] <- 2
table(delta, transplant$event)
times <- pmax(1, ceiling(transplant$futime/7)) ## weeks
##times <- pmax(1, ceiling(transplant$futime/30.5)) ## months
table(times)
typeO <- 1*(transplant$abo=='O')
typeA <- 1*(transplant$abo=='A')
typeB <- 1*(transplant$abo=='B')
typeAB <- 1*(transplant$abo=='AB')
table(typeA, typeO)
x.train <- cbind(typeO, typeA, typeB, typeAB)
x.test <- cbind(1, 0, 0, 0)
dimnames(x.test)[[2]] <- dimnames(x.train)[[2]]
##test BART with token run to ensure installation works
set.seed(99)
post <- crisk2.bart(x.train=x.train, times=times, delta=delta,
x.test=x.test, nskip=1, ndpost=1, keepevery=1)
## Not run:
## run one long MCMC chain in one process
## set.seed(99)
## post <- crisk2.bart(x.train=x.train, times=times, delta=delta, x.test=x.test)
## in the interest of time, consider speeding it up by parallel processing
## run "mc.cores" number of shorter MCMC chains in parallel processes
post <- mc.crisk2.bart(x.train=x.train, times=times, delta=delta,
x.test=x.test, seed=99, mc.cores=8)
K <- post$K
typeO.cif.mean <- apply(post$cif.test, 2, mean)
typeO.cif.025 <- apply(post$cif.test, 2, quantile, probs=0.025)
typeO.cif.975 <- apply(post$cif.test, 2, quantile, probs=0.975)
plot(pfit[4,], xscale=7, xmax=735, col=1:3, lwd=2, ylim=c(0, 0.8),
xlab='t (weeks)', ylab='CI(t)')
points(c(0, post$times)*7, c(0, typeO.cif.mean), col=4, type='s', lwd=2)
points(c(0, post$times)*7, c(0, typeO.cif.025), col=4, type='s', lwd=2, lty=2)
points(c(0, post$times)*7, c(0, typeO.cif.975), col=4, type='s', lwd=2, lty=2)
legend(450, .4, c("Transplant(BART)", "Transplant(AJ)",
"Death(AJ)", "Withdrawal(AJ)"),
col=c(4, 2, 1, 3), lwd=2)
##dev.copy2pdf(file='../vignettes/figures/liver-BART.pdf')
## plot(pfit[4,], xscale=30.5, xmax=735, col=1:3, lwd=2, ylim=c(0, 0.8),
## xlab='t (months)', ylab='CI(t)')
## points(c(0, post$times)*30.5, c(0, typeO.cif.mean), col=4, type='s', lwd=2)
## points(c(0, post$times)*30.5, c(0, typeO.cif.025), col=4, type='s', lwd=2, lty=2)
## points(c(0, post$times)*30.5, c(0, typeO.cif.975), col=4, type='s', lwd=2, lty=2)
## legend(450, .4, c("Transplant(BART)", "Transplant(AJ)",
## "Death(AJ)", "Withdrawal(AJ)"),
## col=c(4, 2, 1, 3), lwd=2)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.