| maple.ph {mable} | R Documentation |
Mable fit of the PH model with given regression coefficients
Description
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.
Usage
maple.ph(
formula,
data,
M,
g,
pi0 = NULL,
p = NULL,
tau = Inf,
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 |
the given |
pi0 |
Initial guess of |
p |
an initial coefficients of Bernstein polynomial model of degree |
tau |
right endpoint of support |
x0 |
a working baseline covariate. See 'Details'. |
controls |
Object of class |
progress |
if |
Details
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.
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.
Value
a class 'mable_reg' object, a list with components
-
Mthe vector(m0, m1), wherem1is the last candidate degree when the search stoped -
mthe selected optimal degreem -
pthe estimate ofp = (p_0, \dots, p_m,p_{m+1}), the coefficients of Bernstein polynomial of degreem -
coefficientsthe given regression coefficients of the PH model -
mloglikthe maximum log-likelihood at an optimal degreem -
lklog-likelihoods evaluated atm \in \{m_0, \ldots, m_1\} -
lrlikelihood ratios for change-points evaluated atm \in \{m_0+1, \ldots, m_1\} -
tau.nmaximum observed time\tau_n -
tauright endpoint of support[0, \tau) -
x0the working baseline covariates -
egx0the value ofe^{\gamma'x_0} -
convergencean integer code. 0 indicates successful completion(the iteration is convergent). 1 indicates that the maximum candidate degree had been reached in the calculation; -
deltathe final convergence criterion for EM iteration; -
chptsthe change-points among the candidate degrees; -
pomthe p-value of the selected optimal degreemas a change-point;
Author(s)
Zhong Guan <zguan@iusb.edu>
References
Guan, Z. (2019) Maximum Approximate Bernstein Likelihood Estimation in Proportional Hazard Model for Interval-Censored Data, arXiv:1906.08882 .
See Also
Examples
## 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(x=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(x=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)