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 cbind. See 'Details'.

data

a dataset

M

a positive integer or a vector (m0, m1). If M = m0 or m0 = m1, then m0 is a preselected degree. If m0 < m1 it specifies the set of consective candidate model degrees m0:m1 for searching an optimal degree, where m1-m0>3.

g

the given dd-vector of regression coefficients

pi0

Initial guess of π(x0)=F(τnx0)\pi(x_0) = F(\tau_n|x_0). Without right censored data, pi0 = 1. See 'Details'.

p

an initial coefficients of Bernstein polynomial model of degree m0, default is the uniform initial.

tau

right endpoint of support [0,τ)[0, \tau) must be greater than or equal to the maximum observed time

x0

a working baseline covariate. See 'Details'.

controls

Object of class mable.ctrl() specifying iteration limit and other control options. Default is mable.ctrl.

progress

if TRUE a text progressbar is displayed

Details

Consider Cox's PH model with covariate for interval-censored failure time data: S(tx)=S(tx0)exp(γT(xx0))S(t|x) = S(t|x_0)^{\exp(\gamma^T(x-x_0))}, where x0x_0 satisfies γT(xx0)0\gamma^T(x-x_0)\ge 0. Let f(tx)f(t|x) and F(tx)=1S(tx)F(t|x) = 1-S(t|x) be the density and cumulative distribution functions of the event time given X=xX = x, respectively. Then f(tx0)f(t|x_0) on [0,τn][0,\tau_n] can be approximated by fm(tx0;p)=τn1i=0mpiβmi(t/τn)f_m(t|x_0; p) = \tau_n^{-1}\sum_{i=0}^m p_i\beta_{mi}(t/\tau_n), where pi0p_i \ge 0, i=0,,mi = 0, \ldots, m, i=0mpi=1pm+1\sum_{i=0}^mp_i = 1-p_{m+1}, βmi(u)\beta_{mi}(u) is the beta denity with shapes i+1i+1 and mi+1m-i+1, and τn\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(tx0)S(t|x_0) on [0,τn][0, \tau_n] by Sm(tx0;p)=i=0m+1piBˉmi(t/τn)S_m(t|x_0; p) = \sum_{i=0}^{m+1} p_i \bar B_{mi}(t/\tau_n), where Bˉmi(u)\bar B_{mi}(u), i=0,,mi = 0, \ldots, m, is the beta survival function with shapes i+1i+1 and mi+1m-i+1, Bˉm,m+1(t)=1\bar B_{m,m+1}(t) = 1, pm+1=1π(x0)p_{m+1} = 1-\pi(x_0), and π(x0)=F(τnx0)\pi(x_0) = F(\tau_n|x_0). For data without right-censored time, pm+1=1π(x0)=0.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 dd columns. The baseline x0 should chosen so that γT(xx0)\gamma^T(x-x_0) is nonnegative for all the observed xx.

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

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

mable.ph

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)


[Package mable version 3.1.3 Index]