mable.aft {mable} | R Documentation |
Mable fit of Accelerated Failure Time Model
Description
Maximum approximate Bernstein/Beta likelihood estimation for accelerated failure time model based on interval censored data.
Usage
mable.aft(
formula,
data,
M,
g = NULL,
p = NULL,
tau = NULL,
x0 = NULL,
controls = mable.ctrl(),
progress = TRUE
)
Arguments
formula |
regression formula. Response must be |
data |
a dataset |
M |
a positive integer or a vector |
g |
a |
p |
an initial coefficients of Bernstein polynomial of degree |
tau |
the right endpoint of the support or truncation interval |
x0 |
a working baseline covariate |
controls |
Object of class |
progress |
if |
Details
Consider the accelerated failure time model with covariate for interval-censored failure time data:
S(t|x) = S(t \exp(\gamma^T(x-x_0))|x_0)
, where x_0
is a baseline covariate.
Let f(t|x)
and F(t|x) = 1-S(t|x)
be the density and cumulative distribution
functions of the event time given X = x
, respectively.
Then f(t|x_0)
on a truncation interval [0, \tau]
can be approximated by
f_m(t|x_0; p) = \tau^{-1}\sum_{i=0}^m p_i\beta_{mi}(t/\tau)
,
where p_i\ge 0
, i = 0, \ldots, m
, \sum_{i=0}^mp_i=1
,
\beta_{mi}(u)
is the beta denity with shapes i+1
and m-i+1
, and
\tau
is larger than the largest observed time, either uncensored time, or right endpoint of interval/left censored,
or left endpoint of right censored time. So we can approximate S(t|x_0)
on [0, \tau]
by
S_m(t|x_0; p) = \sum_{i=0}^{m} p_i \bar B_{mi}(t/\tau)
, where \bar B_{mi}(u)
is
the beta survival function with shapes i+1
and m-i+1
.
Response variable should be of the form cbind(l, u)
, where (l,u)
is the interval
containing the event time. Data is uncensored if l = u
, right censored
if u = Inf
or u = NA
, and left censored data if l = 0
.
The truncation time tau
and the baseline x0
should be chosen so that
S(t|x)=S(t \exp(\gamma^T(x-x_0))|x_0)
on [\tau, \infty)
is negligible for
all the observed x
.
The search for optimal degree m
stops if either m1
is reached or the test
for change-point results in a p-value pval
smaller than sig.level
.
Value
A list with components
-
m
the given or selected optimal degreem
-
p
the estimate ofp = (p_0, ..., p_m)
, the coefficients of Bernstein polynomial of degreem
-
coefficients
the estimated regression coefficients of the AFT model -
SE
the standard errors of the estimated regression coefficients -
z
the z-scores of the estimated regression coefficients -
mloglik
the maximum log-likelihood at an optimal degreem
-
tau.n
maximum observed time\tau_n
-
tau
right endpoint of trucation interval[0, \tau)
-
x0
the working baseline covariates -
egx0
the value ofe^{\gamma^T x_0}
-
convergence
an integer code: 0 indicates a successful completion; 1 indicates that the search of an optimal degree using change-point method reached the maximum candidate degree; 2 indicates that the matimum iterations was reached for calculating\hat p
and\hat\gamma
with the selected degreem
, or the divergence of the last EM-like iteration forp
or the divergence of the last (quasi) Newton iteration for\gamma
; 3 indicates 1 and 2. -
delta
the finaldelta
ifm0 = m1
or the finalpval
of the change-point for searching the optimal degreem
;
and, if m0<m1
,
-
M
the vector(m0, m1)
, wherem1
is the last candidate when the search stoped -
lk
log-likelihoods evaluated atm \in \{m_0, \ldots, m_1\}
-
lr
likelihood ratios for change-points evaluated atm \in \{m_0+1, \ldots, m_1\}
-
pval
the p-values of the change-point tests for choosing optimal model degree -
chpts
the change-points chosen with the given candidate model degrees
Author(s)
Zhong Guan <zguan@iusb.edu>
References
Guan, Z. (2019) Maximum Approximate Likelihood Estimation in Accelerated Failure Time Model for Interval-Censored Data, arXiv:1911.07087.
See Also
Examples
## Breast Cosmesis Data
bcos=cosmesis
bcos2<-data.frame(bcos[,1:2], x=1*(bcos$treat=="RCT"))
g <- 0.41 #Hanson and Johnson 2004, JCGS
aft.res<-mable.aft(cbind(left, right)~x, data=bcos2, M=c(1, 30), g=g, tau=100, x0=1)
op<-par(mfrow=c(1,2), lwd=1.5)
plot(x=aft.res, which="likelihood")
plot(x=aft.res, y=data.frame(x=0), which="survival", model='aft', type="l", col=1,
add=FALSE, main="Survival Function")
plot(x=aft.res, y=data.frame(x=1), which="survival", model='aft', lty=2, col=1)
legend("bottomleft", bty="n", lty=1:2, col=1, c("Radiation Only", "Radiation and Chemotherapy"))
par(op)