Maximum likelihood estimation in proportional odds regression model based on interval-censored event-time data

Abstract

Maximum likelihood estimates of density function and regression coefficients in the proportional odds regression models are proposed and studied based on event-time data that are either completely or partly interval-censored. A smooth estimate of the survival function is then obtained. Theoretical results indicate that the proposed method enjoys an almost parametric $\sqrt{n}$ -consistency. Some simulation studies show that the proposed method not only gives density and smooth survival curve estimates but also outperforms the existing semiparametric method in terms of estimating both the regression coefficients and the survival curves for small and medium sample sizes. The proposed method is illustrated by an application to the HIV infection data.

Keywords

Bernstein polynomial model density estimation interval censoring proportional odds survival curve

1. Introduction

A limited number of parametric families of distributions are available for survival analysis. These parametric models are restrictive due to the small number of parameters. Nonparametric and semiparametric models such as the accelerate-failure-time (AFT) model,^1–4 the proportional hazard (PH) regression model,⁵ and the proportional odds (PO) model^6–9 are much more flexible and thus are preferred.

Let $T$ be an event-time and $X$ be an associated $d$ -dimensional covariate with distribution $H (x)$ on $X$ . We denote the marginal and the conditional survival probability of $T$ , respectively, by $S (t) = 1 - F (t) = P (T > t)$ and $S (t ∣ x) = 1 - F (t | x) = P (T > t | X = x) .$ Define the odds function at time $t > 0$ as $O (t | x) = F (t | x) / S (t | x)$ . As an important alternative to the PH model, the PO model can be defined by

O (t | x) = e^{γ^{⊤} \tilde{x}} O (t | x_{0}),

(1)

where

\tilde{x} = x - x_{0}

and

O (t | x_{0})

is an unknown increasing baseline odds function on

(0, \infty)

and

γ \in G \subset R^{d}

are unknown regression coefficients. This model is a special case of the general linear transformation model:^10,11

μ (T) = - γ^{⊤} X + ε

with

μ (t) = \log O (t | x_{0})

and

ε \sim G (t) = 1 - \exp (- e^{t})

. A very computer-intensive likelihood sampler method and a local approximate method were proposed¹¹ to obtain maximum partial likelihood estimator of

γ

in the general linear transformation model. A semiparametric efficient estimator of

γ

in a generalized odds rate regression model for right-censored data was proposed by Scharfstein et al.¹²

When data are interval-censored,^13–15 the time $T$ is only known to be in an interval $Y = (Y_{1}, Y_{2}]$ . For Case 1 data, that is, the current status or doubly censored data with one examination time $U$ , $T$ is left-censored with $Y = (0, U]$ or right-censored with $Y = (U, \infty)$ . We reduce the cases with more than two examination times to the case with two examination times,¹⁶ that is, the Case $2$ data. In the latter case, in addition to left-censoring with $Y = (0, U_{1}]$ and right-censoring with $Y = (U_{2}, \infty)$ , $T$ could also be interval-censored with $Y = (Y_{1}, Y_{2}] = (U_{1}, U_{2}]$ . For Case 1 data the distribution function of the examination time $U$ given covariates $X = x$ is denoted by $G (u | x)$ . For Case 2 data, the joint distribution function of the observed examination times $U = (U_{1}, U_{2})$ given $X = x$ is denoted by $G (u | x)$ . In order to deal with data which contain uncensored $T$ , we assume that the observable random variables are $Z = (X, Y, Δ)$ , where $Y = (Y_{1}, Y_{2}]$ and $Δ = I (Y_{1} < Y_{2})$ is the censoring indicator, that is, uncensored $T = Y = Y_{1} = Y_{2}$ if $Δ = 0$ , and interval-censored $T \in Y = (Y_{1}, Y_{2}]$ , $0 \leq Y_{1} < Y_{2} \leq \infty$ , if $Δ = 1$ . For convenience we call the uncensored data the Case 0 data.

Traditional treatment of the nonparametric portion of a non- or a semi-parametric model is to approximate the unknown underlying continuous distribution function by a step-function with jumps at the observations and to parameterize it by the unknown jump sizes.^17–19 When data are not censored this method results in efficient estimation. However, especially for small sample size data, the estimated distribution and density functions are unsatisfactory due to its roughness. They are even more unsatisfactory for interval-censored data. Turnbull¹⁸ developed a self-consistency algorithm which can be viewed as an example of EM algorithm for estimating distribution function based on arbitrarily grouped, censored and truncated data. This approach does not work for the AFT regression model with interval-censoring. Although the estimation of the regression coefficients is efficient,¹⁶ the distributional estimation based on interval censored data is typically $n^{1 / 3}$ -consistent and not unique on the so-called Turnbull intervals. In addition to the works reviewed by Huang and Wellner¹⁶ and the references therein, several authors contributed with more advanced methods especially the sieve methods with or without smooth approximation of the survival function for the semiparametric PO model. Huang and Rossini (1997)²⁰ presented a sieve estimation for the PO model with interval censored times. Murphy et al. (1997)²¹ considered maximum likelihood estimation of the regression coefficients and the baseline odds in the PO model with right-censored data. Shen (1998)²² proposed a sieve maximum likelihood estimation for the PO model using spline approximation based on right-censored and Case 2 interval censored data. Rabinowitz et al. (2000)²³ gave an easily implemented method for fitting PO model using conditional logistic regression based on interval-censored data with special settings and current status data. A Semiparametric Bayes’ PO model for Case 1 data with underreporting were studied by Wang and Dunson(2011).²⁴ Pordeli and Lu (2016)²⁵ proposed a B-spline-based method for semiparametric PO models for Case 1 interval-censored data. Spline methods have the advantage in modeling nonparametric functions with sudden changes.

Guan²⁶ proposed to approximate the unknown underlying continuous distribution function defined on a closed interval by a Bernstein polynomial model and use the coefficients as parameters. This method has the advantage over other smoothing techniques such as kernel and spline methods. The Bernstein polynomial model is an actual probability model which is the mixture of beta distributions. This makes it possible and easy to find maximum likelihood estimates based on all kinds of interval-censored data and many other incompletely observed and contaminated data.^27,28 This model was used to successfully obtain maximum likelihood estimation in the PH and AFT models based on interval censored data.^4,29 The distributional estimation has much better than $n^{1 / 3}$ rate of convergence. The Bernstein polynomial model can also improve the small sample performance of the estimates of the regression coefficients.²⁹ Although the PO model shares similarity with AFT and PH models, its specification is unique. The technical treatment of the Bernstein polynomial model applied on the PO model is not obvious and quite different from those on the AFT and PH models. As mentioned by Shen (1998)²² the estimation of $S_{0} (t)$ in PO model is difficult even for right-censored data. Most of recent researches focus on the estimation of the regression coefficients. In this paper, a new and different approach to estimating the regression coefficients together with baseline survival function in the semiparametric PO model as an important alternative to the celebrated PH model will be investigated.

The rest of the paper is organized as follows. Section 2 describes the proposed method and the algorithm to find maximum likelihood estimates. Some asymptotic results are given in Section 3. Simulation results and a real data application are presented in Sections 4 and 5. The proofs of the theoretical results are given in the Supplemental Material.

2. Methodology

2.1. The likelihood

Let $f (t | x)$ denote the conditional density of a continuous event time $T$ . For the PO model (1) we can change the “baseline” covariate $x_{0}$ to any $x_{0}^{*}$ including $x_{0}^{*} = 0$ with the same $γ$ . The PO model (1) is equivalent to

f (t | x) = f (t | x; γ) = \frac{e^{γ^{⊤} \tilde{x}} f_{0} (t)}{D^{2} (t | x; γ)},

(2)

where,

D (t | x; γ) = 1 + (e^{γ^{⊤} \tilde{x}} - 1) F_{0} (t) = e^{γ^{⊤} \tilde{x}} + (1 - e^{γ^{⊤} \tilde{x}}) S_{0} (t)

, with suppressed dependence on

x_{0}

f_{0} (\cdot) = f (\cdot | x_{0})

which is the unknown baseline density, and

F_{0} (t) = 1 - S_{0} (t) = F (t | x_{0})

is the cumulative distribution function of

T

. The equivalent integral versions of (1) are

S (t | x) = S (t | x; γ) = \frac{S_{0} (t)}{D (t | x; γ)}, F (t | x) = F (t | x; γ) = \frac{e^{γ^{⊤} \tilde{x}} F_{0} (t)}{D (t | x; γ)} .

(3)

Another equivalent version is

logit [F (t | x)] = α (t) + γ^{⊤} \tilde{x}

, where

α (t)

is a monotone-increasing function.³⁰

Consider an individual interval-censored data $z = {x, y = (y_{1}, y_{2}], δ}$ , where $δ = I (y_{1} < y_{2})$ . For uncensored $t$ , $δ = 0$ and $t = y = y_{1} = y_{2}$ . The loglikelihood is

\begin{aligned} ℓ (γ; z) & = (1 - δ) \log f (y | x; γ) + δ \log [S (y_{1} | x; γ) - S (y_{2} | x; γ)] \\ = (1 - δ) {γ^{⊤} \tilde{x} + \log f_{0} (y) - 2 \log [e^{γ^{⊤} \tilde{x}} + (1 - e^{γ^{⊤} \tilde{x}}) S_{0} (y)]} \\ + δ \log {\frac{S_{0} (y_{1})}{e^{γ^{⊤} \tilde{x}} + (1 - e^{γ^{⊤} \tilde{x}}) S_{0} (y_{1})} - \frac{S_{0} (y_{2})}{e^{γ^{⊤} \tilde{x}} + (1 - e^{γ^{⊤} \tilde{x}}) S_{0} (y_{2})}} . \end{aligned}

Let

z_{i} = {x_{i}, y_{i} = (y_{1 i}, y_{2 i}], δ_{i}}

i = 1, \dots, n

, be IID observations of

Z = {X, Y = (Y_{1}, Y_{2}], Δ}

Lemma 2.1

The first and second derivatives of $ℓ (γ; z)$ with respect to $γ$ are, respectively,

\begin{aligned} {\dot{ℓ}}_{γ} (γ; z) & = (1 - δ) [2 S (y | x; γ) - 1] \tilde{x} + δ \frac{{\dot{S}}_{γ} (y_{1} | x; γ) - {\dot{S}}_{γ} (y_{2} | x; γ)}{S (y_{1} | x; γ) - S (y_{2} | x; γ)}, \end{aligned}

(4)

\begin{aligned} {\ddot{ℓ}}_{γ γ} (γ; z) & = - 2 (1 - δ) S (t | x; γ) [1 - S (t | x; γ)] {\tilde{x}}^{\otimes 2} \\ - δ {S (y_{1} | x; γ) [1 - S (y_{1} | x; γ)] \\ + S (y_{2} | x; γ) [1 - S (y_{2} | x; γ)]} {\tilde{x}}^{\otimes 2}, \end{aligned}

(5)

where

v^{\otimes 2} = v v^{⊤}

for column vector

v

and

{\dot{S}}_{γ} (t | x; γ) = \frac{\partial}{\partial γ} S (t | x; γ) = - S (t | x; γ) [1 - S (t | x; γ)] \tilde{x} .

(6)

The full likelihood is $ℓ (γ) = \sum ℓ (γ; z_{i})$ . Its gradient and Hessian matrix with respect to $γ$ are denoted by $\nabla ℓ (γ)$ and $H_{n} (γ) = \partial^{2} ℓ (γ) / \partial γ \partial γ^{⊤}$ , respectively.

Theorem 2.2

Let $ρ_{k}$ ( $k = 0, 1, 2$ ) be the probability that $T$ is Case $k$ censored so that $ρ_{0} + ρ_{1} + ρ_{2} = 1$ . At $γ = γ_{0}$ we have, as $n \to \infty$ , with probability one,

\frac{1}{n} H_{n} (γ) = \frac{1}{n} \sum_{i = 1}^{n} \frac{\partial^{2} ℓ (γ; z_{i})}{\partial γ \partial η} \to - \sum_{k = 0}^{2} ρ_{k} \int_{X} J_{k} (x; γ) {\tilde{x}}^{\otimes 2} d H (x),

where

J_{0} (x; γ) = 1 / 3

J_{1} (x; γ) = \int_{0}^{\infty} S (u | x; γ) [1 - S (u | x; γ)] d G (u | x)

, and

\begin{aligned} J_{2} (x; γ) & = \iint_{u_{1} < u_{2}} S (u_{1} | x; γ) [1 - S (u_{2} | x; γ)] [1 - S (u_{1} | x; γ) + S (u_{2} | x; γ)] d G (u ∣ x) . \end{aligned}

Thus for a fixed

f_{0}

the Hessian matrix

H_{n} (γ)

in a neighborhood of

γ = γ_{0}

is asymptotically negative definite if

E [J_{k} (X; γ) {\tilde{X}}^{\otimes 2}]

is positive definite for

k

with

ρ_{k} > 0

Without specifying $f_{0}$ or $S_{0}$ using a finite dimensional model the exact full likelihood $ℓ (γ)$ cannot be maximized. Traditional method approximates $S_{0} (t)$ by step-function and treats the jumps at observations as unknown parameters. For censored or other types of incompletely observed data this parametrization is difficult and complicated. Fortunately, the Bernstein polynomial approximation makes the parametrization of $f_{0}$ much simple and easy.

We assume that $f_{0}$ has a known support or truncation interval $[0, τ]$ so that $P (T > τ | X = x_{0})$ is zero or can be neglected. Then we approximate $f_{0} (t)$ , $F_{0} (t)$ , and $S_{0} (t)$ on $[0, τ]$ , respectively, by $f_{0} (t) \approx f_{m} (t; p) = \frac{1}{τ} \sum_{j = 0}^{m} p_{j} β_{m j} (\frac{t}{τ})$ , $F_{0} (t) \approx F_{m} (t; p) = \sum_{j = 0}^{m} p_{j} B_{m j} (\frac{t}{τ})$ , and $S_{0} (t) \approx S_{m} (t; p)$ $= \sum_{j = 0}^{m} p_{j} {\bar{B}}_{m j} (\frac{t}{τ})$ , $t \in [0, τ],$ where $β_{m j} (t) = (m + 1) (\binom{m}{j}) t^{j} (1 - t)^{m - j}$ , $B_{m j} (t) = \int_{0}^{t} β_{m j} (u) d u,$ ${\bar{B}}_{m j} (t) = 1 - B_{m j} (t)$ , $j = 0, \dots, m$ , $S_{m} (τ; p) = 0$ , and

p \in S_{m} \equiv {p = (p_{0}, \dots, p_{m})^{⊤} \in R^{m + 1} : p_{j} \geq 0, \sum_{j = 0}^{m} p_{j} = 1} .

(7)

Here the dependence of

p

x_{0}

is suppressed. Note that

f_{m} (t; p)

is a mixture of beta densities with shapes

(j + 1, m - j - 1)

j = 0, \dots, m

Then $f (t ∣ x; γ)$ , $F (t ∣ x; γ)$ , and $S (t ∣ x; γ)$ can be approximated, respectively, by

\begin{aligned} f_{m} (t | x; θ) & = \frac{e^{γ^{⊤} \tilde{x}} f_{m} (t; p)}{D_{m}^{2} (t | x, θ)}, \\ F_{m} (t | x; θ) & = \frac{e^{γ^{⊤} \tilde{x}} F_{m} (t; p)}{D_{m} (t | x, θ)}, \\ S_{m} (t | x; θ) & = \frac{S_{m} (t; p)}{D_{m} (t | x, θ)}, \end{aligned}

where

θ = (γ, p)

D_{m} (t | x, θ) = 1 + (e^{γ^{⊤} \tilde{x}} - 1) F_{m} (t; p) = e^{γ^{⊤} \tilde{x}} + (1 - e^{γ^{⊤} \tilde{x}}) S_{m} (t; p) .

It is easy to see that

θ

is identifiable for any choice of baseline

x_{0}

If $τ \neq 1$ we divide all the observed times by $τ$ . Thus we assume in the rest of this article that $τ = 1$ . Let $β_{m} (t) = {β_{m 0} (t), \dots, β_{m m} (t)}^{⊤}$ , $B_{m} (t) = {B_{m 0} (t), \dots, B_{m m} (t)}^{⊤}$ , and ${\bar{B}}_{m} (t) = {{\bar{B}}_{m 0} (t), \dots, {\bar{B}}_{m m} (t)}^{⊤}$ . Then $f_{m} (t; p) = p^{⊤} β_{m} (t)$ , $F_{m} (t; p) = p^{⊤} B_{m} (t) = 1 - S_{m} (t; p)$ , and $S_{m} (t; p) = p^{⊤} {\bar{B}}_{m} (t)$ . Replacing $f_{0}$ by $f_{m}$ in $ℓ (γ; z_{i})$ we have the approximate beta likelihood

\begin{aligned} ℓ_{m} (θ; z) & = (1 - δ) {γ^{⊤} \tilde{x} + \log f_{m} (y; p) - 2 \log [D_{m} (y | x, θ)]} \\ + δ \log [\frac{S_{m} (y_{1}; p)}{D_{m} (y_{1} | x, θ)} - \frac{S_{m} (y_{2}; p)}{D_{m} (y_{2} | x, θ)}], \end{aligned}

(8)

where

θ = (γ, p)

. Alternatively

\begin{aligned} ℓ_{m} (θ; z) & = γ^{⊤} \tilde{x} + (1 - δ) {\log f_{m} (y; p) - 2 \log [D_{m} (y | x, θ)]} \\ + δ \log [\frac{F_{m} (y_{2}; p)}{D_{m} (y_{2} | x, θ)} - \frac{F_{m} (y_{1}; p)}{D_{m} (y_{1} | x, θ)}] . \end{aligned}

(9)

Then the full loglikelihood

ℓ (γ)

can be approximated by

ℓ_{m} (θ) = \sum ℓ_{m} (θ; z_{i}) .

Theorem 2.2 implies that for any fixed

p

, the Hessian matrix of

ℓ_{m} (θ)

with respect to

γ

is negative definite for large sample size

n

under some weak conditions. Because

p_{0} = 1 - \sum_{j = 1}^{m} p_{j}

, for a fixed

γ

ℓ_{m} (θ)

is a function of

p = (p_{1}, \dots, p_{m})^{⊤}

on the simplex

C_{m} := {p \in R^{m} : p_{j} \geq 0, \sum_{j = 1}^{m} p_{j} \leq 1} .

About the asymptotic concavity of

ℓ_{m} (θ)

as a function of

p = (p_{1}, \dots, p_{m})^{⊤}

for a fixed

γ

, we have the following result.

Theorem 2.3

For a fixed $γ$ , if the sample size of Case $k$ data, for $k = 0, 1, 2$ , is either zero or large enough, then the approximate loglikelihood $ℓ_{m} (γ, p) = \sum_{i = 1}^{n} ℓ_{m} (γ, p; z_{i})$ as a function of $p = (p_{1}, \dots, p_{m})^{⊤}$ is concave for some $ϵ_{0} \in (0, 1)$ on the convex set

A_{m} (ϵ_{0}) = {p \in C_{m} : | \frac{f_{m} (t; p)}{f_{0} (t)} - 1 | \leq ϵ_{0} for all t \in [0, 1]} .

2.2. Maximum profile likelihood estimate of p

Let $\tilde{γ}$ be an efficient estimator of $γ$ using Turnbull’s method.¹⁶ For such $\tilde{γ}$ , we define the profile approximate loglikelihood ${\tilde{ℓ}}_{m} (p) = ℓ_{m} (\tilde{γ}, p)$ . With an optimal $m$ obtained by the change-point method,²⁶ the maximizer $\tilde{p} = ({\tilde{p}}_{0}, \dots, {\tilde{p}}_{m})^{⊤}$ of ${\tilde{ℓ}}_{m} (p)$ with respect to $p = (p_{0}, \dots, p_{m})^{⊤}$ subject to constraints in (7) is called the maximum approximate beta (profile) likelihood estimator of $p$ .

We can estimate $f_{0} (t) = f (t | x_{0})$ and $S_{0} (t) = S (t | x_{0})$ , respectively, by ${\tilde{f}}_{0} (t) = f_{m} (t; \tilde{p})$ and ${\tilde{S}}_{0} (t) = S_{m} (t; \tilde{p}) .$ Therefore with $\tilde{θ} = (\tilde{γ}, \tilde{p})$ we can estimate $f (t | x)$ and $S (t | x)$ , respectively, by ${\tilde{f}}_{B} (t | x) = f_{m} (t | x; \tilde{θ})$ and ${\tilde{S}}_{B} (t | x) = S_{m} (t | x; \tilde{θ})$ .

By (8) the derivatives of $ℓ_{m} (θ)$ w.r.t. $p$ are

\begin{aligned} Ψ (γ, p) & = (Ψ_{0} (γ, p), \dots, Ψ_{m} (γ, p))^{⊤} = \nabla_{p} ℓ_{m} (γ, p) \\ = \sum_{i = 1}^{n_{0}} [\frac{β_{m} (y_{i})}{f_{m} (y_{i}; p)} + \frac{2 (e^{γ^{⊤} {\tilde{x}}_{i}} - 1) {\bar{B}}_{m} (y_{i})}{1 + (e^{γ^{⊤} {\tilde{x}}_{i}} - 1) F_{m} (y_{i}; p)}] \\ + \sum_{i = n_{0} + 1}^{n} \frac{\nabla_{p} S_{m} (y_{1 i} | x_{i}; θ) - \nabla_{p} S_{m} (y_{2 i} | x_{i}; θ)}{S_{m} (y_{1 i} | x_{i}; θ) - S_{m} (y_{2 i} | x_{i}; θ)}, \end{aligned}

where

\nabla_{p} S_{m} (t | x; θ) = e^{γ^{⊤} \tilde{x}} {\bar{B}}_{m} (t) / [1 + (e^{γ^{⊤} \tilde{x}} - 1) F_{m} (t; p)]^{2} .

If we choose

x_{0}

so that

{\tilde{γ}}^{⊤} x_{0} \leq min {{\tilde{γ}}^{⊤} x_{i} : 1 \leq i \leq n}

, then

e^{γ^{⊤} {\tilde{x}}_{i}} \geq 1

i = 1, \dots, n

Remark 1
For each $j$ , $Ψ_{j} (γ, p) \geq 0$ provided that $e^{γ^{⊤} {\tilde{x}}_{i}} \geq 1$ , $i = 1, \dots, n$ . $Ψ_{j} (γ, p)$ depends on $p$ through $f_{m} (\cdot; p)$ and $S_{m} (\cdot; p)$ only.

Introducing Lagrange function $Λ (p) = ℓ_{m} (θ) - λ (\sum p_{j} - 1)$ , we have $\partial Λ (p) / \partial p = \partial ℓ_{m} (θ) / \partial p - λ 1 = Ψ (γ, p) - λ 1$ . Then $p^{⊤} \partial Λ (p) / \partial p = 0$ implies
$\begin{aligned} λ & = λ (γ, p) = p^{⊤} Ψ (γ, p) \\ = \sum_{i = 1}^{n_{0}} \frac{2 e^{γ^{⊤} {\tilde{x}}_{i}}}{D_{m} (y_{i} | x_{i}, θ)} - n_{0} + \sum_{i = n_{0} + 1}^{n} \frac{e^{γ^{⊤} {\tilde{x}}_{i}}}{S_{m} (y_{1 i} | x_{i}; θ) - S_{m} (y_{2 i} | x_{i}; θ)} [\frac{S_{m} (y_{1 i} | x_{i}; θ)}{D_{m} (y_{1 i} | x_{i}, θ)} - \frac{S_{m} (y_{2 i} | x_{i}; θ)}{D_{m} (y_{2 i} | x_{i}, θ)}] . \end{aligned}$

Assuming $e^{γ^{⊤} {\tilde{x}}_{i}} \geq 1$ , $i = 1, \dots, n$ , we have $λ (γ, p) > 0$ if $n_{0} > 0$ or $y_{1 i} < y_{2 i}$ for some $i > n_{0}$ .

The equation $\partial Λ (p) / \partial p = 0$ becomes
$\frac{1}{λ (γ, p)} Ψ (γ, p) = 1 .$
(10)

Therefore, starting with an initial guess $p^{[0]} \in S_{m}$ , a fixed-point iteration for solving equation (10) is
$p_{j}^{[s + 1]} = \frac{p_{j}^{[s]}}{λ (γ, p^{[s]})} Ψ_{j} (γ, p^{[s]}), j = 0, \dots, m .$
(11)
Remark 1 implies that $p^{[s]} \in S_{m}$ , $s = 1, \dots .$
2.3. Estimating both p and regression coefficients

For an optimal degree $m$ which can be found by the change-point method,²⁶ we call the maximizer $\hat{θ} = (\hat{γ}, \hat{p})$ of $ℓ_{m} (θ)$ the maximum approximate beta likelihood estimator (MABLE) of $θ = (γ, p)$ . The MABLE’s of $f (t | x)$ and $S (t | x)$ are, respectively,

\hat{f} (t | x) = f_{m} (t | x; \hat{θ}), \hat{S} (t | x) = S_{m} (t | x; \hat{θ}) .

(12)

We can find

(\hat{γ}, \hat{p})

by the following back-fitting algorithm.³¹

Step 0

Let $γ_{0}$ be an initial guess of $γ$ and $p^{(0)}$ in the interior of $S_{m}$ such as $p^{(0)} = u_{m} = (1, \dots, 1) / (m + 1)$ . Choose $x_{0}^{(0)}$ so that $γ_{0}^{⊤} x_{0}^{(0)} \leq min {γ_{0}^{⊤} x_{i} : 1 \leq i \leq n}$ . Use (11) with $\tilde{γ} = γ_{0}$ , $x_{0} = x_{0}^{(0)}$ , and starting point $p^{(0)}$ to get $p_{0} = \tilde{p}$ . Set $k = 0$

Step 1

Find the maximizer $γ_{k + 1}$ of $ℓ_{m} (γ, p_{k})$ using Newton-Raphson method or a quasi-Newton method if the Hessian matrix of the loglikelihood with respect to $γ$ is singular.

Step 2

Choose $x_{0}^{(k + 1)}$ so that $γ_{k + 1}^{⊤} x_{0}^{(k + 1)} \leq min {γ_{k + 1}^{⊤} x_{i} : 1 \leq i \leq n}$ . Use (11) with $\tilde{γ} = γ_{k + 1}$ and $x_{0} = x_{0}^{(k + 1)}$ to get $p_{k + 1} = \tilde{p}$ . Set $k = k + 1$ .

Step 3

Repeat Steps 1 and 2 until convergence. The final $(γ_{k}, p_{k})$ are taken as the MABLE of $θ = (γ, p)$ with baseline $x_{0} = x_{0}^{(k)}$

The use of the working baseline $x_{0}$ rather than the conventional choice $x_{0} = 0$ in the fixed-point iteration (11) ensures $p_{k} \in S_{m}$ without affecting $γ_{k}$ . The survival function at the actual baseline $x_{0}^{*} \neq x_{0}$ including $x_{0}^{*} = 0$ can be estimated based on relation (3). There are semiparametric methods implemented in R such as package icenReg.³² We can use the available efficient estimate of $γ$ as an initial guess for finding the maximum profile likelihood estimate of $p$ . The proposed method is implemented in R as a component of package mable which is available on CRAN.

With a $γ_{k}$ which is close enough to the true $γ$ , Theorem 2.3 ensures that $f_{m} (t; p_{k})$ is close to $f_{0} (t)$ . Theorem 2.2 in turn ensures that $γ_{k + 1}$ would get closer to the true $γ$ for large sample size. The convergence of the above back-fitting algorithm has also been tested in some simulation scenarios. The model flexibility is controlled by the degree $m$ . The change-point method²⁶ for choosing the model degree $m$ is designed to stop the algorithm before over-fitting. This method detects the elbow similar to that in the scree test³³ as a change-point of the log-likelihood ratios.

3. Large-sample properties

We need some assumptions and conditions to state the asymptotic results.

Assumption 1
The support $X$ of the covariate $X$ and $G$ are compact and for each $x_{0} \in X$ , $E (\tilde{X} {\tilde{X}}^{⊤})$ is positive definite, where $\tilde{X} = X - x_{0}$ .
Assumption 2
For each $x_{0} \in X$ the baseline density $f_{0}$ is positive and continuous on $(0, 1)$ , and there exist $ρ > 0$ which depends on $f_{0}$ only and $f_{m} (t; p_{0})$ for all $m \geq 1$ such that, uniformly in $t \in (0, 1)$ ,
$\frac{f_{m} (t; p_{0}) - f_{0} (t)}{f_{0} (t)} = O (m^{- ρ / 2}),$
where $p_{0} = (p_{0}, \dots, p_{m})$ .

Condition 0: Data are not censored.

Condition 1: Data are Case 1 censored. Examination time $U$ has density $g (u | x)$ with support $[τ_{l}, τ_{u}]$ , $0 < τ_{l} < τ_{u} < 1$ , and $0 < F_{0} (τ_{l}) < F_{0} (τ_{u}) < 1$ .

Condition 2: Data are Case 2 censored. Observed examination times $(U_{1}, U_{2})$ has density $g (u_{1}, u_{2} | x)$ with support $T = {(u_{1}, u_{2}) : τ_{l} \leq u_{1} < u_{2} \leq τ_{u}}$ , where $0 < τ_{l} < τ_{u} < 1$ , $0 < F_{0} (τ_{l}) < F_{0} (τ_{u}) < 1$ , and $P (U_{2} - U_{1} > η) = 1$ for some $η > 0$ .

The conditions on the examination times are similar to those used by other authors, for example, Huang and Wellner.¹⁶ Assumption 2 is fulfilled with $ρ = r + α$ by a large family of density functions(see Lemma 3.1 of Wang and Guan³⁴): $f_{0} (t) = t^{a} (1 - t)^{b} φ (t)$ for some nonnegative integers $a$ and $b$ , $φ (t) \geq b_{0} > 0$ , $φ \in C^{(r)} [0, 1]$ , $r \geq 0$ , and $φ^{(r)}$ is $α$ -Hölder continuous with $α \in (0, 1)$ . This is a generalization of a result of Lorentz³⁵ with $a = b = 0$ .
3.1. Some statistical distances

Define statistical distances

\begin{aligned} Δ_{00}^{2} (p; x_{0}) & = χ^{2} (f_{m} (\cdot; p) ‖ f_{0} (\cdot)) = \int_{0}^{1} {[\frac{f_{m} (t; p)}{f_{0} (t)} - 1]}^{2} f_{0} (t) d t, \\ Δ_{01}^{2} (p; x_{0}) & = \int_{0}^{1} [F_{m} (t; p) - F_{0} (t)]^{2} f_{0} (t) d t, \\ Δ_{0}^{2} (θ; x_{0}) & = \int_{X} \int_{0}^{1} {[\frac{f_{m} (t | x; θ)}{f (t | x; γ_{0}, f_{0})} - 1]}^{2} f (t | x; γ_{0}) d t d H (x), \\ Δ_{1}^{2} (θ; x_{0}) & = \int_{X} \int_{τ_{l}}^{τ_{u}} \frac{{[F_{m} (u | x; θ) - F (u | x; γ_{0})]}^{2}}{F (u | x; γ_{0}) S (u | x; γ_{0})} g (u | x) d u d H (x), \\ Δ_{10}^{2} (p; x_{0}) & = \int_{X} \int_{τ_{l}}^{τ_{u}} \frac{[F_{m} (u; p) - F_{0} (u)]^{2}}{F_{0} (u) S_{0} (u)} d G (u | x) d H (x), \\ Δ_{2}^{2} (θ; x_{0}) & = \int_{X} \iint_{T} \frac{[V_{m} (u | x)]^{⊤} M (u | x) V_{m} (u | x)}{F (u_{2} | x; γ_{0}) - F (u_{1} | x; γ_{0})} g (u | x) d u_{1} d u_{2} d H (x), \\ Δ_{20}^{2} (p; x_{0}) & = \int_{X} \iint_{T} {\frac{[F_{m} (u_{1}; p) - F_{0} (u_{1})]^{2}}{F_{0} (u_{1})} + \frac{[F_{m} (u_{2}; p) - F_{0} (u_{2})]^{2}}{S_{0} (u_{2})}} \\ \cdot g (u | x) d u_{1} d u_{2} d H (x), \end{aligned}

where

u = (u_{1}, u_{2})

V_{m} (u | x) = [F_{m} (u_{1} | x; θ) - F (u_{1} | x; γ_{0}), F_{m} (u_{2} | x; θ) - F (u_{2} | x; γ_{0})]^{⊤}

and

It is easy to see that

M (u | x)

is positive definite whenever

0 < F (u_{1} | x; γ_{0}) < F (u_{2} | x; γ_{0}) < 1

Theorem 3.1
Suppose that Assumptions 1 and 2 are satisfied. About the relations of these distances we have the following results.
$\begin{aligned} Δ_{01}^{2} (p, x_{0}) & \leq 0.5 Δ_{00}^{2} (p, x_{0}), \end{aligned}$
(13)

$\begin{aligned} Δ_{k}^{2} (θ; x_{0}) & \leq Δ_{0}^{2} (θ; x_{0}), k = 1, 2. \end{aligned}$
(14)
Under Condition $k$ , $k = 0, 1, 2$ ,
$\begin{aligned} Δ_{k}^{2} (θ; x_{0}) & \leq C_{1 k} [‖ γ - γ_{0} ‖^{2} + Δ_{k 0}^{2} (p, x_{0})], \end{aligned}$
(15)

$\begin{aligned} Δ_{k}^{2} (θ; x_{0}) & \geq C_{2 k} ‖ γ - γ_{0} ‖^{2} - C_{2 k}^{'} Δ_{k 0}^{2} (p, x_{0}), \end{aligned}$
(16)

$\begin{aligned} Δ_{k}^{2} (θ; x_{0}) & \geq C_{3 k} Δ_{k 0}^{2} (p, x_{0}) - C_{3 k}^{'} ‖ γ - γ_{0} ‖^{2}, \end{aligned}$
(17)
where $C_{i k} > 0$ and $C_{i k}^{'} > 0$ , $i = 1, 2, 3$ , are constants and additionally (16) requires that $‖ γ - γ_{0} ‖^{2} \leq 1 / λ^{}$ , where
$λ^{} = max {λ (x) : λ (x) is the largest eigenvalue of \tilde{x} {\tilde{x}}^{⊤}, x \in X} .$

3.2. Consistency

Define $r_{n} (ϵ) = (\log \log n)^{1 + ϵ} / n$ , $ϵ > 0$ , $B_{d} (n, ϵ) = {γ \in G : ‖ γ - γ_{0} ‖^{2} \leq r_{n} (ϵ)}$ , and $D_{m}^{(k)} (n, ϵ) = {p \in A_{m} (ϵ_{0}) : Δ_{k 0}^{2} (p; x_{0}) \leq r_{n} (ϵ)}, k = 0, 1, 2.$

Lemma 3.2
The sets $B_{d} (n, ϵ)$ and $D_{m}^{(k)} (n, ϵ)$ , $k = 0, 1, 2$ , are convex, so are $B_{d} (n, ϵ) \times D_{m}^{(k)} (n, ϵ)$ for $k = 0, 1, 2$ .
Theorem 3.3
Suppose that Assumptions 1 and 2 are satisfied with degree $m \sim C (n / \log \log n)^{1 / ρ}$ and Condition $k$ is fulfilled for a $k = 0, 1, 2$ . (i) If $‖ γ - γ_{0} ‖^{2} \leq r_{n} (ϵ)$ for an $ϵ \in (0, 1 / 2)$ , then for any $ϵ^{'} \in (ϵ, 1 / 2)$ and large $n$ the maximizer $\tilde{p} (γ)$ of $ℓ_{m} (γ, p)$ almost surely satisfies $Δ_{k 0}^{2} (p; x_{0}) < r_{n} (ϵ^{'})$ . (ii) Conversely, if $Δ_{k 0}^{2} (p; x_{0}) \leq r_{n} (ϵ)$ for an $ϵ \in (0, 1 / 2)$ , then for any $ϵ^{'} \in (ϵ, 1 / 2)$ and large $n$ the maximizer $\tilde{γ} (p)$ of $ℓ_{m} (γ, p)$ almost surely satisfies $‖ γ - γ_{0} ‖^{2} < r_{n} (ϵ^{'})$ .
Remark 2
If $\tilde{γ}$ is an efficient estimator of $γ$ such as the Turnbull’s estimator, then the maximum approximate beta profile likelihood estimator $\tilde{p} = \tilde{p} (\tilde{γ})$ satisfies $Δ_{k 0}^{2} (\tilde{p}; x_{0}) < r_{n} (ϵ^{'})$ . Moreover, if $\hat{θ} = (\hat{γ}, \hat{p})$ can be obtained using the algorithm of Section 2.3, then for each $k = 0, 1, 2$ , and any $ϵ \in (0, 1 / 2)$ , under Condition $k$ , we have $‖ \hat{γ} - γ ‖ = O (r_{n} (ϵ))$ and $Δ_{k 0}^{2} (\hat{p}; x_{0}) = O (r_{n} (ϵ))$ .
Remark 3
Theorem 3.3 implies that the proposed estimators posses an almost parametric $\sqrt{n}$ -consistency, while Turnbull’s approach is typically $n^{1 / 3}$ -consistent.
3.3. Asymptotic normality

Theorem 3.4
Suppose that Assumptions 1 and 2 are fulfilled, and $m = m_{n}$ satisfies $n^{1 / 2} m^{- ρ / 2} = o (1)$ . Let $\tilde{γ} = \tilde{γ} (p_{0})$ be the maximizer of $ℓ_{m} (γ, p_{0})$ . For $k = 0, 1, 2$ , under Condition k, $\sqrt{n} (\tilde{γ} - γ_{0})$ converges in distribution to $N ((0, I_{k}^{- 1})$ as $n \to \infty$ , where $x_{0}$ satisfies $γ_{0}^{⊤} x_{0} < min {γ_{0}^{⊤} x_{i} : 1 \leq i \leq n}$ , and $I_{k} = E [J_{k} (X; γ_{0}) \tilde{X} {\tilde{X}}^{⊤}]$ .

From the proof of this theorem we have approximation
$I_{k} \approx \frac{1}{n} \sum_{i = 1}^{n} {(F_{m} S_{m}) (y_{1 i} | x_{i}; \hat{γ}, \hat{p}) + (F_{m} S_{m}) (y_{2 i} | x_{i}; \hat{γ}, \hat{p})} {\tilde{x}}_{i} {\tilde{x}}_{i}^{⊤},$
where all $δ_{i} = 0$ if $k = 0$ and all $δ_{i} = 1$ if $k = 1, 2$ .

The asymptotic normality of $\hat{γ}$ can also be obtained under some additional regularity conditions imposed on the mixture model $f_{m} (t | x; p)$ which requires a known $m$ and a known positive subvector of $p$ . Clearly, the variation of $\hat{p}$ also contributes to the variance of $\hat{γ}$ .
4. Simulation

In this simulation we shall use exponential and gamma distributions as the baseline $F_{0}$ . The covariates $X = (X_{1}, X_{2})^{⊤}$ with coefficients $γ = (γ_{1}, γ_{2})^{⊤} = (0.5, - 0.5)^{⊤}$ and baseline $x = 0$ are generated as independent $X_{1}$ and $X_{2}$ , where $X_{1} = (U_{1} + \dots + U_{k}) / k$ with $U_{1}, \dots, U_{k}$ being IID uniform $(- a_{1}, a_{1})$ and $X_{2} = a_{2} V + a_{3}$ with $V$ being binomial $(ν, 0.5)$ . Two models were used to generate $X$ . Model I: $k = 1, ν = 1, a_{1} = 1, a_{2} = 2, a_{3} = - 1$ ; Model II: $k = 3, ν = 2, a_{1} = 3, a_{2} = 2, a_{3} = - 2$ .

Table 1.
Simulation results with exponential baseline, sample size $n = 50, 100, 500$ , and case $k$ censoring ( $k = 1, 2, 5$ ): (i) the bias( $\times 10$ ) and the standard deviation( $\times 10$ ) in parentheses of parametric estimator $\overset{ˇ}{γ}$ , the proposed estimator $\hat{γ}$ , and the semiparametric estimator $\tilde{γ}$ ; (ii) the root mean integrated squared errors( $\times 10$ ) of parametric estimators ${\overset{ˇ}{f}}_{0}$ and ${\overset{ˇ}{S}}_{0}$ , the proposed estimators ${\hat{f}}_{0}$ and ${\hat{S}}_{0}$ , and the semiparametric estimator ${\tilde{S}}_{0}$ ; and (iii) the convergence percentage C(%). For case 1 censoring, left-censoring (lc) rates are 25%, 50%, and 75%.

$k$ (lc) $n$ ${\overset{ˇ}{γ}}_{1}$ ${\overset{ˇ}{γ}}_{2}$ ${\hat{γ}}_{1}$ ${\hat{γ}}_{2}$ ${\tilde{γ}}_{1}$ ${\tilde{γ}}_{2}$ ${\overset{ˇ}{f}}_{0}$ ${\hat{f}}_{0}$ ${\overset{ˇ}{S}}_{0}$ ${\hat{S}}_{0}$ ${\tilde{S}}_{0}$ C(%)

Model I

1(25) 50 0.41(7.24) −0.31(4.73) 1.32(8.56) −0.76(6.18) 4.99(25.55) −7.15(19.85) 1.69 2.66 2.57 4.68 4.96 94.1

1(25) 100 0.12(4.87) −0.01(2.79) 0.65(5.28) −0.30(3.13) 1.79(6.32) −1.75(4.95) 1.09 1.99 1.15 3.10 3.44 96.3

1(25) 500 0.01(2.05) 0.00(1.20) 0.06(2.12) 0.04(1.28) 0.49(2.24) −0.47(1.37) 0.48 1.05 0.47 1.70 2.14 100

1(50) 50 0.38(6.90) −0.07(3.73) 1.03(7.95) −0.35(4.34) 4.89(27.86) −6.97(18.65) 1.26 2.60 1.21 1.81 2.58 94.1

1(50) 100 0.20(4.29) −0.02(2.50) 0.57(4.64) −0.11(2.70) 1.91(5.59) −1.70(4.35) 0.82 1.92 0.80 1.08 1.81 96.3

1(50) 500 0.02(1.87) 0.10(1.11) −0.21(2.03) 0.36(1.24) 0.51(2.09) −0.37(1.31) 0.35 1.06 0.34 0.60 1.01 100

1(75) 50 0.60(9.78) −0.31(5.11) 1.79(10.68) −0.80(5.84) 37.59(431.71) −30.81(226.15) 1.54 5.09 1.26 1.62 3.22 94.1

1(75) 100 0.22(5.95) −0.05(3.27) 0.74(6.27) −0.12(3.55) 2.50(8.84) −3.13(9.86) 0.94 2.51 0.84 1.11 2.13 96.3

1(75) 500 0.08(2.29) 0.00(1.39) −0.35(2.42) 0.40(1.54) 0.71(2.61) −0.59(1.66) 0.40 1.33 0.38 0.57 1.11 100

2 50 0.51(6.39) −0.59(3.84) 0.77(6.57) −0.68(4.07) 1.49(8.28) −2.56(10.14) 1.03 2.81 0.92 1.34 2.17 94.1

2 100 0.09(4.09) −0.19(2.39) 0.09(4.07) −0.11(2.46) 0.66(4.66) −0.72(2.90) 0.65 1.99 0.61 0.92 1.58 96.3

2 500 0.07(1.81) 0.01(1.06) 0.11(1.84) −0.01(1.08) 0.21(1.89) −0.13(1.12) 0.29 1.22 0.28 0.47 0.85 100

5 50 0.34(5.07) −0.23(2.92) 0.52(5.13) −0.20(2.99) 0.74(5.52) −0.59(3.23) 0.88 2.02 0.84 1.22 1.70 94.1

5 100 0.27(3.22) −0.06(1.92) 0.15(3.38) 0.05(2.02) 0.53(3.45) −0.28(2.07) 0.61 1.53 0.58 0.82 1.26 96.3

5 500 0.02(1.41) 0.01(0.85) 0.05(1.45) −0.01(0.86) 0.09(1.46) −0.05(0.88) 0.26 0.78 0.25 0.41 0.69 100

Model II

1(25) 50 0.24(3.56) −0.31(3.09) 1.13(4.56) −1.21(4.16) 5.32(20.70) −7.79(32.00) 1.68 3.15 1.95 4.20 4.43 98.8

1(25) 100 0.08(2.32) −0.06(2.01) 0.59(2.83) −0.53(2.44) 2.08(3.68) −2.04(3.28) 1.18 2.49 1.28 3.58 3.52 100

1(25) 500 −0.01(0.96) 0.04(0.82) 0.09(1.08) −0.04(0.95) 0.53(1.15) −0.49(1.02) 0.49 1.23 0.48 1.58 2.04 100

1(50) 50 0.04(3.18) −0.15(2.76) 0.80(4.11) −0.83(3.60) 4.79(16.51) −5.40(12.53) 1.21 3.10 1.17 1.71 2.49 98.8

1(50) 100 0.02(2.14) 0.00(1.76) 0.42(2.58) −0.34(2.15) 1.87(3.21) −1.85(2.83) 0.85 2.33 0.83 1.20 1.86 100

1(50) 500 −0.02(0.89) 0.05(0.79) 0.05(1.03) 0.01(0.91) 0.55(1.10) −0.49(0.99) 0.36 1.21 0.35 0.54 1.00 100

1(75) 50 0.64(13.83) −0.54(12.64) 1.40(5.87) −1.16(5.13) 33.88(181.12) −38.05(202.43) 2.12 4.49 1.30 1.78 3.05 98.8

1(75) 100 0.12(2.71) −0.19(2.38) 0.68(3.40) −0.65(2.95) 2.79(5.14) −3.15(6.42) 0.94 3.04 0.84 1.26 2.21 100

1(75) 500 0.03(1.09) 0.03(1.00) 0.02(1.32) 0.06(1.16) 0.66(1.38) −0.58(1.29) 0.40 1.55 0.38 0.60 1.15 100

2 50 0.57(3.24) −0.32(2.68) 0.95(3.62) −0.59(2.99) 1.88(4.86) −1.66(6.45) 1.09 3.29 0.97 1.39 2.17 98.8

2 100 0.14(2.10) −0.23(1.89) 0.35(2.26) −0.37(2.01) 0.68(2.46) −0.77(2.27) 0.70 2.42 0.66 0.96 1.61 100

2 500 0.00(0.88) −0.02(0.76) 0.06(0.95) −0.05(0.79) 0.15(0.94) −0.17(0.81) 0.30 1.41 0.29 0.48 0.86 100

5 50 0.15(2.44) −0.23(2.17) 0.47(2.66) −0.47(2.38) 0.58(2.80) −0.65(2.51) 0.92 2.43 0.87 1.14 1.66 98.8

5 100 0.00(1.69) −0.12(1.44) 0.17(1.78) −0.25(1.53) 0.23(1.82) −0.35(1.57) 0.62 1.78 0.58 0.77 1.23 100

5 500 0.02(0.71) −0.01(0.64) 0.07(0.74) −0.04(0.66) 0.08(0.74) −0.07(0.67) 0.28 1.03 0.27 0.39 0.67 100

$k$ (lc)	$n$	${\overset{ˇ}{γ}}_{1}$	${\overset{ˇ}{γ}}_{2}$	${\hat{γ}}_{1}$	${\hat{γ}}_{2}$	${\tilde{γ}}_{1}$	${\tilde{γ}}_{2}$	${\overset{ˇ}{f}}_{0}$	${\hat{f}}_{0}$	${\overset{ˇ}{S}}_{0}$	${\hat{S}}_{0}$	${\tilde{S}}_{0}$	C(%)
Model I
1(25)	50	0.41(7.24)	−0.31(4.73)	1.32(8.56)	−0.76(6.18)	4.99(25.55)	−7.15(19.85)	1.69	2.66	2.57	4.68	4.96	94.1
1(25)	100	0.12(4.87)	−0.01(2.79)	0.65(5.28)	−0.30(3.13)	1.79(6.32)	−1.75(4.95)	1.09	1.99	1.15	3.10	3.44	96.3
1(25)	500	0.01(2.05)	0.00(1.20)	0.06(2.12)	0.04(1.28)	0.49(2.24)	−0.47(1.37)	0.48	1.05	0.47	1.70	2.14	100
1(50)	50	0.38(6.90)	−0.07(3.73)	1.03(7.95)	−0.35(4.34)	4.89(27.86)	−6.97(18.65)	1.26	2.60	1.21	1.81	2.58	94.1
1(50)	100	0.20(4.29)	−0.02(2.50)	0.57(4.64)	−0.11(2.70)	1.91(5.59)	−1.70(4.35)	0.82	1.92	0.80	1.08	1.81	96.3
1(50)	500	0.02(1.87)	0.10(1.11)	−0.21(2.03)	0.36(1.24)	0.51(2.09)	−0.37(1.31)	0.35	1.06	0.34	0.60	1.01	100
1(75)	50	0.60(9.78)	−0.31(5.11)	1.79(10.68)	−0.80(5.84)	37.59(431.71)	−30.81(226.15)	1.54	5.09	1.26	1.62	3.22	94.1
1(75)	100	0.22(5.95)	−0.05(3.27)	0.74(6.27)	−0.12(3.55)	2.50(8.84)	−3.13(9.86)	0.94	2.51	0.84	1.11	2.13	96.3
1(75)	500	0.08(2.29)	0.00(1.39)	−0.35(2.42)	0.40(1.54)	0.71(2.61)	−0.59(1.66)	0.40	1.33	0.38	0.57	1.11	100
2	50	0.51(6.39)	−0.59(3.84)	0.77(6.57)	−0.68(4.07)	1.49(8.28)	−2.56(10.14)	1.03	2.81	0.92	1.34	2.17	94.1
2	100	0.09(4.09)	−0.19(2.39)	0.09(4.07)	−0.11(2.46)	0.66(4.66)	−0.72(2.90)	0.65	1.99	0.61	0.92	1.58	96.3
2	500	0.07(1.81)	0.01(1.06)	0.11(1.84)	−0.01(1.08)	0.21(1.89)	−0.13(1.12)	0.29	1.22	0.28	0.47	0.85	100
5	50	0.34(5.07)	−0.23(2.92)	0.52(5.13)	−0.20(2.99)	0.74(5.52)	−0.59(3.23)	0.88	2.02	0.84	1.22	1.70	94.1
5	100	0.27(3.22)	−0.06(1.92)	0.15(3.38)	0.05(2.02)	0.53(3.45)	−0.28(2.07)	0.61	1.53	0.58	0.82	1.26	96.3
5	500	0.02(1.41)	0.01(0.85)	0.05(1.45)	−0.01(0.86)	0.09(1.46)	−0.05(0.88)	0.26	0.78	0.25	0.41	0.69	100
Model II
1(25)	50	0.24(3.56)	−0.31(3.09)	1.13(4.56)	−1.21(4.16)	5.32(20.70)	−7.79(32.00)	1.68	3.15	1.95	4.20	4.43	98.8
1(25)	100	0.08(2.32)	−0.06(2.01)	0.59(2.83)	−0.53(2.44)	2.08(3.68)	−2.04(3.28)	1.18	2.49	1.28	3.58	3.52	100
1(25)	500	−0.01(0.96)	0.04(0.82)	0.09(1.08)	−0.04(0.95)	0.53(1.15)	−0.49(1.02)	0.49	1.23	0.48	1.58	2.04	100
1(50)	50	0.04(3.18)	−0.15(2.76)	0.80(4.11)	−0.83(3.60)	4.79(16.51)	−5.40(12.53)	1.21	3.10	1.17	1.71	2.49	98.8
1(50)	100	0.02(2.14)	0.00(1.76)	0.42(2.58)	−0.34(2.15)	1.87(3.21)	−1.85(2.83)	0.85	2.33	0.83	1.20	1.86	100
1(50)	500	−0.02(0.89)	0.05(0.79)	0.05(1.03)	0.01(0.91)	0.55(1.10)	−0.49(0.99)	0.36	1.21	0.35	0.54	1.00	100
1(75)	50	0.64(13.83)	−0.54(12.64)	1.40(5.87)	−1.16(5.13)	33.88(181.12)	−38.05(202.43)	2.12	4.49	1.30	1.78	3.05	98.8
1(75)	100	0.12(2.71)	−0.19(2.38)	0.68(3.40)	−0.65(2.95)	2.79(5.14)	−3.15(6.42)	0.94	3.04	0.84	1.26	2.21	100
1(75)	500	0.03(1.09)	0.03(1.00)	0.02(1.32)	0.06(1.16)	0.66(1.38)	−0.58(1.29)	0.40	1.55	0.38	0.60	1.15	100
2	50	0.57(3.24)	−0.32(2.68)	0.95(3.62)	−0.59(2.99)	1.88(4.86)	−1.66(6.45)	1.09	3.29	0.97	1.39	2.17	98.8
2	100	0.14(2.10)	−0.23(1.89)	0.35(2.26)	−0.37(2.01)	0.68(2.46)	−0.77(2.27)	0.70	2.42	0.66	0.96	1.61	100
2	500	0.00(0.88)	−0.02(0.76)	0.06(0.95)	−0.05(0.79)	0.15(0.94)	−0.17(0.81)	0.30	1.41	0.29	0.48	0.86	100
5	50	0.15(2.44)	−0.23(2.17)	0.47(2.66)	−0.47(2.38)	0.58(2.80)	−0.65(2.51)	0.92	2.43	0.87	1.14	1.66	98.8
5	100	0.00(1.69)	−0.12(1.44)	0.17(1.78)	−0.25(1.53)	0.23(1.82)	−0.35(1.57)	0.62	1.78	0.58	0.77	1.23	100
5	500	0.02(0.71)	−0.01(0.64)	0.07(0.74)	−0.04(0.66)	0.08(0.74)	−0.07(0.67)	0.28	1.03	0.27	0.39	0.67	100

Table 2.

Simulation results with gamma baseline, sample size $n = 50, 100, 500$ , and case $k$ censoring ( $k = 1, 2, 5$ ): (i) the bias( $\times 10$ ) and the standard deviation( $\times 10$ ) in parentheses of parametric estimator $\overset{ˇ}{γ}$ , the proposed estimator $\hat{γ}$ , and the semiparametric estimator $\tilde{γ}$ ; (ii) the root mean integrated squared errors( $\times 10$ ) of parametric estimators ${\overset{ˇ}{f}}_{0}$ and ${\overset{ˇ}{S}}_{0}$ , the proposed estimators ${\hat{f}}_{0}$ and ${\hat{S}}_{0}$ , and the semiparametric estimator ${\tilde{S}}_{0}$ ; and (iii) the convergence percentage C(%). For case 1 censoring, left-censoring (lc) rates are 25%, 50%, and 75%.

$k$ (lc)	$n$	${\overset{ˇ}{γ}}_{1}$	${\overset{ˇ}{γ}}_{2}$	${\hat{γ}}_{1}$	${\hat{γ}}_{2}$	${\tilde{γ}}_{1}$	${\tilde{γ}}_{2}$	${\overset{ˇ}{f}}_{0}$	${\hat{f}}_{0}$	${\overset{ˇ}{S}}_{0}$	${\hat{S}}_{0}$	${\tilde{S}}_{0}$	C(%)
Model I
1(25)	50	1.26(9.39)	−1.09(6.53)	0.57(8.90)	−0.06(5.51)	5.59(35.89)	−9.02(24.40)	1.09	0.86	5.77	7.21	8.15	98.1
1(25)	100	0.60(5.40)	−0.47(3.35)	0.31(5.74)	0.06(3.24)	1.94(6.50)	−2.02(5.91)	0.69	0.58	3.30	3.95	5.39	99.5
1(25)	500	0.04(2.21)	−0.10(1.35)	−0.03(2.33)	0.58(1.26)	0.50(2.37)	−0.54(1.44)	0.31	0.59	1.50	3.54	3.13	100
1(50)	50	1.16(8.19)	−0.90(4.83)	0.27(7.75)	0.31(4.09)	4.81(27.80)	−7.23(19.10)	0.91	0.74	3.37	3.65	5.81	98.1
1(50)	100	0.54(4.67)	−0.34(2.84)	0.00(4.75)	0.38(2.63)	1.83(5.55)	−1.67(4.37)	0.57	0.53	2.24	2.38	4.04	99.5
1(50)	500	0.17(2.02)	−0.04(1.24)	−0.02(2.04)	0.30(1.20)	0.60(2.15)	−0.45(1.33)	0.24	0.30	0.99	1.22	2.25	100
1(75)	50	1.64(14.20)	−1.77(8.38)	−0.38(8.85)	0.55(4.74)	22.25(215.05)	−22.88(119.00)	1.35	1.09	4.34	4.32	6.98	98.1
1(75)	100	0.63(5.80)	−0.47(3.57)	−0.37(5.43)	0.64(3.10)	2.31(7.39)	−2.56(7.55)	0.71	0.72	2.68	2.90	4.82	99.5
1(75)	500	0.03(2.34)	−0.10(1.43)	−0.32(2.32)	0.30(1.37)	0.51(2.53)	−0.57(1.57)	0.28	0.35	1.13	1.32	2.56	100
2	50	0.81(6.61)	−0.47(3.66)	1.38(7.14)	−0.38(3.37)	1.79(8.04)	−1.44(5.84)	0.65	0.69	2.80	4.05	4.87	98.1
2	100	0.18(4.00)	−0.24(2.39)	1.17(4.56)	−0.20(2.16)	0.62(4.44)	−0.61(2.63)	0.43	0.60	1.92	3.49	3.69	99.5
2	500	−0.05(1.75)	−0.07(1.00)	0.08(1.74)	−0.01(0.98)	0.09(1.82)	−0.20(1.04)	0.18	0.22	0.82	1.01	2.03	100
5	50	0.40(5.58)	−0.34(3.06)	1.15(6.18)	−0.50(3.04)	0.72(6.07)	−0.67(3.38)	0.54	0.49	2.58	3.26	6.03	98.1
5	100	0.26(3.60)	−0.16(2.08)	1.10(3.93)	−0.30(2.05)	0.42(3.76)	−0.31(2.20)	0.37	0.38	1.79	2.39	4.73	99.5
5	500	−0.01(1.51)	−0.04(0.90)	0.21(1.49)	−0.05(0.90)	0.05(1.53)	−0.09(0.91)	0.15	0.21	0.77	1.16	2.97	100
Model II
1(25)	50	1.37(5.35)	−1.41(4.68)	0.79(4.62)	−0.73(4.07)	6.33(23.85)	−8.74(33.36)	1.12	0.97	5.61	7.11	7.83	99.8
1(25)	100	0.55(2.87)	−0.46(2.50)	0.41(2.94)	−0.23(2.44)	2.24(3.82)	−2.17(3.91)	0.69	0.68	3.44	4.48	5.49	100
1(25)	500	0.19(1.19)	−0.12(1.00)	0.12(1.23)	0.12(1.02)	0.72(1.30)	−0.64(1.11)	0.29	0.41	1.43	2.57	3.03	100
1(50)	50	0.92(4.37)	−1.08(3.99)	0.23(3.96)	−0.24(3.41)	4.80(16.54)	−5.42(12.62)	0.90	0.84	3.39	3.63	5.61	99.8
1(50)	100	0.37(2.51)	−0.32(2.12)	0.08(2.51)	0.06(2.04)	1.88(3.20)	−1.86(2.83)	0.55	0.61	2.29	2.49	4.18	100
1(50)	500	0.05(0.99)	−0.09(0.92)	−0.07(1.00)	0.06(0.91)	0.54(1.07)	−0.58(1.00)	0.23	0.37	0.97	1.19	2.25	100
1(75)	50	2.71(19.71)	−5.22(97.79)	0.14(4.78)	0.06(3.97)	21.75(154.96)	−26.09(155.60)	1.34	1.19	4.34	4.39	6.81	99.8
1(75)	100	0.57(3.06)	−0.50(2.71)	−0.04(3.01)	0.21(2.58)	2.33(4.17)	−2.33(4.10)	0.69	0.83	2.71	2.99	4.81	100
1(75)	500	−0.02(1.21)	−0.05(1.08)	−0.26(1.23)	0.21(1.07)	0.53(1.35)	−0.59(1.20)	0.29	0.50	1.18	1.45	2.57	100
2	50	0.69(3.27)	−0.52(2.75)	0.91(3.25)	−0.55(2.67)	1.59(4.20)	−1.41(3.53)	0.66	0.63	2.87	3.25	5.04	99.8
2	100	0.21(2.15)	−0.28(1.90)	0.40(2.15)	−0.33(1.87)	0.67(2.46)	−0.72(2.17)	0.45	0.48	2.06	2.38	3.77	100
2	500	0.06(0.87)	−0.05(0.77)	0.12(0.86)	−0.07(0.77)	0.21(0.91)	−0.20(0.80)	0.20	0.30	0.92	1.23	2.08	100
5	50	0.33(2.79)	−0.51(2.50)	0.64(2.75)	−0.68(2.48)	0.64(3.06)	−0.81(2.76)	0.58	0.57	2.85	3.79	6.48	99.8
5	100	0.13(1.78)	−0.23(1.61)	0.36(1.79)	−0.35(1.61)	0.31(1.89)	−0.40(1.71)	0.38	0.41	1.90	2.51	4.88	100
5	500	0.07(0.78)	0.02(0.68)	0.10(0.78)	0.00(0.69)	0.12(0.79)	−0.04(0.70)	0.16	0.31	0.82	1.95	3.05	100

We shall compare the proposed method with the parametric method according to the known simulation models and the semiparametric method based on Turnbull’s method. We used ic_par() and ic_sp() of R package icenReg,³² respectively, for the implementation of the parametric and the semiparametric methods. Specifically, for the regression coefficients $γ$ we compare the proposed estimator $\hat{γ}$ with the parametric estimator $\overset{ˇ}{γ}$ and the semiparametric estimator $\tilde{γ}$ in terms of the root mean square errors; for the baseline density we compare the proposed estimator ${\hat{f}}_{0} (\cdot) = f_{m} (\cdot | 0; \hat{θ})$ with the parametric estimator ${\overset{ˇ}{f}}_{0}$ with parameters calculated by ic_par() in terms of the square root of the mean integrated squared errors; for the baseline survival function we compare the proposed estimator ${\hat{S}}_{0} (\cdot) = S_{m} (\cdot | 0; \hat{θ})$ with parametric estimator ${\overset{ˇ}{S}}_{0}$ with parameters calculated by ic_par(), and the semiparametric estimator ${\tilde{S}}_{0}$ returned by ic_sp() in terms of the square root of the mean integrated squared errors. Because the step-function estimate ${\tilde{S}}_{0}$ is not necessarily unique we estimate $S_{0}$ as in R package interval³⁶ using the interpolation method for the purpose of comparison.

The event times are subject to Case $k$ interval censoring, $k = 1, 2, 5$ . The examination times of the Case $k$ interval censored data ( $k > 1$ ) are generated in the similar way as does the R function simIC_web of package icenReg. Sample sizes are 50 and 100. The number of Monte Carlo runs is 1000.

For a given covariate $X = x$ , we generate independent times $T$ and $T^{*}$ from the same conditional distribution $F (t | x)$ . We use $T$ as the actual event time and use $E = r T^{*}$ as the examination time with a chosen $r > 0$ so that the probability of left censoring is $P (T < r T^{*}) = p$ . In this simulation the left-censoring rates are 25%, 50%, and 75%. Because the loglikelihood is only asymptotically concave, the algorithm may not converge. The last columns of Tables 1 and 2 show the percentage the runs which converged.

For the exponential baseline with rate $β = 1$ and truncation interval $[0, 15]$ , we choose duration=5/exams for exams=2,5 as the arguments of simCS_weib(). We obtain the simulation results as shown in Table 1.

For the gamma baseline with shape $α = 5$ and rate $β = 0.5$ and truncation interval $[0, 50]$ we choose duration=20/exams for exams=2,5 as the arguments of simCS_weib(). We obtain the simulation results as shown in Table 2.

An examplified comparison of the proposed method with the parametric and semiparametric methods is given in Figure 1.

Figure 1.

Estimated survival(left) and density(right) curves based on a simulated current status dataset using model II with gamma baseline as shown in Table 2 and 75% left-censoring rate.

The simulation shows that (i) the proposed method performs much better than the semiparametric one; (ii) although the parametric method outperforms the proposed one the latter is closer to the former than to the semiparametric one in most cases; (iii) in some cases of the current status data the proposed method even outperforms the parametric method possibly due to the inflexibility of the latter; (iv) the current status data with 50% left-censoring produce better estimation for the specific examination time as generated in this simulation; (v) as sample size increases the percentage of convergent Monte Carlo runs increases.

5. Application to the HIV infection data

Several authors^37–39 studied the times to HIV infection for Danish homosexual men from two cities in Denmark. In this study among 297 who were tested at least once for HIV-antibody positivity on six different dates 65 have been tested positive. The dataset is contained in R package Epi.⁴⁰ The covariates that we used are us: whether the person has or not visited the USA, age: the age of the subject at entry, and pyr: the annual number of sexual partners. The event time is measured in days since start of the study. The estimated survival curves with us=0 and us=1, and average values of age and pyr using the proposed method and the semiparemetric method implemented in R package icenReg are shown in Figure 2. The semiparametric estimates for PH and PO models and the proposed estimates are shown in the top panels of Figure 2. The proposed estimates are also compared with those obtained using PH and AFT models^4,29 in the lower panels. We chose $τ = 5500$ for baseline $x_{0} = (us, age, pyr) = (0, 72, 0)$ , where 72 is the largest value of age. The proposed method gave coefficient estimates (0.596, -0.021, 0.012) with standard errors (0.297, 0.005, 0.005) and model degree $m = 7$ chosen by the change-point method.²⁶ Using the method of Guan²⁹ for PH model we obtained coefficient estimates (0.486, -0.024, 0.012) with standard errors (0.26, 0.004, 0.004) and model degree $m = 22$ . These are very close to the semiparametric estimates for PH and PO models as shown in Figure 2. Guan⁴ obtained the estimated coefficients in the AFT model (0.542, 0.015, 0.004) with standard errors (0.193, 0.006, 0.003) and model degree $m = 12$ . All the models agree the most on the effect of US visit.

Figure 2.

Estimated survival curves for HIV infection time data with us=0 and us=1, and average values of age and pyr. MABLE PO model is compared with semiparametric PH model (upper-left), semiparametric PO model (upper-right), and MABLE AFT model (lower-right).

Using the proposed method the ${\hat{p}}_{j}$ ’s that are greater than $10^{- 5}$ are $({\hat{p}}_{0}, {\hat{p}}_{1}, {\hat{p}}_{7})$ = $(0.05355, 0.03978, 0.90666)$ . The estimated baseline cumulative distribution function has a simple expression

{\hat{F}}_{0} (t) = 0.05355 B_{7, 0} (t / τ) + 0.03978 B_{7, 1} (t / τ) + 0.90666 B_{7, 7} (t / τ) .

Thus the survival function given

x = (us, age, pyr)

\hat{S} (t ∣ x) = \frac{1 - {\hat{F}}_{0} (t)}{1 + [\exp {0.5962 us - 0.0207 (age - 72) + 0.0122 pyr} - 1] {\hat{F}}_{0} (t)} .

6. Concluding remarks

As pointed out by a referee, the Bernstein polynomial model as a global approximation may be difficult to capture features such as sudden changes. The degree of approximation of the generalized Bernstein polynomial for any continuous function depends on its smoothness. Much larger model degree $m$ may be needed to achieve a good model fit for a survival function with sudden changes. The concerns about this issue and discontinuous baseline density warrants further generalization of the Bernstein polynomial model in the future communication.

Footnotes

Acknowledgments

The author is grateful for the three anonymous reviewers for their careful reviews and insightful and useful comments.

ORCID iD

Zhong Guan

Funding

The author received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interest

The author have no conflicting interests to declare with respect to the research, authorship, and/or publication of this article.

Supplemental material

Supplemental material for this article is available online.

References

Buckley

James

. Linear regression with censored data. Biometrika 1979; 66: 429–436.

Wei

. The accelerated failure time model: a useful alternative to the Cox regression model in survival analysis. Stat Med 1992; 11: 1871–1879.

Tian

Cai

. On the accelerated failure time model for current status and interval censored data. Biometrika 2006; 93: 329–342.

Guan

. Maximum approximate likelihood estimation in accelerated failure time model for interval-censored data. Stat Med 2023; 42: 4886–4896.

Cox

. Regression models and life-tables. J Roy Statist Soc Ser B 1972; 34: 187–220.

McCullagh

. Regression models for ordinal data. J R Stat Soc Ser B Methodol 1980; 42: 109–142.

Bennett

. Analysis of survival data by the proportional odds model. Stat Med 1983; 2: 273–277.

Bennett

. Log-logistic regression models for survival data. Appl Stat 1983; 32: 165–171.

Pettitt

. Proportional odds model for survival data and estimates using ranks. J R Stat Soc Ser C (Appl Stat) 1984; 33: 169–175.

10.

Doksum

. An extension of partial likelihood methods for proportional hazard models to general transformation models. Ann Stat 1987; 15: 325–345.

11.

Dabrowska

Doksum

. Partial likelihood in transformation models with censored data. Scand J Stat Theory Appl 1988; 15: 1–23.

12.

Scharfstein

Tsiatis

Gilbert

. Semiparametric efficient estimation in the generalized odds-rate class of regression models for right-censored time-to-event data. Lifetime Data Anal 1998; 4: 355–391.

13.

Finkelstein

Wolfe

. A semiparametric model for regression analysis of interval-censored failure time data. Biometrics 1985; 41: 933–945.

14.

Finkelstein

. A proportional hazards model for interval-censored failure time data. Biometrics 1986; 42: 845–854.

15.

Schick

. Consistency of the GMLE with mixed case interval-censored data. Scand J Stat 2000; 27: 45–55.

16.

Huang

Wellner

. Interval censored survival data: a review of recent progress. In: DY L and TR F (eds.) Proceedings of the first seattle symposium in biostatistics, Lecture Notes in Statistics, Vol. 123, pp.123–169. New York, NY: Springer.

17.

Kaplan

Meier

. Nonparametric estimation from incomplete observations. J Am Stat Assoc 1958; 53: 457–481.

18.

Turnbull

. The empirical distribution function with arbitrarily grouped, censored and truncated data. J R Stat Soc Ser B (Methodol) 1976; 38: 290–295.

19.

Owen

. Empirical likelihood ratio confidence intervals for a single functional. Biometrika 1988; 75: 237–249.

20.

Huang

Rossini

. Sieve estimation for the proportional-odds failure-time regression model with interval censoring. J Am Stat Assoc 1997; 92: 960–967.

21.

Murphy

Rossini

van der Vaart

. Maximum likelihood estimation in the proportional odds model. J Am Stat Assoc 1997; 92: 968–976.

22.

Shen

. Propotional odds regression and sieve maximum likelihood estimation. Biometrika 1998; 85: 165–177.

23.

Rabinowitz

Betensky

Tsiatis

. Using conditional logistic regression to fit proportional odds models to interval censored data. Biometrics 2000; 56: 511–518.

24.

Wang

Dunson

. Semiparametric bayes’ proportional odds models for current status data with underreporting. Biometrics 2011; 67: 1111–1118.

25.

Pordeli

. A proportional odds model for regression analysis of case I interval-censored data. Singapore: Springer Singapore, 2016.

26.

Guan

. Efficient and robust density estimation using Bernstein type polynomials. J Nonparametr Stat 2016; 28: 250–271.

27.

Guan

. Bernstein polynomial model for grouped continuous data. J Nonparametr Stat 2017; 29: 831–848.

28.

Guan

. Fast nonparametric maximum likelihood density deconvolution using Bernstein polynomials. Stat Sin 2021; 31: 891–908.

29.

Guan

. Maximum approximate Bernstein likelihood estimation in proportional hazard model for interval-censored data. Stat Med 2021; 40: 758–778.

30.

Rossini

Tsiatis

. A semiparametric proportional odds regression model for the analysis of current status data. J Am Stat Assoc 1996; 91: 713–721.

31.

Hastie

Tibshirani

. Generalized additive models, Monographs on Statistics and Applied Probability, Vol. 43. Chapman and Hall, Ltd., London, 1990.

32.

Anderson-Bergman

. icenReg: regression models for interval censored data in R. J Stat Softw 2017; 81: 1–23.

33.

Cattell

. The scree test for the number of factors. Multivariate Behav Res 1966; 1: 245–276.

34.

Wang

Guan

. Bernstein polynomial model for nonparametric multivariate density. Statistics 2019; 53: 321–338.

35.

Lorentz

. The degree of approximation by polynomials with positive coefficients. Mathematische Annalen 1963; 151: 239–251.

36.

Fay

Shaw

. Exact and asymptotic weighted logrank tests for interval censored data: the interval R package. J Stat Softw 2010; 36: 1–34.

37.

Melbye

Biggar

Ebbesen

, et al. Seroepidemiology of HTLV-III antibody in Danish homosexual men: prevalence, transmission, and disease outcome. Br Med J (Clin Res Ed) 1984; 289: 573–575.

38.

Carstensen

. Regression models for interval censored survival data: application to HIV infection in Danish homosexual men. Stat Med 1996; 15: 2177–2189.

39.

Lawless

Babineau

. Models for interval censoring and simulation-based inference for lifetime distributions. Biometrika 2006; 93(3): 671–686.

40.

Carstensen

Plummer

Laara

, et al. Epi: a package for statistical analysis in epidemiology, 2022. https://CRAN.R-project.org/package=Epi. R package version 2.47.