View source: R/mable-ph-model.r
maple.ph | R Documentation |
Maximum approximate profile likelihood estimation of Bernstein polynomial model in Cox's proportional hazards regression based on interal censored event time data with given regression coefficients which are efficient estimates provided by other semiparametric methods.
maple.ph(
formula,
data,
M,
g,
pi0 = NULL,
p = NULL,
tau = Inf,
x0 = NULL,
controls = mable.ctrl(),
progress = TRUE
)
formula |
regression formula. Response must be |
data |
a data frame containing variables in |
M |
a positive integer or a vector |
g |
the given |
pi0 |
Initial guess of |
p |
an initial coefficients of Bernstein polynomial model of degree |
tau |
right endpoint of support |
x0 |
a data frame specifying working baseline covariates on the right-hand-side of |
controls |
Object of class |
progress |
if |
Consider Cox's PH model with covariate for interval-censored failure time data:
S(t|x) = S(t|x_0)^{\exp(\gamma^T(x-x_0))}
, where x_0
satisfies \gamma^T(x-x_0)\ge 0
,
where x
and x_0
may
contain dummy variables and interaction terms. The working baseline x0
in arguments
contains only the values of terms excluding dummy variables and interaction terms
in the right-hand-side of formula
. Thus g
is the initial guess of
the coefficients \gamma
of x-x_0
and could be longer than x0
.
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 [0,\tau_n]
can be approximated by
f_m(t|x_0; p) = \tau_n^{-1}\sum_{i=0}^m p_i\beta_{mi}(t/\tau_n)
,
where p_i \ge 0
, i = 0, \ldots, m
, \sum_{i=0}^mp_i = 1-p_{m+1}
,
\beta_{mi}(u)
is the beta denity with shapes i+1
and m-i+1
, and
\tau_n
is 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_n]
by
S_m(t|x_0; p) = \sum_{i=0}^{m+1} p_i \bar B_{mi}(t/\tau_n)
, where
\bar B_{mi}(u)
, i = 0, \ldots, m
, is the beta survival function with shapes
i+1
and m-i+1
, \bar B_{m,m+1}(t) = 1
, p_{m+1} = 1-\pi(x_0)
, and
\pi(x_0) = F(\tau_n|x_0)
. For data without right-censored time, p_{m+1} = 1-\pi(x_0) = 0.
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 associated covariate contains d
columns. The baseline x0
should chosen so that
\gamma^T(x-x_0)
is nonnegative 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
.
a class 'mable_reg
' object, a list with components
M
the vector (m0, m1)
, where m1
is the last candidate degree when the search stoped
m
the selected optimal degree m
p
the estimate of p = (p_0, \dots, p_m,p_{m+1})
, the coefficients of Bernstein polynomial of degree m
coefficients
the given regression coefficients of the PH model
mloglik
the maximum log-likelihood at an optimal degree m
lk
log-likelihoods evaluated at m \in \{m_0, \ldots, m_1\}
lr
likelihood ratios for change-points evaluated at m \in \{m_0+1, \ldots, m_1\}
tau.n
maximum observed time \tau_n
tau
right endpoint of support [0, \tau)
x0
the working baseline covariates
egx0
the value of e^{\gamma'x_0}
convergence
an integer code. 0 indicates successful completion(the iteration is
convergent). 1 indicates that the maximum candidate degree had been reached in the calculation;
delta
the final convergence criterion for EM iteration;
chpts
the change-points among the candidate degrees;
pom
the p-value of the selected optimal degree m
as a change-point;
Zhong Guan <zguan@iu.edu>
Guan, Z. (2019) Maximum Approximate Bernstein Likelihood Estimation in Proportional Hazard Model for Interval-Censored Data, arXiv:1906.08882 .
mable.ph
## Simulated Weibull data
require(icenReg)
set.seed(123)
simdata<-simIC_weib(70, inspections = 5, inspectLength = 1)
sp<-ic_sp(cbind(l, u) ~ x1 + x2, data = simdata)
res0<-maple.ph(cbind(l, u) ~ x1 + x2, data = simdata, M=c(2,20),
g=sp$coefficients, tau=7)
op<-par(mfrow=c(1,2))
plot(res0, which=c("likelihood","change-point"))
par(op)
res1<-mable.ph(cbind(l, u) ~ x1 + x2, data = simdata, M=res0$m,
g=c(.5,-.5), tau=7)
op<-par(mfrow=c(1,2))
plot(res1, y=data.frame(x1=0, x2=0), which="density", add=FALSE, type="l",
xlab="Time", main="Desnity Function")
lines(xx<-seq(0, 7, len=512), dweibull(xx, 2,2), lty=2, col=2)
legend("topright", bty="n", lty=1:2, col=1:2, c("Estimated","True"))
plot(res1, y=data.frame(x1=0, x2=0), which="survival", add=FALSE, type="l",
xlab="Time", main="Survival Function")
lines(xx, 1-pweibull(xx, 2, 2), lty=2, col=2)
legend("topright", bty="n", lty=1:2, col=1:2, c("Estimated","True"))
par(op)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.