mable.ph {mable}

Mable fit of Cox's proportional hazards regression model

Description

Maximum approximate Bernstein/Beta likelihood estimation in Cox's proportional hazards regression model based on interal censored event time data.

Usage

mable.ph(
  formula,
  data,
  M,
  g = NULL,
  p = NULL,
  pi0 = 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 = m 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

initial guess of d-vector of regression coefficients. See 'Details'.

p

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

pi0

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

tau

right endpoint of support [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(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'(x-x_0) is nonnegative for all the observed x and all \gamma in a neighborhood of its true value.

A missing initial value of g is imputed by ic_sp() of package icenReg.

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. This process takes longer than maple.ph to select an optimal degree.

Value

A list with components

• m the selected/preselected 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 estimated regression coefficients of the PH 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 degree m

• 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, 1 indicates either the EM-like iteration for finding maximum likelihood reached the maximum iteration for at least one m or the search of an optimal degree using change-point method reached the maximum candidate degree, 2 indicates both occured, and 0 indicates a successful completion.

• delta the final delta if m0 = m1 or the final pval of the change-point for searching the optimal degree m;

and, if m0<m1,

• M the vector (m0, m1), where m1 is the last candidate degree when the search stoped

• 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\}

• 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 <[e-mail]>

References

Guan, Z. Maximum Approximate Bernstein Likelihood Estimation in Proportional Hazard Model for Interval-Censored Data, Statistics in Medicine. 2020; 1–21. https://doi.org/10.1002/sim.8801.

See Also

maple.ph

Examples

   # Ovarian Cancer Survival Data
   require(survival)
   futime2<-ovarian$futime
   futime2[ovarian$fustat==0]<-Inf
   ovarian2<-data.frame(age=ovarian$age, futime1=ovarian$futime, 
        futime2=futime2)
   ova<-mable.ph(cbind(futime1, futime2) ~ age, data = ovarian2, 
        M=c(2,35), g=.16, x0=35)
   op<-par(mfrow=c(2,2))
   plot(ova, which = "likelihood")
   plot(ova, which = "change-point")
   plot(ova, y=data.frame(age=60), which="survival", add=FALSE, type="l", 
         xlab="Days", main="Age = 60")
   plot(ova, y=data.frame(age=65), which="survival", add=FALSE, type="l", 
         xlab="Days", main="Age = 65")
   par(op)