maple.po {mable} | R Documentation |
Mable fit of the PO model with given regression coefficients
Description
Maximum approximate profile likelihood estimation of Bernstein polynomial model in proportional odds rate regression based on interal censored event time data with given regression coefficients and select an optimal degree m if coefficients are efficient estimates provided by other semiparametric methods.
Usage
maple.po(
formula,
data,
M,
g,
tau,
x0 = NULL,
controls = mable.ctrl(),
progress = TRUE
)
Arguments
formula |
regression formula. Response must be |
data |
a data frame containing variables in |
M |
a positive integer or a vector |
g |
the given |
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 |
Details
Consider Generalized PO model with covariate for interval-censored failure time data:
S(t|x) = S(t|x_0)^{\exp(\gamma'(x-x_0))}
, where x_0
satisfies
\gamma'(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^{-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 the right endpoint of support interval. 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)
, i = 0, \ldots, m
, 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 are
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'(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
-
M
the vector(m0, m1)
, wherem1
is the last candidate degree when the search stoped -
m
the selected optimal degreem
-
p
the estimate ofp = (p_0, \dots, p_m,p_{m+1})
, the coefficients of Bernstein polynomial of degreem
-
coefficients
the given regression coefficients of the PH model -
mloglik
the maximum log-likelihood at an optimal degreem
-
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\}
-
tau.n
maximum observed time\tau_n
-
tau
right endpoint of support[0, \tau)
-
x0
the working baseline covariates -
egx0
the value ofe^{\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 degreem
as a change-point;
Author(s)
Zhong Guan <zguan@iu.edu>
References
Guan, Z. et al. (???) Maximum Approximate Bernstein Likelihood Estimation in Generalized Proportional Odds Regression Model for Interval-Censored Data
See Also
Examples
## Simulated Weibull data
require(icenReg)
set.seed(111)
simdata<-simIC_weib(100, model = "po", inspections = 2,
inspectLength = 2.5, prob_cen=1)
sp<-ic_sp(cbind(l, u) ~ x1 + x2, data = simdata, model="po")
gt<--sp$coefficients
res0<-maple.po(cbind(l, u) ~ x1 + x2, data = simdata, M=c(1,20), g=gt, tau=6)
op<-par(mfrow=c(1,2))
plot(res0, which=c("likelihood","change-point"))
par(op)
res1<-mable.po(cbind(l, u) ~ x1 + x2, data = simdata, M=c(1,20),
g=gt, tau=6, x0=data.frame(x1=max(simdata$x1),x2=-1))
op<-par(mfrow=c(2,2))
plot(res1, which=c("likelihood","change-point"))
plot(res0, y=data.frame(x1=0,x2=0), which="density", add=FALSE, type="l",
xlab="Time", main="Desnity Function")
plot(res1, y=data.frame(x1=0,x2=0), which="density", add=TRUE, lty=2, col=4)
lines(xx<-seq(0, 7, len=512), dweibull(xx, 2,2), lty=3, col=2, lwd=1.5)
legend("topright", bty="n", lty=1:3, col=c(1,4,2), c(expression(hat(f)[0]),
expression(tilde(f)[0]), expression(f[0])))
plot(res0, y=data.frame(x1=0,x2=0), which="survival", add=FALSE, type="l",
xlab="Time", main="Survival Function")
plot(res1, y=data.frame(x1=0,x2=0), which="survival", add=TRUE, lty=2, col=4)
lines(xx, 1-pweibull(xx, 2, 2), lty=2, col=2)
legend("topright", bty="n", lty=1:3, col=c(1,4,2), c(expression(hat(S)[0]),
expression(tilde(S)[0]), expression(S[0])))
par(op)