Tsukuba Economics Working Papers
No. 2019-003

Accelerated Failure Time Models with Log-concave Errors
byRuixuan Liu and Zhengfei Yu
November 2019

UNIVERSITY OF TSUKUBA
Department of Economics
1-1-1 Tennodai
Tsukuba, Ibaraki 305-8571
JAPAN

Manuscript submitted to The Econometrics Journal , pp. 1–20.

Accelerated Failure Time Models with Log-concave
ErrorsRuixuan Liu† and Zhengfei Yu‡†

Department of Economics, Emory University, 201 Dowman Drive, Atlanta, GA, 30322

E-mail: ruixuan.liu@emory.edu‡

Faculty of Humanities and Social Sciences, University of Tsukuba, 1-1-1, Tennodai, Tsukuba, Ibaraki, Japan.

E-mail: yu.zhengfei.gn@u.tsukuba.ac.jp
Summary
We study accelerated failure time (AFT) models in which the survivor functionof the additive error term is log-concave. The log-concavity assumption covers large
families of commonly-used distributions and also represents the aging or wear-out phenomenon of the baseline duration. For right-censored failure time data, we constructsemi-parametric maximum likelihood estimates of the finite dimensional parameterand establish the large sample properties. The shape restriction is incorporated via a
nonparametric maximum likelihood estimator (NPMLE) of the hazard function. Our
approach guarantees the uniqueness of a global solution for the estimating equations
and delivers semiparametric efficient estimates. Simulation studies and empirical applications demonstrate the usefulness of our method.
Keywords: Accelerate Failure Time Models, NPMLE, Weighted Rank Estimation,
Shape Restriction.

1. INTRODUCTION
The accelerated failure time (AFT) model provides an attractive alternative to the popular proportional hazards model (Cox, 1972) for analyzing censored duration/failure
time data. Let Y denote the logarithm of the duration T , C be the corresponding (logtransformed) censoring time, ∆ = 1(Y ≤ C), and V = min(Y, C). The model of interestisYi = Xi0 β0 + εi ,i = 1, · · · , n,(1.1)

where Xi stands for d-dimensional covariates and the independent error term εi has some
unknown distribution F . We denote exponential transform of ε by T0 , which represents
the baseline duration variable in the absence of the covariate effect. The AFT model
directly examines the effect of covariates on the length of survival, in contrast to the
proportional hazards model which focuses on the hazard rate. In many applications, it iseasier to visualize the concept that a treatment intervention or exposure to certain environment increases or decreases the length of survival itself by some amount, as comparedto the notion that the hazard rate is changed.1
In the presence of censoring and unknown error distribution, weighted log-rank estimators (Prentice, 1978) have been proposed to estimate the unknown regression coefficient
β0 . The large sample properties of parameter estimates are established by Tsiatis (1990),
1 This natural type of regression relationship led Sir David Cox himself (Reid, 1994) to acknowledge“the physical or substantive basis for...proportional hazards models...is one of its weakness...accelerated
life models are in many ways more appealing because of their quite direct physical interpretation.”

2

Liu and Yu

Ritov (1990), Lai and Ying (1991), and Ying (1993). Semi-parametric efficient estimators are first discussed in Lai and Ying (1991) and subsequently studied by Zeng and Lin(2007) and Ding and Nan (2011) based on kernel and sieve methods, respectively. Despite
the aforementioned theoretical advance, the practical applications of these estimators of
the AFT model remain limited. A fundamental issue is whether the estimating equation
has a global unique solution at the true parameter value, as commented by Kalbfleischand Prentice (2002) [p.219] “when viewed as a function of β0 , these test statistics are
step functions, and furthermore they are, in general, not monotone in β0 . This gives riseto the possibility of multiple solutions to estimating equations.” For some special andinefficient weighting schemes such as the Gehan’s weight (Fygenson and Ritov, 1994),
the resulting estimating equation is monotone and thus delivers the global uniquenessand consistency properties. Jin et al. (2003) develop a feasible iterative algorithm thatuses the rank estimator based on Gehan’s weights as the pilot estimate. Instead of narrowing down the weighting scheme, Ying (1993) explore another path that utilizes shape
restrictions on the error term to show the global properties of weighted rank estimators.
We study the AFT model under the restriction that the survivor function of the error
term ε is log-concave. Throughout this paper, a function g(x) is said to be log-concave
if the function log(g(x)) is strictly concave. It is well-known that the log-concavity ofthe survival function is equivalent to the increasingness of the corresponding hazard
function2 . The main advantage of such a shape restriction is that the increasing hazard
function guarantees the uniqueness of global solution to the population-level estimating
equations for various weighted rank estimators (Ying, 1993). To the best of our knowledge, this key insight has not yet been converted to an efficient estimation method thatformally takes into account of the log-concave restriction. For that purpose, we propose
semi-parametric maximum likelihood estimation methods that incoporate such a restriction for the AFT model. Specifically, for any given β, we estimate the hazard rate of(Y − X 0 β) using nonparametric maximum likelihood estimator (NPMLE) and enforce
the monotonicity restriction as in Huang and Wellner (1995). Note that the general NPMLE for the hazard function does not exist without such shape restrictions; see Zeng and
Lin (2007). Thereafter, we obtain the estimate of β that solves the efficient score function(see Theorem 4.3 in Ding and Nan, 2011) by plugging in the NPMLE of the hazard rate
and the kernel smoothed estimate for the derivative of the hazard rate. We establish
the asymptotic distribution theory of the resulting estimators for the finite dimensionalparameter β0 and show that our estimators achieve the semi-parametric efficiency bound(Bickel et al., 1993). Note that the efficiency bound under our log-concave restriction
remains unchanged as the one under the standard smoothness assumption. We establish
this result by showing the tangent set under the shape restriction can be well approximated by the one under the smoothness restriction and vice versa. A formal proof can befound in Section S1.3 of the online supplement. In sum, the main appealing property of
our approach is that it achieves the same semiparametric efficiency as the kernel or sieve
based estimators (Lai and Ying, 1991; Ding and Nan, 2011) and at the same time guarantees the uniqueness of a global solution for the estimation equations (and the resulting
computational convenience), which previously was only shown for the inefficient Gehan’sweighting scheme. In principle, a shape-constrained sieve MLE can serve as an alterna2 See Proposition C.1.c. in Chapter 4 of Marshall and Olkin (2007) and the proof therein. Throughout this paper, the increasingness (decreasingness) of a function refers to the strictly increasingness(decreasingness).

Logconcave-error-AFT

3tive for incorporating the log-concave restriction. However, it usually requires imposing
additional inequality restrictions on the coefficients of basis functions, which may further
complicate the numerical optimization procedure; see Section 2.3.5 of Chen (2007). In
contrast, our NPMLE-based approach is made possible by the underlying log-concave
restriction and automatically satisfies it. In this sense, the shape restriction is a blessing
instead of a burden for our estimation approach.
The framework of log-concave errors comprises numerous parametric duration models
where the failure time T (conditional on covariates) can follow exponential, Weibull,
Gamma, Log-normal, Log-logistic, and many other distributions (Bagnoli and Bergstrom,
2005; Marshall and Olkin, 2007). It is worth clarifying the relationship between the logconcavity restriction on ε and the shape of the baseline hazard function of T0 = exp(ε).
On one hand, if one believes that baseline duration T0 exhibits some aging or wear-out
phenomenon, as encoded by an increasing baseline hazard function, then its logarithm
transform ε will also have an increasing hazard rate. This feature often arises from a
non-stationary job search model where the increasing exit rate of unemployment status
is caused by a declining reservation wage and/or an increasing search intensity (Burdett
and Vishwanath, 1988; Van den Berg, 1990). On the other hand, the increasing hazardof ε itself does not necessarily restrict the shape of the hazard of T0 . Our Example 2.1
provides an encompassing parametric family in which the hazard shape of T0 can beincreasing, decreasing, or non-monotonic, under the log-concavity restriction on ε.
In establishing the large sample properties of our estimators, our main technical contribution is a formal analysis of the asymptotics of the NPMLE of the hazard rate andits kernel smoothed derivative, adapting the recent development from Groeneboom et al.
(2010), Groeneboom and Hendrickx (2018), and Westling and Carone (2019). This is
not a standard problem, as the NPMLE of the hazard rate is a piece-wise constantfunction with random jump locations determined by the data, and the kernel smoothed
estimate of NPMLE is not a linear functional of the empirical measure (Groeneboom
et al., 2010). Compared with the binary choice model in Groeneboom and Hendrickx(2018), the characterization of NPMLE is more complicated in our setting, because its
min-max representation in equation (3.15) involves random denominators, whereas theone in Groeneboom and Hendrickx (2018) does not; see Westling and Carone (2019)
for more discussions. Our proofs that combine empirical process theory and the characterization of the (smoothed) NPMLE for the monotone hazard are also of independentinterest. The shape restricted estimation and inference constitute a rich and evolving literature in econometrics and statistics, as reviewed by Groeneboom and Jongbloed (2014)
and Chetverikov et al. (2018). The rest of the paper is organized as follows: Section 2
discusses the log-concavity restriction. Section 3 proposes efficient semiparametric maximum likelihood estimation (SMLE) methods subject to the log-concave restriction on
the error term. Section 4 derives the asymptotic properties of the parameter estimates,
proving the consistency, asymptotic normality, and semi-parametric efficiency. Section 5
conducts simulation studies and applies the proposed methods to a real data set. The
final section concludes. Main results are proved in Appendix. Other proofs, technicallemmas, and an additional real data example are relegated to the online supplement.
2. LOG-CONCAVITY OF ERROR TERMS IN AFT MODELS
We begin with a clarification about the log-concavity restriction used in this paper.
Proposition 2.1 collects some important properties related to the log-concavity restric-

4

Liu and Yution and states the relationship between the hazard functions of ε and T0 . Let f (u)be the (absolutely continuous) density function of the error term ε in (1.1), F (u) be
the cumulative distribution function, F̄ (u) = 1 − F (u) be the survivor function andλ(u) = f (u)/F̄ (u) be the hazard function. Let h(u) be the hazard function of the baseline duration T0 = exp(ε). The log-concave restriction in this paper refers to the survivor
function F̄ (u) being log-concave.
Proposition 2.1. (i). If f (u) is log-concave, then both F (u) and F̄ (u) are log-concave.
(ii). F̄ (u) is log-concave if and only if λ(u) is increasing. (iii). If h(u) is increasing, thenλ(u) is also increasing.
As shown in part (i) of Proposition 2.1, the log-concavity restriction on F̄ can be derived
from the log-concave restriction on the density function f . For parametric models, it isoften easier to directly verify the log-concavity of the density function. Example 2.1 below
shows that a large family of commonly-used parametric duration models falls into ourmodeling framework. From the perspective of semiparametric estimation for the AFTmodel, the log-concave restriction yields nice properties for weighted rank estimation.
Specifically, Lemma 4.1 shows that an increasing hazard λ is sufficient for the populationlevel estimating equation to have a unique global solution.
Example 2.1. (Generalized F -distribution). Various frequently used continuous parametric failure time models are special cases of the model (1.1) in which the error term
ε is distributed as the logarithm of the F -variate with the degrees of freedom equal to(2m1 , 2m2 ) on [p.38] of Kalbfleisch and Prentice (2002). The density function of ε writes:m1

f (u) = (m1 /m2 )em1 u (1 + (m1 /m2 )eu )−(m1 +m2 )/B(m1 , m2 ),(2.2)

where B(m1 , m2 ) is Beta function. The density function (2.2) reduces to the logistic
density when (m1 , m2 ) = (1, 1); it reduces to the extreme value density (that leads to the
Weibull model) when (m1 , m2 ) = (1, ∞); it generates the Gamma model when m2 = ∞;it approaches to the standard normal density as (m1 , m2 ) → (∞, ∞). It is easy to showthat the second order derivative of the logarithm of f (u) is200(log f (u)) = −(m1 + m2 )(m2 /m1 )e−u / 1 + (m2 /m1 )e−u < 0,
(2.3)
for any m1 , m2 > 0. Hence, the density function (2.2) is log-concave and the model isnested by our setup.
Note that the converse of part (iii) of Proposition 2.1 is not true. In other words, ourlog-concave restriction on ε does not restrict the hazard function of T0 to be monotone.
One can easily see this by inspecting the following expression:h(t) = [f (log(t))/t] /F̄ (log(t)) = λ(log t)/t.

(2.4)

Therefore, our shape restriction can still be plausible in the cases where empirical evidence
rejects a monotonic hazard function of T0 ; see Christofides and McKenna (1996) for suchan example3 . Figure 1 plots hazard functions of the error term ε and the baseline duration
T0 for parametric models encompassed by Example 2.1. Despite an increasing hazard rate
3 Christofides and McKenna (1996) model the baseline duration by a generalized Gamma distributionthat is log-concave, in order to accommodate a non-monotonic hazard of the unemployment duration.

Logconcave-error-AFT

5

Figure 1. Hazard functions of the error term ε (left) and the baseline duration T0 (right)

for various parametric models.

for ε, the hazard rate of its exponential transform T0 can be increasing, decreasing or
non-monotonic. We observe the following: (i). The log-normal failure time model, in
which ε has a standard normal distribution, yields a increasing hazard function for ε,but an inverse U-shaped hazard function for T0 . (ii). The Weibull model, in which thebaseline hazard function takes the form h(t) = γtγ−1 , has increasing hazard functions
for both ε and T0 when the shape parameter γ > 1; it implies a increasing hazard for εand a decreasing hazard for T0 when 0 < γ < 1. (iii). The Log-logsitic model, in which
ε has a logistic density eu/σ /σ(1 + eu/σ )2 , produces a increasing hazard function for εand an inverse U-shaped function for T0 , when 0 < σ < 1 (σ = 0.7 in the figure). (iv).
The Gamma model, in which T0 has a density tk−1 e−t /Γ(k), yields increasing hazard
functions for both ε and T0 when the shape parameter k > 1 (k = 3 in the figure).
Inspired by Abbring (2012) which studies durations defined by threshold-crossing rules,
we provide a model that yields an increasing hazard function for T0 . Conditional on
X = x, the duration T is defined as the survival time after a random number of shocks
and the arrival of shocks is governed by a Poisson process Nx (t), so that
Pr {T > t|X = x} =∞Xk

(λx t) e−λx t /k!P k ,

(2.5)

k=0

where λx ≡ exp(−x0 β) is the failure rate of Nx (t) and P k is the probability of surviving
after k shocks. The following two examples are from Esary and Marshall (1973) and spellout the interpretative content of the underlying failure mechanism based on different
specifications of the sequence P k .
Example 2.2. (Cumulative Shock Models) The failure is caused by cumulative shocks
U1 , U1 + U2 , · · · reaching some positive random threshold Z, as in Section 5 of Esary
and Marshall (1973), with the following specification on P̄k for k ≥ 1:
Z ∞(k)
P k = Pr {U1 + ... + Uk ≤ Z} =FU (z)dGZ (z),0

6

Liu and Yu(k)where GZ is the cumulative distribution function of Z, and FU stands for the k-th foldconvolution of the function FU , which is the distribution function of i.i.d. shocks Uk fork ≥ 1.
Example 2.3. (Maximum Shock Models) The failure is caused by the maximum of independent shocks U1 , U2 , · · · reaching a fixed threshold z, as in Section 6 of Esary andMarshall (1973), with the following specification on P̄k for k ≥ 1:
P k = Pr {max(U1 , · · · , Uk ) ≤ z} = Πki=1 FUi (z),where the distribution function FUk is not necessarily identical for k ≥ 1.
We show in the following proposition that the construction via equation (2.5) generates
an AFT model with some baseline duration T0 that has an increasing hazard function
under reasonable assumptions on P k for both Examples 2.2 and 2.3. Therefore, thecorresponding additive error term ε also has an increasing hazard rate, according to part(iii) of Proposition 2.1.
Proposition 2.2. The duration T defined by equation (2.5) generates an AFT modelwith a log-concave baseline duration T0 if the following conditions hold: (i) For the cumulative shock models, both the shock U and threshold Z have log-concave densities; (ii)
For the maximum shock models, FUk (u) is decreasing in k for any u (so that subsequent
shocks are stronger).
3. SEMI-PARAMETRIC EFFICIENT ESTIMATION WITH SHAPERESTRICTIONS
The data consists of i.i.d. observations of {Vi , Xi , ∆i ; i = 1, · · · , n}, where Vi = min(Yi , Ci )and ∆i = 1(Yi ≤ Ci ). The error term εi in (1.1) is assumed to be independent of Xi andthe right-censoring time Ci . Let the marginal distribution function of covariates X be
H(x) for x in the support X . Following Ying (1993), we denote the conditional distribution and density functions of the censoring variable C as Gx (·) and gx (·), respectively.
For a given β, define β = V − X 0 β and 0 = V − X 0 β0 . We use the standard empirical
process notations as follows. For a function f of a random vector Z = (V, X, ∆) that
follows distribution P ,ZnX
P f = f (z)dP (z), Pn f = n−1
f (Zi ), Gn f = n1/2 (Pn − P ) f.
(3.6)
i=1

Function f can be replaced by a random function z 7→ fˆn (z; Z1 , · · · , Zn ) (Nan andWellner, 2013), so we also write P fˆn , Pn fˆn , or Gn fˆn . Furthermore, defineDn(0) (t, β) = Pn {1(β ≥ t)},Dn(1) (t, β) = Pn {X1(β ≥ t)},D(0) (t, β) = P {1(β ≥ t)},D

(1)(t, β) = P {X1(β ≥ t)},

(3.7)
(3.8)

and
η̂n (t, β) = Dn(1) (t, β)/Dn(0) (t, β).

(3.9)

The probability limit of η̂n (t, β) for any given β is
η0 (t, β) = D(1) (t, β)/D(0) (t, β) = E[X|V − X 0 β ≥ t].

(3.10)

Logconcave-error-AFT

7

In addition, we defineNn (t, β) = Pn {1(β ≤ t, ∆ = 1)},N (t, β) = P {1(β ≤ t, ∆ = 1)}.

(3.11)The likelihood function (omitting terms irrelevant for our parameters of interest) is
Z Vi
n n
o
X
−1
0
Ln (β, λ) = n∆i log λ(Vi − Xi β) −λ(t − Xi0 β)dt ;(3.12)−∞

i=1

also see Zeng and Lin (2007) and Ding and Nan (2011). The unknown parameters consistof β ∈ B and λ(t) ∈ A, where B is a compact set in Rd and A is the space for increasing
and non-negative functions. We first derive λ̂n (·, β) that maximizes (3.12) for any fixed
β. It is an estimator of the increasing function λ0 (t, β) = (dN (t, β)/dt) /D(0) (t, β) andcan be obtained by applying NPMLE for the monotone hazard function (Huang and
Wellner, 1995). Recall that i,β = Vi − Xi0 β, and the corresponding indicator is ∆i,β fori = 1, · · · , n. Let (1),β ≤ · · · ≤ (n),β be the order statistics, and ∆(i),β be the indicatorassociated with the i-th order statistic (i),β . Define the weights
Zwi,β = n{(i+1),β − (i),β }dPn (u), for i = 1, ..., n − 1.
(3.13)u≥(i+1),βAssuming that the function λ0 (t, β) ≤ Mλ for some constant Mλ and for all β ∈ B,4
the NPMLE λ̂n (·, β) can be characterized as the left derivative of the greatest convex
minorant of the following “cumulative sum diagram” with points


i
i
X
X
P0 = (0, 0), and Pi = 
wj,β ,
wj,β ∆(i),β /wi,β  , for i = 1, · · · , n − 1. (3.14)
j=1

j=1

By the min-max characterization (Groeneboom and Jongbloed, 2014), the NPMLE is a
step-wise constant function uniquely determined at i,β for i = 1, · · · , n − 1 asPt
j=s ∆(j),β,(3.15)λ̂n ((i),β , β) = max min Pt
1≤s≤i i≤t≤nj=s (n − j + 1)((j),β − (j−1),β )and λ̂n ((n),β , β) = Mλ . In practice, it is enough to compute λ̂n ((i),β , β) at the firstn − 1 points as in (3.15) and thus there is no need to specify Mλ . The established pooladjacent-violators algorithm (PAVA) can be used to compute NPMLE in a very efficientway; see Groeneboom and Jongbloed (2014). The probability limit of λ̂n (·, β) is λ0 (t, β).
We also denote λ0 (t) ≡ λ0 (t, β0 ). It is worthwhile emphasizing that for any given β, thelikelihood function cannot be meaningfully maximized without the shape restriction on
λ (Zeng and Lin, 2007).
Referring to the semi-parametric efficient estimator presented below, smoothing isinevitable, because the efficient score function involves the derivative of the hazard function. Thus, we adopt the smoothed maximum likelihood estimator in Groeneboom et al.
4 This boundedness assumption is required for maximizing the likelihood in (3.12). Otherwise,
λ̂n ((n),β , β) can be chosen arbitrarily large. After obtaining the solution under the upper bound Mλ ,one can allow Mλ to grow arbitrarily large. This type of argument is common in isotonic optimization
problems, see [p338-339] of Robertson et al. (1988) and [p38-39] of Groeneboom and Jongbloed (2014).
In computation, one can set λ̂n ((n),β , β) = λ̂n ((n−1),β , β).

8Liu and Yu(2010) and obtain˙
λ̂n (t, β) =Z
Kh (t − u)dλ̂n (u, β),(3.16)

for a kernel density function K(·) and Kh (·) ≡ K(·/h)/h, with the bandwidth equal toh. It is straightforward to find that the probability limit of the kernel smoothed estimateisλ̇0 (t, β) = ∂λ0 (t, β)/∂t.

(3.17)

To estimate β0 , we consider the following random mapnΨn (β, η, ρ) ≡1X
∆i ρ(Vi − Xi0 β, β){Xi − η(Vi − Xi0 β, β)},n i=1(3.18)

where β ∈ B is the d-dimensional Euclidean parameter of interest with an unknown true
value of β0 , η and ρ are functions that can be viewed as infinite dimensional nuisanceparameters. Intuitively speaking, η provides the correct centering term for the regressor
X under the censoring, whereas the weight ρ influences the efficiency. In this paper,
we consider Ψn (β, η̂n , ρ̂n ) which approximates the efficient score function and belongs
to the general weighted rank estimating equation as in Tsiatis (1990), with an efficientweighting scheme. Specifically, these nuisance components are estimated by˙η̂n (t, β) = Dn(1) (t, β)/Dn(0) (t, β), ρ̂n (t, β) = I{t ≥ τU }λ̂n (t, β)/λ̂n (t, β).

(3.19)

Following Tsiatis (1990) (see also Lai and Ying, 1991; Ding and Nan, 2011), the trimming(0)
constant τU is introduced to avoid the instability of Dn and λ̂n near the right tail. The
corresponding probability limit of ρ̂n (t, β) is denoted by ρ0 (t, β).
Our estimator β̂n is the root of the estimating function Ψn (β, η̂n (·, β), ρ̂n (·, β)), which
is discontinuous in β. Therefore, Ψn (β, η̂n (·, β), ρ̂n (·, β)) = 0 may not hold exactly. Following Groeneboom and Hendrickx (2018), we define β̂n as a zero-crossing point suchthat each coordinate value of Ψn (β, η̂n (·, β), ρ̂n (·, β)) changes sign in the left and rightneighborhoods of β̂n .
Definition 3.1. (Zero-crossing) We say that β∗ is a zero-crossing of a function
C : B 7→ R if each open neighborhood of β∗ contains points β1 , β2 ∈ B such thatC(β1 )C(β2 ) ≤ 0. Moreover, we say that a function C̃ : B 7→ Rd crosses zero at point β∗
if β∗ is zero-crossing in each component C̃j for j = 1, · · · , d.
We propose another estimation procedure which uses a smoothed version of λ̂n (t, β).
It is well documented in the literature (Groeneboom et al., 2010) that smoothed NPMLE
often performs better in finite samples. For the kernel density function K(·), we defineRt
K(t) = −∞ K(u)du and Kh̃ (u) = K(u/h̃). The smoothed estimator for λ0 (t, β) isZ
λ̃n (t, β) = Kh̃ (t − u)dλ̂n (u, β).
(3.20)
Therefore, the alternative estimating equation isnΨ̃n (β, η, ρ) ≡1X∆i ρ̃n (Vi − Xi0 β, β){Xi − η̂n (Vi − Xi0 β, β)},n i=1(3.21)

Logconcave-error-AFT

9

with
˙
ρ̃n (t, β) = I{t ≥ τU }λ̂n (t, β)/λ̃n (t, β).

(3.22)

The corresponding solution is denoted by β̃n .
Figure 2 illustrates typical shapes of estimating functions Ψn and Ψ̃n for a single
parameter β using the unsmoothed and smoothed estimates of λ0 (t, β). The sample size
is 500 and is generated from a N (2, 0.52 ) covariate, a standard normal error term, anda uniform censoring variable on [0, 4]. The true parameter is β0 = 1. The bandwidths
˙
are 5 × n−1/5 for λ̃n (t, β) and 5 × n−1/7 for λ̂n (t, β), respectively. Both methods deliver
unique zero-crossing points, whereas Ψ̃n exhibits smoother variation.

Figure 2. The estimating equations Ψn (left) and Ψ̃n (right) as functions of the single
parameter β.

Remark 3.1. Lai and Ying (1991) are among the first who propose to use the kernelsmoothing methods to estimate the underlying hazard and its derivative functions in theefficient score function. Their specific construction consists of (i) dividing the sample into
two disjoint subsets and evaluating a preliminary consistent estimate b̂j of β from thejth subsample (j = 1, 2), (ii) finding from the uncensored residuals in the jth subsample
a smooth consistent estimate λ̂j of the hazard function, (iii) smoothing the estimatedhazard to obtain a smooth consistent estimate ρ̂j of λ̇/λ, and (iv) using ρ̂1 (respectivelyρ̂2 ) as the weight function for the linear rank statistic of the second (respectively first)
sample of residuals Yi − Xi0 b. The sum S(b) of these two linear rank statistics is used
to define the rank estimator as the minimizer of k S(b) k. However, there are practicaldifficulties in carrying out this procedure, as reviewed by Kim and Lai (2000). First, thisefficient rank estimator is difficult to compute when the regressor X is multidimensional.

10

Liu and Yu

Second, the kernel smoothing approach employed does not give good results unless the
sample size is very large.
Remark 3.2. Besides the guarantee of the unique solution, another advantage of our efficient weighted rank estimator over the sieve MLE in Ding and Nan (2011) stems fromits computational convenience. In the former, only the regression parameter β needs tobe solved from the estimating equations, whereas both β and the spline coefficients (thenumber of such coefficients also diverges with the sample size to eliminate the approximation bias) need to be solved in the latter. For example, the sieve approach using cubic
splines and two internal knots has to solve six more parameters from the optimization
procedure than our approach. Meanwhile, the pool-adjacent-violators algorithm (PAVA)
used to compute the NPMLE of the hazard rate (Groeneboom and Jongbloed, 2014) isfaster than the iterative Newton-Raphson method for the sieve MLE estimate.
4. ASYMPTOTIC PROPERTIESThe following regularity conditions, adapted from Tsiatis (1990), Lai and Ying (1991),
and Ying (1993), are imposed throughout. We let T to be the support of ε trimmed fromthe right tail at some given τU (Tsiatis, 1990).
Assumption 4.1. The covariates Xi are uniformly bounded, i.e., max1≤i≤n kXi k ≤ Mx
for some finite constant Mx .
Assumption 4.2. The error term ε is independent of (X, C) and has a finite mean. The2
R∞ density f and its derivative f˙ are bounded. Moreover,f˙(t)/f (t) f (t)dt < ∞.
−∞Assumption 4.3. The conditional density gx (t) of the censoring time C is uniformly
bounded; that is, supt∈C,x∈X |gx (t)| ≤ Mc for some finite constant Mc , where X is the
support of covariates X and C is the support of censoring variable C.
Assumption 4.4. The hazard rate λ0 of error term ε is increasing and is continuously
third-order differentiable. Moreover, λ0 (t, β) is uniformly bounded away from zero forany β in the parameter space and t ∈ T . Its derivative λ̇0 (t, β) is also uniformly bounded
away from 0 and ∞.
Assumption 4.5. For the finite positive constant τU , there exists a value of ξ such that
P r(V − X 0 β0 ≥ τU ) ≥ ξ > 0.
Assumption 4.6. R(i) The kernel function
K(u) is a kernelRdensity function with compact
Rsupport such that K(u)du = 1, uK(u)du = 0, and u2 K(u)du < ∞. Moreover,its derivative k(u) is uniformly continuous over the support. The kernel function andits derivative can be written in the form of ϕ (p (x)), with some function ϕ (·) being ofbounded variation and p (x) a real polynomial on R. (ii) The bandwidths satisfy h  n−1/7
and h̃  n−1/5 .
Assumptions 4.1 to 4.3 are from Ying (1993). The monotonicity restriction in Assumption 4.4 results from our log-concave restriction on the survival function of ε throughProposition 2.1 (ii). The class of hazard rate functions that are bounded away from zero

Logconcave-error-AFT

11

when t approaches to the left boundary is discussed extensively in Ridder and Woutersen(2003). With log-concave errors, the moment bound in Condition 4 of Ying (1993) is automatically satisfied. When the error term ε has a finite mean and an increasing hazard
rate, it has a sub-exponential tail, cf. Theorem 4.1 in Barlow and Marshall (1964), which
gives rise to finite moments of all orders. In accordance with Assumption 4.5, we trim
the right tail at some large constant τU ; see equation (3.1) of Tsiatis (1990). This is astandard practice in analyzing censored regression models (Lai and Ying, 1991; Ritov,
1990; Ding and Nan, 2011). In fact, the calculation of the semiparametric information
bound takes this trimming parameter τU as a fixed constant; see Example 4 on [p.284] ofBickel et al. (1993). It is possible to allow τU to be infinity through more technical proofs(see Lemma 2 of Ying, 1993), but it will not be attempted here. In Assumption 4.6, wecollect the standard requirements on kernel smoothing methods, as in Groeneboom et al.
(2010). Under the restriction, we have the following two VC-type functional classes due
to Nolan and Pollard (1987):
K1 = K h−1 (x − ·) : x ∈ R, h > 0
and K2 = k h−1 (x − ·) : x ∈ R, h > 0 ,ˆwhich are needed to establish the convergence rates for λ̇n and λ̃n . The population level
estimating equations can be written asΨ(β, η0 (·, β), ρ0 (·, β)) = P [∆ρ0 (V − X 0 β, β){X − η0 (V − X 0 β, β)}] .

(4.23)To show the existence of a solution of estimating equations (3.18) and (3.21) with probability tending to one and its consistency, we verify that the population level estimating
function Ψ(β, η0 (·, β), ρ0 (·, β)) has a unique root under the shape restriction given by
Assumption 4.4. This result essentially follows from the discussion in Section 5 of Ying
(1993) and we restate it as Lemma 4.1 for completeness.
Lemma 4.1. If the hazard function λ of the error term ε is increasing, as assumed inAssumption 4.4, then β0 ∈ B is the unique solution to the population level estimating
functions Ψ(β, η0 (·, β), ρ0 (·, β)).
In particular, Ying (1993) shows that the population criterion function has a uniquezero-crossing point globally when the error term has an increasing hazard function andthe weights do not have alternating signs (either all positive or all negative). In our case,the weighting function ρ̂n (·) defined in (3.19) is automatically non-negative, so is its
probability limit. In sharp contrast, standard kernel or sieve type estimators for ρ̂n (·)
without such shape constraint fail to deliver a unique global solution; see Kim and Lai(2000). Based on Lemma 4.1, the consistency of our estimators can be obtained from theuniform convergence of Ψn (β, η̂n (·, β), ρ̂n (·, β)) and the Glivenko-Cantelli theorem.
Theorem 4.1. (Consistency) Supposethat Assumptions 4.1 to 4.6 hold. Then, for all
large n, a zero-crossing β̂n for Ψn β̂n , η̂n (·, β̂n ), ρ̂n (·, β̂n ) exists with probability tending
to one and is a consistent estimatorof β0 . The same conclusion holds for β̃n , which is azero-crossing point for Ψ̃n β̃n , η̂n (·, β̃n ), ρ̃n (·, β̃n ) .
Because our estimation procedure belongs to the general Z-estimation with bundledparameter and nonparametric nuisance components, we prove the root-n rate and asymptotic normality of β̂n following the route in Nan and Wellner (2013). Note that under our12

Liu and Yu

assumptions, ρ0 (β , β) and η0 (β , β) are both continuously differentiable with derivatives
denoted by ρ̇0β and ρ̇0β , respectively. As shown in Theorem 3.3 in Nan et al. (2009),
this implies that Ψ(β, η0 (·, β), ρ0 (·, β)) is differentiable in β with the bounded derivative Ψ̇β (β, η0 (·, β), ρ0 (·, β)) in B. Given the explicit form of Ψ(β, η(·, β), ρ(·, β)), the pathderivatives (Ichimura and Lee, 2010) with respect to both η and ρ exist, and we denotethem by Ψ̇η (β, η(·, β), ρ(·, β)) and Ψ̇ρ (β, η(·, β), ρ(·, β)), respectively.
Theorem 4.2. (Asymptotic Normality) Suppose that Assumptions 4.1 to 4.6 hold.
Then, we have the following linear representationZ
1/2−Ψ̇β (β0 , η0 (·, β0 ), ρ0 (·, β0 ))n (β̂n − β0 ) = Gn ρ0 (t, β0 ){X − η0 (t, β0 )}dM (t) + op (1),where M (t) is the martingale:
M (t) = ∆1(V − X 0 β0 ≤ t) −Z t

1(V − X 0 β0 ≥ s)λ0 (s)ds.

(4.24)−∞

Therefore, we obtainn1/2 (β̂n − β0 ) → N(0, I −1 (β0 ))in distribution. The information matrix I(β0 ) is the semi-parametric efficient information
matrix given by Lemma S1 in the online supplement. The same conclusion holds for β̃n .
The calculation of semi-parametric information bound embedding shape restrictions is
not trivial; see Tripathi (2000), and Kuchibhotla et al. (2017). We formally verify in the
online supplement that for the increasing hazard rate, the information bound remains
unchanged as in Bickel et al. (1993). In order to determine the bound, the tangent set andthe projection of the (parametric) score function to the tangent set need to be calculated.
The score function is not affected by the shape restriction. The crux in our proof is toshow the tangent set remains unchanged by showing that scores for smooth sub-models
lie in the set and by exhibiting a family of smooth sub-models with scores that canapproximate any element of the set arbitrarily well are dense in the set.
Remark 4.1. The standard errors of β̂n can be obtained based on the following result:
n
i⊗2
1X h
∆i ρ̂n (Vi − Xi0 β̂n , β̂n ){Xi − η̂n (Vi − Xi0 β̂n , β̂n )}→p I(β0 ).
(4.25)n i=1
In the proofs of Theorem 4.1 and Theorem 4.2, we have Ψn β̂n , η̂n (·, β̂n ), ρ̂n (·, β̂n ) =Pn {lβ∗0 (V, ∆, X)} + op (1), where lβ∗0 is the efficient score function such that I(β0 ) =E[lβ∗0 (Vi , Xi , ∆i )⊗2 ]. The experssion of lβ∗0 is given by Lemma S1 in the online supplement. An adaption to the result in equation (4.25) by the Glivenko-Cantelli property isstraightforward; see Theorem 4.3 in Ding and Nan (2011).
Remark 4.2. Our asymptotic analysis also produces a uniform consistent estimator forthe hazard function λ0 with a cubic-root rate (modulo some logarithm term):k λ̂n (t, β̂n ) − λ0 (t) k2 = Op (n−1/3 log n).
Given the natural bound via the triangular inequality:k λ̂n (t, β̂n ) − λ0 (t) k2 ≤k λ̂n (t, β̂n ) − λ0 (t, β̂n ) k2 + k λ0 (t, β̂n ) − λ0 (t) k2 ,(4.26)

Logconcave-error-AFT

13(4.26) follows from Lemma S8 in the online supplement and the root-n consistency ofβ̂n .
Remark 4.3. The foundation to incorporate the time-dependent covariates in AFT models is laid down by Robins and Tsiatis (1992) and further developed by Lin and Ying
(1995) and Zeng and Lin (2007). In this setup, it is more convenient to work with theduration T without making the logarithmic transformation. Recall in the setting with time
independent covariates X, the essence of the AFT model is that there exists some inde0
pendent baseline duration T0 such that T0 = T e−X β0 . With time-dependent covariates
X(t), the following specification generalizes the AFT model:
Z Ti
exp(Xi0 (s)β0 )ds, i = 1, · · · , n,(4.27)
T0,i =
0where T0,i is independent of the covariates. Lin and Ying (1995) have derived the efficientestimating function extending equation (3.18), whereas Zeng and Lin (2007) formallydevelop a semiparametric efficient estimator for the AFT model with time-dependent
covariates. It is interesting to extend our methodology to this scenario in future work.
5. NUMERICAL RESULTS
5.1. Monte Carlo Simulations
We conduct Monte Carlo simulations to evaluate the finite sample performances of ourestimators. The logarithm of the duration Y = log T is generated from the model:
Y = 2 + β1 X1 + β2 X2 + ε, where β1 = β2 = 1. The covariate X1 is Bernoulli witha success probability 0.5, and the covariate X2 is normal with mean zero and standard
deviation 0.5. X1 and X2 are independent. We consider two types of log-concave errors:(I). The error term ε is the logarithm of a Weibull-distributed variable. The density of
ε is f (u) = γeγu exp (−eγu ), where γ = 1.2. In this case, the duration T conditional on (X1 , X2 ) has a Weibull distribution and its hazard function is increasing. (II).
The error term ε is the logarithm of a Gamma-distributed variable. The density of ε isf (u) = exp (ku − eu ) /Γ(k), where k = 0.7. In this case, the duration T conditional on(X1 , X2 ) has a Gamma distribution and its hazard function is decreasing. Recall thatV = min(Y, C) and ∆ = 1{Y ≤ C}. The censoring time C is drawn from a uniform distribution supported on [0, τ ], where τ is chosen to produce a 25% censoring rate.
The data in each replication is an i.i.d. sample of {Vi , ∆i , X1i , X2i }ni=1 for n = 500 and
1000. We obtain semi-parametric maximum likelihood estimates of (β1 , β2 )0 using an unsmoothed hazard rate (SMLE-u) and a smoothed hazard rate (SMLE-s), by repetitively
solving the estimating equations (3.18) and (3.21). We use the quartic kernel function
2
15
K(u) = 16
1 − u2 1{|u| ≤ 1}, which satisfies Assumption 4.6. The bandwidths for smoothing the estimated hazard rate and its derivative are set to hc ×n−1/7 or hc ×n−1/7 ,respectively. The constant hc = 2 and 5. Regarding the right-tail trimming factor τU in
Assumption 4.5, we trim at the 98% quantile of Y − β1 X1 − β2 X2 for each (β1 , β2 ).
As comparison, we also include the Gehan-weight estimator (Jin et al., 2003), which isinefficient, and the sieve ML estimator (Ding and Nan, 2011, using cubic splines and twointerior knots). In our simulations, SMLE methods are computationally much faster than
the sieve estimator. On a 3.0 GHz Intel Core i5 processor and with 10 GB of RAM, SMLE
methods on average take about 25 seconds (elapsed time) while the sieve estimator takes

14

Liu and Yu

about 300 seconds for each replication with sample size equal to 500. This is mainly because SMLE methods solve two-dimensional parameter β from the estimating equations,whereas the sieve estimator has to solve both β and additional six spline coefficients.
Table 1 presents the bias, empirical standard error (SE) and square root of the meansquared error (RMSE) for SMLE-u, SMLE-s, Genhan and sieve estimators. Results are
based on 1000 simulated samples. We make the following observations when comparing these estimators. First, the performances of the proposed SMLE methods are stable
across two sets of bandwidths. SMLE-s using the smoothed hazard rate yields smaller
empirical standard errors than its unsmoothed version SMLE-u in all scenarios. Second,
SMLE methods, especially SMLE-s, turn out to be more efficient (having smaller empirical standard error) than Gehan-weighted estimator in all scenarios. This substantiates
the efficiency gain of SMLE methods over the rank estimator using inefficient weights.
Third, the standard errors of SMLE-s are similar to those of the efficient sieve estimatorin all scenarios. Fourth, the SMLE methods and the Gehan-weighted estimator yieldsmaller bias for β1 than the sieve estimator. Bias performances are comparable for β2 .
In sum, our SMLE methods achieve the same accuracy as the sieve ML estimator andenjoy two advantages: the guarantee of a unique global solution and faster computation.
5.2. Real Data Example

We re-visit the 4-month panel of the Current Population Survey (CPS) from September
to December 1993, previously studied by Romeo (1999). We analyze the unemployment
duration using the AFT model and focus on the sample with a positive employment
duration of less than a year. The sample size is 399 and the censoring rate is 42%.
We regress the natural logarithm of the unemployment duration (in weeks) on threecovariates: race, gender, and a re-entrant indicator that indicates whether the individual
is a labor force re-entrant. Table 2 reports coefficient estimates and standard errors. All
methods suggest that on average, white, female, and re-entrants to the labor market haveshorter unemployment durations. The negative sign of the coefficient on the re-entrant
indicator suggests that people not in the labor force are more likely to return when they
predict that there is a good chance of finding a job (Romeo, 1999). The magnitudes of
estimates do not vary much across different methods. The proposed SMLE methods givesmaller standard errors than other methods.
Figure 3 plots the estimated hazard of error term ε and the baseline hazard for theunemployment duration T0 . The left panel displays an increasing hazard function of ε, asa result of our log-concavity restriction. The right panel shows that the estimated baseline
unemployment hazard function roughly takes a U-shape. The initial negative duration
dependence can be a result of employers perceiving unemployment duration as a signal
about the potential productivity of the worker, because the unemployed may lose valuable
skills. On the other hand, the later positive duration dependence can be explained byfalling reservation wage and/or increasing search intensity (Burdett and Vishwanath,
1988; Van den Berg, 1990). Such U-shaped duration dependence was also documented
in empirical literature (Christofides and McKenna, 1996; McCall, 1996; Addison andPortugal, 2003).

Logconcave-error-AFT

15

Table 1. Finite sample performances of the SMLEs, sieve estimator, and Gehan-weight

estimator.
n
500

(I)
Weibull

β1

β2

1000

β1

β2

(II)
Gamma

500

β1

β2

1000

β1

β2

Bias
SE
RMSE
Bias
SE
RMSE
Bias
SE
RMSE
Bias
SE
RMSE
Bias
SE
RMSE
Bias
SE
RMSE
Bias
SE
RMSE
Bias
SE
RMSE

SMLE-uhc = 2 hc = 5-.0007 -.0017
.0960
.0944
.0960
.0944
.0075
.0055
.0970
.0954
.0973
.0956
.0021
.0035
.0703
.0704
.0703
.0705
.0049
.0030
.0696
.0698
.0698
.0698
-.0051 -.0020
.1421
.1368
.1422
.1368
.0032
.0042
.1428
.1418
.1429
.1419
.0017
.0012
.0973
.0994
0974
.0994
-.0003 -.0009
.0994
.1014
.0994
.1014

SMLE-shc = 2 hc = 5-.0007 -.0010
.0927
.0920
.0927
.0920
.0053
.0050
.0931
.0936
.0932
.0937
.0018
.0030
.0691
.0698
.0691
.0698
.0041
.0030
.0689
.0691
.0690
.0692
-.0020 -.0010
.1368
.1337
.1369
.1338
.0042
.0040
.1396
.1395
.1396
.1396
.0012
.0012
.0966
.0984
.0966
.0984
-.0004 -.0018
.0973
.0999
.0973
.0999

Gehan

Sieve

-.0041
.1021
.1021
.0056
.1043
.1044
.0016
.0762
.0762
.0067
.0752
.0755
.0036
.1551
.1551
-.0014
.1577
.1577
.0014
.1098
.1098
.0014
.1138
.1138

-.0213
.0959
.0982
.0006
.0914
.0914
-.0243
.0693
.0734
.0050
.0672
.0674
.0181
.1366
.1378
.0133
.1388
.1394
-.0253
.1065
.1094
.0031
.1053
.1054

Note: SMLE-u and SMLE-s use the unsmoothed and smoothed hazard rates, respectively. Scenario(I): The error term is the logarithm of a Weibull-distributed variate; Scenario (II): The error termis the logarithm of a Gamma-distributed variate. The bandwidths for smoothing the hazard rate arehc × n−1/5 (only used for SMLE-s) and hc × n−1/7 for its derivative (used for SMLE-u and SMLE-s).

Table 2. AFT model estimates for unemployment duration.

Covariates
Race
GenderRe-entrant

SMLE-uEst.
SE-.296 .099-.281 .098-.274 .096

SMLE-sEst.
SE-.211 .098-.164 .093-.313 .096

Gehan
Est.
SE
-.280 .185-.154 .174-.259 .188

Sieve
Est.
SE
-.366 .147-.259 .136-.280 .148

Note: The bandwidths for smoothing the estimated hazard rate is 1 × n−1/5 and for smoothing its
derivative is 1 × n−1/7 .

6. CONCLUSION
In this paper, we study the AFT model under the constraint that the error term has a
log-concave survivor function. Exploiting this shape restriction, we proposed two semi-

16

Liu and Yu

Figure 3. Estimated hazard functions of the error term ε (left) and the baseline duration
T0 (right).

parametric MLE estimators for the model coefficients. SMLE-u directly uses the NPMLEof the hazard function, which is piece-wise constant, and SMLE-s delivers smoothed estimates by resorting to the kernel smoothed version of the NPMLE. Both estimators obtain the semiparametric efficiency bound asymptotically. Simulation exercises show that
SMLE-s has better finite sample performances than SMLE-u, even though the former
introduces an additional bandwidth. Therefore, we recommend the SMLE-s to appliedresearchers. Our work makes the AFT model a more viable tool for analyzing duration
data.
ACKNOWLEDGMENT
We thank two anonymous referees for their valuable comments. Ruixuan Liu acknowledges the financial support from the Program to Enhance Research and Scholarship(PERS) at Emory University. Zhengfei Yu acknowledges the support of JSPS KAKENHI Grant Number 19K13666.
REFERENCES
Abbring, J. H. (2012). Mixed hitting-time models. Econometrica 80 (2), 783–819.
Addison, J. T. and P. Portugal (2003). Unemployment duration competing and defective
risks. Journal of Human Resources 38 (1), 156–191.
Bagnoli, M. and T. Bergstrom (2005). Log-concave probability and its applications.
Economic theory 26 (2), 445–469.
Barlow, R. and A. Marshall (1964). Bounds for distributions with monotone hazard rate
I. The Annals of Mathematical Statistics 35 (3), 1234–1257.
Bickel, P. J., C. A. Klaassen, Y. Ritov, and J. A. Wellner (1993). Efficient and adaptive
estimation for semiparametric models. Johns Hopkins University Press.
Burdett, K. and T. Vishwanath (1988). Declining reservation wages and learning. The
Review of Economic Studies 55 (4), 656–665.
Chen, h. (2007). Large sample sieve estimation of semi-nonparametric models. Handbookof Econometrics 6, 5549–5632.

Logconcave-error-AFT

17

Chetverikov, D., A. Santos, and A. M. Shaikh (2018). The econometrics of shape restrictions. Annual Review of Economics 10, 31–63.
Christofides, L. N. and C. J. McKenna (1996). Unemployment insurance and job durationin Canada. Journal of Labor Economics 14, 286–312.
Cox, D. (1972). Regression models and life-tables. Journal of Royal Statistical Society,
Series B. 24, 187–220.
Ding, Y. and B. Nan (2011). A sieve m-theorem for bundled parameters in semiparametricmodels, with application to the efficient estimation in a linear model for censored data.
The Annals of Statistics 39, 3032–3061.
Esary, J. and A. Marshall (1973). Shock models and wear processes. The Annals ofProbability 1, 627–649.
Fygenson, M. and Y. A. Ritov (1994). Monotone estimating equations for censored data.
The Annals of Statistics 22, 732–746.
Groeneboom, P. and K. Hendrickx (2018). Current status linear regression. The Annalsof Statistics 46 (4), 1415–1444.
Groeneboom, P. and G. Jongbloed (2014). Nonparametric estimation under shape constraints. Cambridge University Press.
Groeneboom, P., G. Jongbloed, and B. I. Witte (2010). Maximum smoothed likelihood
estimation and smoothed maximum likelihood estimation in the current status model.
The Annals of Statistics 38, 352–387.
Huang, J. and J. A. Wellner (1995). Estimation of a monotone density or monotonehazard under random censoring. Scandinavian Journal of Statistics, 3–33.
Ichimura, H. and S. Lee (2010). Characterization of the asymptotic distribution of semiparametric m-estimators. Journal of Econometrics 159, 252–266.
Jin, Z., D. Lin, L. Wei, and Z. Ying (2003). Rank-based inference for the accelerated
failure time model. Biometrika 90 (2), 341–353.
Kalbfleisch, J. D. and R. L. Prentice (2002). The statistical analysis of failure time data.
John Wiley & Sons.
Kim, C. and T. Lai (2000). Efficient score estimation and adaptive m-estimators incensored and truncated regression models. Statistica Sinica 10, 731–749.
Kuchibhotla, A. K., R. K. Patra, and B. Sen (2017). Efficient estimation in convex single
index models. working paper .
Lai, T. L. and Z. Ying (1991). Rank regression methods for left-truncated and rightcensored data. The Annals of Statistics 19, 531–556.
Lin, D. Y. and Z. Ying (1995). Semiparametric inference for the accelerated life modelwith time-dependent covariates. Journal of Statistical Planning and Inference 44, 47–
63.
Marshall, A. W. and I. Olkin (2007). Life distributions. Springer.
McCall, B. P. (1996). Unemployment insurance rules, joblessness, and part-time work.
Econometrica, 647–682.
Nan, B., J. D. Kalbfleisch, and M. Yu (2009). Asymptotic theory for the semiparametricaccelerated failure time model with missing data. The Annals of Statistics 37, 2351–
2376.
Nan, B. and J. A. Wellner (2013). A general semiparametric Z-estimation approach forcase-cohort studies. Statistica Sinica 23, 1155–1180.
Nolan, D. and D. Pollard (1987). U-processes: Rates of convergence. The Annals of
Statistics 15, 780–799.
Prentice, R. L. (1978). Linear rank tests with right censored data. Biometrika 65,
167–179.

18

Liu and Yu

Reid, N. (1994). A conversation with Sir David Cox. Statistical Science 9, 439–455.
Ridder, G. and T. M. Woutersen (2003). The singularity of the information matrix ofthe mixed proportional hazard model. Econometrica 71 (5), 1579–1589.
Ritov, Y. (1990). Estimating in a linear regression model with censored data. The Annalsof Statistics 18, 303–328.
Robertson, T., F. T. Wright, and R. L. Dykstra (1988). Order restricted statisticalinference. John Wiley & Sons.
Robins, J. and A. A. Tsiatis (1992). Semiparametric estimation of an accelerated failure
time model with time-dependent covariates. Biometrika 79, 311–319.
Romeo, C. (1999). Conducting inference in semiparametric duration models under inequality restrictions on the shape of the hazard implied by job search theory. Journalof Applied Econometrics 14, 587–605.
Tripathi, G. (2000). Local semiparametric efficiency bounds under shape restrictions.
Econometric Theory 16, 729–739.
Tsiatis, A. A. (1990). Estimating regression parameters using linear rank tests for censored data. The Annals of Statistics 18, 354–372.
Van den Berg, G. J. (1990). Nonstationarity in job search theory. The Review of Economic
Studies 57 (2), 255–277.
Westling, T. and M. Carone (2019). A unified study of nonparametric inference formonotone functions. Annals of Statistics, forthcoming.
Ying, Z. (1993). A large sample study of rank estimation for censored regression data.
The Annals of Statistics 21, 76–99.
Zeng, D. and D. Lin (2007). Efficient estimation for the accelerated failure time model.
Journal of the American Statistical Association 102, 1387–1396.

APPENDIX: PROOFS OF MAIN RESULTS
The Appendix proves the asymptotic results in Section 4. The cited technical lemmas
are collected in the online supplement. We first introduce some notations. Let |β| denotethe Euclidean norm for a vector β. For functional nuisance parameters ρ and η, we
consider the semi-norm: k ρ k≡ supβ∈B k ρ(t, β) k2 and k η k≡ supβ∈B k η(t, β) k2 withk · k2 denoting the L2 norm for the underlying function. Let A be the class of monotonefunctions with values in [0, M ] and C be the class of functions of bounded variation withvalues in [0, M ]. We further define functional classes F0 and F1 as follows:noF0 ≡ (y, x) 7→ I{y − x0 β ≥ s} : β ∈ B, s ∈ T ,no
F1 ≡ (y, x) 7→ xI{y − x0 β ≥ s} : β ∈ B, s ∈ T .
The corresponding convex hull of F0 and F1 are denoted as convF0 and convF1 .
Proof of Theorem 4.1: We first show the following uniform convergence resultk Ψn (β, η̂n , ρ̂n ) − Ψ(β, η0 , ρ0 ) k→p 0.

(A.1)

Observe thatk Ψn (β, η̂n , ρ̂n ) − Ψ(β, η0 , ρ0 ) k≤k (Pn − P )[∆ρ̂n (β , β){X − η̂n (β , β)}] k(A.2)+ k P [∆(ρ̂n (β , β) − ρ0 (β , β))X] k + k P [∆(ρ̂n (β , β)η̂n (β , β) − ρ0 (β , β)η0 (β , β))] k .
Note that for any β, the NPMLE estimator λ̂n (u, β) is monotonically increasing andLogconcave-error-AFT

19˙
λ̂n (u, β) is of bounded variation (see the proof of Lemma S6). Moreover, we prove inLemma S6 the following functional classnoλ̇(v − x0 β, β)0G ≡ (v, x, δ) 7→ δ
(x
−
η(v
−
x
β,
β))
:
β
∈
B,
η
∈
F̄,
λ
∈
A,
λ̇
∈
C
λ(v − x0 β, β)˙
is Glivenko-Cantelli, given that η̂n ∈ F̄ ≡ convF1 /convF0 , λ̂n ∈ A, and λ̂n ∈ C. Therefore, the first term on the right hand side of (A.2) converges to zero in outer probability.
By Lemmas S8 and S12, we have k ρ̂n (β , β) − ρ0 (β , β) k→p 0. Hence the second termon the right hand side of (A.2) is bounded byk ρ̂n (β , β) − ρ0 (β , β) k P |∆X| →p 0.
Lemma S7 implies that k η̂n (β , β) − η0 (β , β) k→p 0. Then the third term on the righthand side of (A.2) is bounded byk ρ̂n (β , β)η̂n (β , β) − ρ0 (β , β)η0 (β , β) k P |∆| ≤k ρ̂n (β , β) − ρ0 (β , β) kk η̂n (β , β) k P |∆|+ k η̂n (β , β) − η0 (β , β) kk ρ0 (β , β) k P |∆| →p 0.
We now prove the consistency of β̂n . By Lemma 4.1, β0 is the unique solution of
Ψ(β, η0 , ρ0 ) = 0. This implies that for any ε > 0, there is a δ > 0 such that
h
i
h


i
Pr | β̂n − β0 |> ε ≤ Pr Ψn β̂n , η0 (·, β̂n ), ρ0 (·, β̂n ) > δ .
Observe thatΨn β̂n , η0 (·, β̂n ), ρ0 (·, β̂n ) ≤ Ψn β̂n , η̂n (·, β̂n ), ρ̂n (·, β̂n )+ Ψn β̂n , η0 (·, β̂n ), ρ0 (·, β̂n ) − Ψn β̂n , η̂n (·, β̂n ), ρ̂n (·, β̂n ) = op (1),where the last equality follows from the definition of β̂n and (A.1). Hence β̂n →p β0 .
The proof of β̃n →p β0 is almost the same, with the only deviation of verifying
k ρ̃n (β , β) − ρ0 (β , β) k→p 0, where ρ̃n (β , β) is defined in (3.22). This immediately
follows from Lemma S12. The proof of β̂n ’s existence follows closely from Theorem 4.1
in Groeneboom and Hendrickx (2018) and is provided in S1.2 of the online supplement. 
Proof of Theorem 4.2: We present a detailed proof for β̂n and highlight the necessary changes for β̃n . First define Ψn (β̂n , η̂n (·, β̂n ), ρ̂n (·, β̂n )) = 0 following the proof of
Theorem 4.1 in Groeneboom and Hendrickx (2018): this can be done by taking convex
combinations of Ψn evaluated at the left and right limit of β̂n , given the zero-crossing
nature of β̂n . We then verify that the functional estimates η̂ and ρ̂ converge fast enoughask η̂n − η0 k= Op (n−1/2 ),k ρ̂n − ρ0 k= Op (log2 n × n−2/7 ),(A.3)

by using Lemmas S7 and S12. Note that both rates are faster than the critical rate n−1/4
in Nan and Wellner (2013). By Lemma S6, the following stochastic equicontinuity result
holds:h 
iGn ∆ ρ̂n (β̂n , β̂n ){X − η̂n (β̂n , β̂n )} − ρ0 (0 , β0 ){X − η0 (0 , β0 )} = op (1). (A.4)

20

Liu and Yu

Since Ψ(β0 , η0 (·, β0 ), ρ0 (·, β0 )) = 0, we gethin1/2 Ψ(β̂n , η̂n (·, β̂n ), ρ̂n (·, β̂n ) − Ψ(β0 , η0 (·, β0 ), ρ0 (·, β0 ))= n1/2 (Ψ − Ψn )(β̂n , η̂n (·, β̂n ), ρ̂n (·, β̂n )) + op (1)= −n1/2 (Ψn − Ψ)(β0 , η0 (·, β0 ), ρ0 (·, β0 )) + op (1),where we make use of (A.4) in the second equality. Then we take the second-order Taylor
expansion of the left hand side of the first equality, which leads to− Ψ̇β (β0 , η0 (·, β0 ), ρ0 (·, β0 ))n1/2 (β̂n − β0 )no= n1/2 (Ψn − Ψ)(β0 , η0 (·, β0 ), ρ0 (·, β0 )) + Ψ̇η (β0 , η0 (·, β0 ), ρ0 (·, β0 ))[(η̂n − η0 )(·, β0 )] + op (1).
(A.5)
The second order terms are negligible due to (A.3). The first order effect of estimationfor ρ0 is canceled:Ψ̇ρ (β0 , η0 (·, β0 ), ρ0 (·, β0 ))[(ρ̂n − ρ0 )(·, β0 )] = 0.
To prove this assertion, note that for any function g(0 , β0 ),ZP1(0 ≥ s)[X − η0 (0 , β0 )]ρ̇0β (0 , β0 )g(0 , β0 )dΛ0 (s) = 0,

(A.6)

(A.7)

because the conditional expectation of the term inside the brackets is zero by first conditioning on 1{V − X 0 β0 ≥ t}. Moreover, the following local martingale has a zero mean:ZP[X − η0 (0 , β0 )]ρ̇0β (0 , β0 )g(0 , β0 )dM (s) = 0,
(A.8)
where M (t) is the martingale associated with the counting process:
Z t
0
M (t) = ∆1(V − X β0 ≤ t) −1(V − X 0 β0 ≥ s)dΛ0 (s).

(A.9)

−∞

Therefore, the functional derivative with respect to ρ is equal to zero by summing up(A.7) and (A.8) and letting g(·, β0 ) = [(ρ̂n − ρ0 )(·, β0 )]:P {[X − η0 (0 , β0 )]ρ̇0β (0 , β0 )∆[(ρ̂n − ρ0 )(·, β0 )]} = 0.
Further simplifying (A.5) leads tonon1/2 (Ψn − Ψ)(β0 , η0 (·, β0 ), ρ0 (·, β0 )) + Ψ̇η (β0 , η0 (·, β0 ), ρ0 (·, β0 ))[(η̂n − η0 )(·, β0 )]Z= Gn ρ0 (t, β0 ){X − η0 (t, β0 )}dM (t) + op (1),R
along the same line of Nan et al. (2009). The asymptotic covariance matrix of Gn ρ0 (t, β0 ){X−
η0 (t, β0 )}dM (t) coincides with the efficient information matrix given by Lemma S1. Regarding the asymptotic analysis of β̃n , the only change occurs to the convergence rate of
the smoothed hazard function ask λ̃n (t, β) − λ0 (t, β) k= Op (n−2/5 log2 n).

(A.10)˙
Nevertheless, this does not alter the rate of ρ̃n as the convergence rate of λ̂n remains theslower and determining one.


Manuscript submitted to The Econometrics Journal , pp. S1–S16.

Online supplement of “Accelerated Failure Time
Models with Log-concave Errors”
Ruixuan Liu† and Zhengfei Yu‡†

Department of Economics, Emory University, 201 Dowman Drive, Atlanta, GA, 30322

E-mail: ruixuan.liu@emory.edu‡

Faculty of Humanities and Social Sciences, University of Tsukuba, 1-1-1, Tennodai, Tsukuba, Ibaraki, Japan.

E-mail: yu.zhengfei.gn@u.tsukuba.ac.jp
S1. PROOFS AND TECHNICAL LEMMAS
We first introduce some notations. Let |β| denote the Euclidean norm for a vector β.
For the functional nuisance parameters ρ and η, we consider the following semi-norm:k ρ k≡ supβ∈B k ρ(t, β) k2 and k η k≡ supβ∈B k η(t, β) k2 with k · k2 denoting the L2
norm for the underlying function. We denote a large positive constant by M whose value
might change line by line. We add subscripts for M if there is potential confusion. For
two sequences an , bn , an . bn , if an ≤ bn for some large M independent of n. Moreover,
we denote an  bn if an . bn and bn . an simultaneously. For any set E, linE denotes
the linear subspace generated by E and E stands for the L2 closure of E.
S1.1. Log-concavity
Proof of Proposition 2.1: (i). See Theorems 1 and 3 of Bagnoli and Bergstrom (2005).
(ii). See Corollary 2 of Bagnoli and Bergstrom (2005).
(iii). Let fT0 (t) and F̄T0 (t) be the density and the survival function of the baseline duration T0 , respectively. Since ε = log(T0 ), one immediately obtains
λ(u) =fT (exp(u)) exp(u)
f (u)
= 0
= h(exp(u)) exp(u),F̄ (u)
F̄T0 (exp(u))

which is increasing in u as both h(·) and exp(u) are increasing functions.



Proof of Proposition 2.2: It is straightforward to observe that the duration defined by
the shock process can also be written as T = T0 exp(X 0 β0 ), where the baseline duration
T0 is defined as the survival time after a random number of shocks and the arrival ofshocks is governed by a homogeneous Poisson process N (t):Pr {T0 > t} =∞kX(t) e−t
k=0

k!

P k.

The conditional distribution of T is

k
−x0 β
−x0 β
t
∞
n
o
e
t
e−e
X
−x0 β
P k.
Pr {T > t|X = x} = Pr T0 > e
t|X = x =
k!
k=0

0

The covariate effect exp(x β) enters the model by changing the arrival rate of the null Poisson process N (t). This feature is in accordance with the influence of covariates

S2

Liu and Yu

on the conditional hazard function of the duration T in AFT models, as its role is toaccelerate or decelerate the time to failure. According to Theorem 3.1 in Esary andMarshall (1973), T0 has an increasing hazard function if P k /P k−1 is decreasing in k.
The requirement stated in the proposition imply such restriction on the sequence P̄k for
Examples 2.2 and 2.3; see Theorem 5.5 and Theorem 6.1 in Esary and Marshall (1973). S1.2. Global Unique Solution

Proof of Lemma 4.1: For notation simplicity, we use Ψ(β0 , η0 , ρ0 ) to denote Ψ(β, η0 (·, β), ρ0 (·, β))whenever no confusion is caused.
It suffices to show that Ψ(β0 , η0 , ρ0 ) = 0 and (β − β0 )0 Ψ(β, η0 , ρ0 ) > 0 for all β ∈B \ {β0 }. Let Ḡx (·) ≡ 1 − GC|X (·|X = x), where GC|X is the CDF of C conditional on
X, and let H(·) be the CDF of X. Using equation (5.4) of Ying (1993), Ψ(β, η0 , ρ0 ) canbe written asZ

λ̇(t, β)λ(t, β)Z

Ḡx (t + β 0 x)f (t + (β − β0 )0 x)dH(x) q0 (t, β)dt,

whereRxḠx (t + β 0 x)f (t + (β − β0 )0 x)dH(x)q0 (t, β) ≡ R
Ḡx (t + β 0 x)f (t + (β − β0 )0 x)dH(x)RxḠx (t + β 0 x)F̄ (t + (β − β0 )0 x)dH(x)− R.
Ḡx (t + β 0 x)F̄ (t + (β − β0 )0 x)dH(x)Obviously q0 (t, β0 ) = 0. Hence Ψ(β0 , η0 , ρ0 ) = 0. Then following equation (5.5) of Ying
(1993), we obtain(β − β0 )0 Ψ(β, η0 , ρ0 ) =Zλ̇(t, β)λ(t, β)Z

Ḡx (t + β 0 x)f (t + (β − β0 )0 x)dH(x) q(t, β)dt,

where
R

(β − β0 )0 xḠx (t + β 0 x)f (t + (β − β0 )0 x)dH(x)RḠx (t + β 0 x)f (t + (β − β0 )0 x)dH(x)R(β − β0 )0 xḠx (t + β 0 x)F̄ (t + (β − β0 )0 x)dH(x)R.
Ḡx (t + β 0 x)F̄ (t + (β − β0 )0 x)dH(x)q(t, β) ≡−For x ∈ Rd , let R1 , R2 , ...., Rd be a sequence of orthogonal vectors satisfying R1 = β − β0 .
Define x̃β ≡ [x̃β,1 , x̃β,2 , ..., x̃β,d ]0 ≡ R0 x, where R ≡ [R1 , R2 , ..., Rd ]. For any β 6= β0 , were-write q(t, β) using x̃β,1 = (β − β0 )0 x and x = (R0 )−1 x̃β .
Zq(t, β) =Zx̃β,1 q1 (t, x̃β,1 )dH̃(x̃β,1 ) −x̃β,1 q2 (t, x̃β,1 )dH̃(x̃β,1 ),Online supplement

S3

where
q1 (t, x̃β,1 ) ≡ R

f (t + x̃β,1 ),f (t + x̃β,1 )dH̃(x̃β,1 )q2 (t, x̃β,1 ) ≡ RF̄ (t + x̃β,1 ),F̄ (t + x̃β,1 )dH̃(x̃β,1 )Z
dH̃(x̃β,1 ) ≡Ḡx (t + β 0 (R0 )−1 x̃β )dH((R0 )−1 x̃β ).

x̃β,2 ,...,x̃β,pBecause the log-concavity of F (·) implies an increasing hazard function f (·)/F̄ (·), the ratio of two density functions (with respect to dH̃(x̃β,1 )) q1 (t, x̃β,1 )/q2 (t, x̃β,1 ) is increasingin x̃β,1 . Since the monotone likelihood ratio implies first order stochastic dominance (seeChapter 2 of Marshall and Olkin, 2007), it follows that q(t, β) > 0 for all β 6= β0 . Again
by the log-concavity of F (·), λ̇(t, β) > 0. Hence we conclude that (β −β0 )0 Ψ(β, η0 , ρ0 ) > 0
for all β ∈ B \ {β0 }.

Proof of the existence of zero-crossing points in Theorem 4.1 The uniform
convergence in (A.1) of the Appendix leads toΨn (β, η̂n (β), ρ̂n (β)) = Ψ̇β0 (β − β0 ) + rn (β),(S.1)

where rn (β) = op (1) + o(β − β0 ). We now define for h > 0,the function
Ψn,h (β, η̂n (β), ρ̂n (β)) = Ψ̇β0 (β − β0 ) + r̃n,h (β),

(S.2)

with
−d

r̃n,h (β) = hZKh (u1 − β1 ) · · · Kh (u1 − β1 )rn (u1 , · · · , ud )du1 · · · dud ,(S.3)

where K is standard kernel density function supported on [−1, 1] and β 0 = (β1 , · · · , βd )0 .
Note that limh→0 r̃n,h (β) = rn (β). We reparameterize by defining γ = Ψ̇β0 β and γ0 =Ψ̇β0 β0 . This gives
Ψn,h (β, η̂n (β), ρ̂n (β)) = γ − γ0 + r̃n,h (Ψ̇−1β0 γ).

(S.4)

Given the result in (S.1), the mapping γ 7→ γ0 − r̃n (Ψ̇−1β0 γ) maps, for each δ > 0, theball Bδ (γ0 ) = {γ : |γ − γ0 | ≤ δ} into Bδ/2 (γ0 ) = {γ : |γ − γ0 | ≤ δ/2} with probability
approaching to 1. Therefore by Brouwer’s fixed point theorem (Groeneboom and Hendrickx (2018)), the mapping γ 7→ γ0 − r̃n,h (Ψ̇−1β0 γ) has a fixed point which we denote by−1γn,h . Let βn,h ≡ Ψ̇β0 γn,h , then we haveΨn,h (βn,h , η̂n (βn,h ), ρ̂n (βn,h )) = 0(S.5)

By compactness of B, the sequence (βn,1/k )∞k=1 must have a subsequence (βn,1/kl ) witha limit point β̄n as l → ∞. Finally, we prove as in Groeneboom and Hendrickx (2018)that Ψn (β, η̂n (β), ρ̂n (β)) has a crossing of zero at β̄n by contradiction.
Suppose that the j-th component Ψjn of Ψn does not have a crossing of zero at β̄n .
Then there must be an open ball Bδ (β̄n ) = {β : |β − β̄n | < δ} of β̄n such that Ψjn has
a constant sign in Bδ (β̄n ), say Ψjn (β, η̂n (β), ρ̂n (β)) ≥ c > 0 for all β ∈ Bδ (β̄n ) and someconstant c > 0. Arguing as Groeneboom and Hendrickx (2018), the j-th component ofΨjn,h of Ψn,h satisfiesc(S.6)
Ψjn,h (β, η̂n (β), ρ̂n (β)) ≥ ,2

S4

Liu and Yu

for sufficiently small h and all β ∈ Bδ (β̄n ), which contradicting (S.5), since βn,h for
h = 1/kl belongs to Bδ (β̄n ) for large kl .
S1.3. Semi-parametric Efficiency
The calculation of semi-parametric information bound embedding shape restrictions is
not a trivial task. Adapting the argument in Kuchibhotla et al. (2017), we formallyverify that for the increasing hazard rate as assumed in our paper, the information
bound remains unchanged as in Bickel et al. (1993). In order to determine the bound,the tangent set and the projection of the (parametric) score function to the tangent set
need to be calculated. The score function is not affected by the shape restriction. The
crux in our proof is to show the tangent set remains unchanged by showing that scoresfor smooth submodels lie in the set and by exhibiting a family of smooth submodels withscores that can approximate any element of the set arbitrarily well are dense in the set.
Note that the calculation of the semiparametric information bound takes this trimming
parameter τU as fixed; see Example 4 on page 284 of Bickel et al. (1993). Thus, the rangeof the (stochastic) integral lβ∗0 in Lemma S1 is taken to be T throughout. We suppress
T without confusion for notational simplicity.
Lemma S1. (Semiparametric efficiency bound) In model (??) with right censorship, the efficient score function for estimating β0 is!
Z
λ̇0 (t)∗0dM (t),(S.7)
lβ0 (Vi , Xi , ∆i ) = (Xi − E[Xi |Yi − β0 Xi ≥ t])
λ0 (t)
where
M (t) = ∆i 1{Vi − β00 Xi ≤ t} −Z t1{Vi − β00 Xi ≥ s}dΛ0 (s),

(S.8)

−∞

and Λ0 (·) is the cumulative hazard function of the error term ε.
Proof of Lemma S1: The complication in our setting is that the true hazard rate λ0 (t)is increasing. Recall that we assume λ0 ≥ c > 0 over its support and its derivative isbounded and strictly positive, say λ̇0 ≥ c > 0. We consider the following parametricsub-family that contains the true model:
βh1 = β0 + h1 α,log λh2 (t) = log λ0 (t) + h2 θ(t),(S.9)

where h1 and h2 are both scalars converging to zero and α is a d−dimensional vector
as in Lai and Ying (1991). The local perturbation function θ(t) is a univariate functionthat is uniformly bounded and has uniformly bounded first-order derivative θ̇(t). These
requirements are made to guarantee the monotonicity of λh2 (t) for small enough h2 .
Therefore, the score functions areZ ˙
λ0 (t)dM (t),(S.10)
lβ0 = Xλ0 (t)Zlλ0 [θ] = θ(t)dM (t),(S.11)

Online supplement

S5

where M (t) is the natural martingale in our Theorem 3.2. The efficient score function
is lβ∗0 = lβ0 − lλ0 [θ∗ ] for which lβ0 − lλ0 [θ∗ ] is orthogonal to the nuisance tangent spacedenoted by ΛS :
ΛS ≡ lin{lλ0 [θ] : λh2 (t)is increasing for small enoughh2

and θ ∈ L2 (P )} (S.12)

Obviously the L2 closure of the linear span of lλ0 [θ] where the function θ has bounded
first-order derivative (Bickel et al., 1993) is a subset of ΛS . As shown in Section 4 ofKuchibhotla et al. (2017), we getlin{lλ0 [θ] : θ

has a bounded first-order derivative} = lin{lλ0 [θ] : θ ∈ L2 (P )}, (S.13)

leading to the conclusion that ΛS is also equal to lin{lλ0 [θ] : θ ∈ L2 (P )}. The linear span
generated by those nuisance functional components agree with each other and the monotonicity of the nonparametric component does not change the efficient score calculation;
also see Tripathi (2000) for a similar setting of where the functional nuisance parameteris an increasing function in the partial linear model. Now the efficient score function is
the one obtained by Ritov and Wellner (1988):!
Z
λ̇0 (t)∗0dM (t),(S.14)
lβ0 (Vi , Xi , ∆i ) = (Xi − E[Xi |Yi − β0 Xi ≥ t])
λ0 (t)and the resulting information I(β0 ) immediately follows from Lemma S1:
I(β0 ) = E[lβ∗0 (Vi , Xi , ∆i )⊗2 ]= E [Xi Xi0 J(Ci − β00 Xi )]Z− E[Xi |Ci − β00 Xi ≥ t]E[Xi |Ci − β00 Xi ≥ t]0 E[1{Ci − β00 Xi ≥ t}]dJ(t),(S.15)
where
Z t
J(t) =−∞λ̇0 (s)λ0 (s)!2
dF0 (s).

(S.16)
S1.4. Empirical Processes
We first restate some necessary definitions and Theorem 2.4.1 in Van Der Vaart andWellner (1996) that will be used repeatedly in the sequel. Let F be the class of functionsand L2 (Q) be the L2 -norm defined by a probability measure Q . For any probability
measure Q, let N (ε, F, L2 (Q)) be the minimal number of balls of radius ε needed tocover the class F. The entropy integral J(δ, F) is defined as
Z δpJ(δ, F) ≡ sup
1 + log N (ε, F, L2 (Q))dε.
Q0

An envelope function of a functional class F is a function F such that |f (x)| ≤ F (x) for
all x and f ∈ F.

S6

Liu and Yu

Lemma S2. (Theorem 2.14.1 in Van Der Vaart and Wellner (1996)) Let P0 be
the distribution of the underlying observation and let F be a P0 -measurable class with anenvelope function F . We have
E sup |Gn f | . J(1, F) k F kP0 ,2(S.17)f ∈F

Lemma S3. (Theorem 2.7.5 in Van Der Vaart and Wellner (1996)) Let A be theclass of monotone functions with values in [0, M ], then for all δ > 0,√
(S.18)J(δ, A) . δ.
Let C be the class of functions of bounded variation with values in [0, M ], then for all
δ > 0,√J(δ, C) . δ.
(S.19)
Lemma S4. (Entropy for Convex Hull of VC Classes) Let F be a VC-class with
VC index equal to v, and denote convF be the convex hull of F, then2

J(δ, convF) . δ v+2 .

(S.20)

In our context, the following two VC classes play important roles in analyzing η̂n :noF0 ≡ (y, x) 7→ {1(y − x0 β ≥ s) : β ∈ B, s ∈ T ,(S.21)andnoF1 ≡ (y, x) 7→ x{1(y − x0 β ≥ s) : β ∈ B, s ∈ T .

(S.22)We denote the corresponding convex hull of two classes as convF0 and convF1 , respective(k)ly. Note that Dn (s, β) (k = 0, 1) are in convF0 and convF1 , see Nan and Wellner (2013).
The next lemma provides the entropy bound for an important functional class in our
remaining proofs; see also Lemma 10.1 in Groeneboom and Hendrickx (2018).
Lemma S5. Consider the following function class:2FKδ,0= {λ(x0 β) : sup |λ| ≤ K1 , |β − β0 | ≤ K2 },(S.23)

where the function λ(·) belongs to the class of monotone functions, then the following
entropy bound holds:
 M K1
2,(S.24)
H[] , FKδ,0, k·k2 ≤
for some finite constant M .
Proof of Lemma S5: For any small β , the compact neighborhood of β0 can be coveredby Nβ neighborhoods with diameters no larger than β , where Nβ ≤ M −q
β . Thus, for any
β, we can find i ∈ {1, · · · , Nβ } such that |β − βi | ≤ . For the monotone function λ, weU
can find brackets [λL
j , λj ] with size  covering the class of monotone functions with range
restricted to [−K1 , K1 ]. Moreover, the number of brackets Nλ is bounded by exp(K1 −1 )up to some finite constant.

Online supplement

S7

2Consider any function f (x) in FKδ,0, one hasf (x) ≡ λ(x0 β) = λ(x0 βi + x0 (β − βi )),(S.25)

λ(x0 βi − M β ) ≤ f (x) ≤ λ(x0 βi + M β ),(S.26)

which leads togiven that the covariates X have compact support. Therefore, we can cover the element2
in FKδ,0
by
0
U
0
λL
j (x βi − M β ) ≤ f ≤ λj (x βi + M β ),(S.27)

U
for a pair [λLj , λj ] that covers λ.
0
U0
Now we verify the size of new bracket [λL
j (x βi − M β ), λj (x βi + M β )] is less than
 up to some finite constant with a proper choice of β . We start with the following
decomposition:
0
L 0
U
0
0
k λU
j (x βi + M β ) − λj (x βi − M β ) k2 ≤ k λj (x βi + M β ) − λ(x βi + M β ) k2

+ k λ(x0 βi + M β ) − λ(x0 βi − M β ) k20
+ k λ(x0 βi − M β ) − λL
j (x βi − M β ) k2 .

Apparently, the first and third terms are bounded up by  by the construction ofU[λL
j , λj ]. Considering the second term, one get
Z 2M2k λ(x0 βi + M β ) − λ(x0 βi − M β ) k22 ≤ M(λ(t) − λ(t − 2M β )) dt−2M

by the change of variable. Now given the monotonicity of λ and the fact that it is boundedin absolute value by K, we have
Z 2M
Z 2M
2
(λ(t) − λ(t − 2M β )) dt ≤M(λ(t − 2M β ) − λ(t)) dt
−2M
−2M
"Z
#
Z
−2M

2M

λ(t − 2M )dt −=M−2M −2β Mλ(t)dt2M −2β M.β .
Then we take β =  , we get k λ(x0 βi + M β ) − λ(x0 βi − M β ) k2 . . Thus, the overallbracketing entropy number is bounded by:

M K1
M K1
2
.
.
H[] , FKδ,0, k·k2 ≤ log Nβ + log Nλ ≤ 2q log(−1 ) +
2

Now we obtain the entropy bounds for the key functional class in our context andprove the asymptotic characterizations for terms in showing consistency and asymptotic
normality. Note that the uniform convergence result in Lemma S8 implies that the estimated hazard rate function is uniformly bounded away from zero, given that the true
hazard rate function is uniformly bounded away from zero. Therefore, we can restrict
our attention to bounded monotone functions for the classes involving 1/λ(·, β).
Lemma S6. The functional class G defined bynoG ≡ (y, x, δ) 7→ δρ(y − x0 β)[x − η(y − x0 β)] : β ∈ B, η(t, β) ∈ F̄, ρ ∈ H(S.28)

S8

Liu and Yu

has bounded entropy integral, where F̄ ≡ convF1 /convF0 and H ≡ C/A. Therefore, wehave the following Glivenko-Cantelli resultihPn ∆ρ̂n (β̂n , β̂n ){X − η̂n (β̂n , β̂n )} − P [ρ0 (0 , β0 ){X − η0 (0 , β0 )}] = op (1),and the stochastic equicontinuity ash 
iGn ∆ ρ̂n (β̂n , β̂n ){X − η̂n (β̂n , β̂n )} − ρ0 (0 , β0 ){X − η0 (0 , β0 )} = op (1).
Proof of Lemma S6: We first verify that the uniform entropy integral J(1, G) isbounded. We consider the following three subclasses:noG1 ≡ (y, x, δ) 7→ 1/λ(y − x0 β, β) : β ∈ B, λ(t, β) ∈ A ,(S.29)noG2 ≡ (y, x, δ) 7→ λ̇(y − x0 β, β) : β ∈ B, λ̇(t, β) ∈ C ,(S.30)noG3 ≡ (y, x, δ) 7→ δ[x − η(y − x0 β)] : β ∈ B, η(t, β) ∈ F̄ .
(S.31)
Because the NPMLE estimator λ̂n (t, β) is an increasing function for any given β, G1 is theclass involving monotonically decreasing functions. Hence, the uniform entropy integral
J(1, G1 ) is bounded. A similar argument applies to G2 because the kernel estimator of its
˙
density function is λ̂n (t, β) isZ˙λ̂n (t, β) = Kh (t − u)dλ̂n (u, β),for a continuous kernel function K and a function of bounded variation λ̂n (since it
˙
is increasing). Therefore, λ̂n (t, β) is of bounded variation by Theorem I.5.c in Widder
˙(1941). The composition of λ̂n or λ̂n with the linear index y − x0 β still has the bounded(k)entropy integral as shown in our Lemma S5. Moreover, Dn (t, β) and D(k) (t, β) are in
the convex hull of certain VC classes for k = 0, 1. Thus, the resulting uniform entropy
integral is bounded. Thus, the conclusion for the whole functional class G = G1 · G2 · G3
follows from Example 2.10.8 Van Der Vaart and Wellner (1996).
Thereafter, we show
k P [δ{ρ[X − η] − ρ0 [X − η0 ]}] k≤k P [δX(ρ − ρ0 )] k + k P [δ(ρη − ρ0 η0 )] k→ 0,
as |β − β0 | → 0, k η − η0 k→ 0, and k ρ − ρ0 k→ 0.
Hence, by (S.17) we get the desired stochastic equicontinuity as in (??) of Appendix A
in the paper.
We then present the linear expansion of η̂n in Lemma S7, which has been establishedby Nan and Wellner (2013). We refer readers to [page 1172] of their paper for the proof.
Lemma S7. Suppose that Assumptions 4.1 to 4.6 hold. Then we have the following expansion
n1/2 [η̂n (t, β) − η0 (t, β)] = D(0) (t, β)−1 Gn [1(β ≥ t){X − η0 (t, β)}] + op (1).
Therefore, k η̂n − η0 k= Op (n−1/2 ).

(S.32)

Online supplement

S9

S1.5. Convergence Rate of Estimators of the Hazard Function and Its Derivative
This section is targeted to show the rate of convergence of ρ̂n as k ρ̂n − ρ0 k= Op (log2 n ×
˙
n−2/7 ). Given the ratio form of ρ̂n , the crux is to determine the rate of λ̂n . We introduce
additional notations. Let λ̂n be the NPMLE and {τ1 , ..., τm } with m ≤ n be its jump
points. Denote intervals Ji ≡ [τi , τi+1 ) for 0 ≤ i ≤ m − 1. For u ∈ Ji , define
if ∀t ∈ Ji : λ0 (t, β) > λ̂n (τi , β)τi ,Ân (u, β) = s,if ∃s ∈ Ji : λ0 (s, β) = λ̂n (s, β)τi+1 ,
if ∀t ∈ Ji : λ0 (t, β) < λ̂n (τi , β).
We first state three lemmas regarding the NPMLE of the monotone hazard function.
Their proofs are in analog with the ones in Groeneboom et al. (2010) and are hence omitted. Lemma S8 concerning the properties of the NPMLE for the increasing hazard rate
λ0 (t, β). The proof is in analog with Equations (2.1), (A.20), and (A.21) in Groeneboom
et al. (2010).
Lemma S8. Suppose that Assumptions 4.1 to 4.6 hold. Then we havePr lim k λ̂n (u, β) − λ0 (u, β) k= 0 = 1,n→∞k λ̂n (u, β) − λ0 (u, β) k= Op log n × n−1/3 ,k Ân (u, β) − u k= Op log n × n−1/3 .

(S.33)

The following Lemma is in analogy with Lemma 4.1 of Groeneboom et al. (2010) .
Lemma S9. Suppose that Assumptions 4.1 to 4.6 hold. For any given t and β, the following representations hold:
Z
Z(S.34)Kh̃ (t − u)d(λ̂n − λ0 )(u, β) = − ψh̃,t (λ̂n (u, β) − λ0 (u, β))du,andZ

Z
Kh (t − u)d(λ̂n − λ0 )(u, β) = −φh,t (λ̂n (u, β) − λ0 (u, β))du,(S.35)

φh,t (u) = kh (t − u)/p(v, x, δ),(S.36)p(v, x, δ) = f δ (v − x0 β0 )Ḡδx (v)F̄ 1−δ (v − x0 β0 )gx1−δ (v)dH(x),(S.37)whereψh̃,t (u) = Kh̃ (t − u)/p(v, x, δ),andwith kh (·) ≡ K 0 (·/h)/h2 andis the joint density function of (V, X, ∆).
Consider the piece-wise constant version of φh,t which is constant on the same intervals
where the NPMLE λ̂n (·, β) remains constant. For u ∈ Ji , define
φ̄h,t (u) = φh,t (Ân (u, β))and ψ̄h̃,t (u) = ψh̃,t (Ân (u, β)).

(S.38)

S10

Liu and Yu

The next lemma spells out the difference between the functions ψh̃,t , φh,t and their piecewise constant approximations ψ̄h̃,t , φ̄h,t . This lemma is in analogy with Lemma A.4 of
Groeneboom et al. (2010).
Lemma S10. Suppose that Assumptions 4.1 to 4.6 hold. Then we have|ψh̃,t (u) − ψ̄h̃,t (u)| .

|λ̂n (u, β) − λ0 (u, β)|1{|t − u| ≤ h̃},h̃2(S.39)|φh,t (u) − φ̄h,t (u)| .

|λ̂n (u, β) − λ0 (u, β)|1{|t − u| ≤ h}.
h3

(S.40)

andAn immediate consequence of Lemma S10 is the characterization of some negligible termsin Lemma S11, which in turn will be used in the proof of Lemma S12. Recall that T is
the support of ε trimmed from the right tail at come given τU as in Tsiatis (1990).
Lemma S11. Suppose that Assumptions Assumptions 4.1 to 4.6 hold. Then we haveZ(ψh̃,t − ψ̄h̃,t (u))(λ̂n (u, β) − λ0 (u, β))dP = op (log2 n × n−2/5 ),(S.41)
andZ(φh,t − φ̄h,t (u))(λ̂n (u, β) − λ0 (u, β))dP = op (log2 n × n−2/7 ),(S.42)

uniformly over t ∈ T , β ∈ B.
Proof of Lemma S11: We only prove the second claim to avoid repetition. We apply(S.40) in Lemma S10 and the Cauchy-Schwarz inequality to obtain:Z(φh,t − φ̄h,t (u))(λ̂n (u, β) − λ0 (u, β))dPZ−3.h(λ̂n (u, β) − λ0 (u, β))2 dP
u[∈t−h,t+h]−3+1/2.hk λ̂n (u, β) − λ0 (u, β) k2 = Op (n5/14−2/3 log2 n) = op (log2 n × n−2/7 ),which leads to the desired result.

˙
Now we state and prove the key lemma concerning the rate of convergence of λ̂n and
λ̃n . The convolution type estimation in (3.16) and (3.20) (of the paper) appear similar to kernel density estimation. However, the asymptotic analysis is more complicated(Groeneboom et al., 2010), because λ̂n is not linear, in contrast with the empirical distribution function.
Lemma S12. Suppose that Assumptions 4.1 to 4.6 hold. Then, for h  n−1/7 ,˙k λ̂n (t − x0 β, β) − λ̇0 (t − x0 β, β) k= Op n−2/7 log2 n .

(S.43)

Moreover, for h̃  n−1/5 ,k λ̃n (t − x0 β, β) − λ0 (t − x0 β, β) k= Op n−2/5 log2 n .

(S.44)

Online supplement

S11

Proof of Lemma S12: We only prove (S.43) which is in fact the harder one. The
deterministic bias term can be easily dealt with based on standard kernel smoothingarguments:
Z
Kh (t − u)dλ0 (u, β) − λ0 (t, β) = O(h2 ) = O(n−2/7 ).
supβ,tWhen it comes to the stochastic part, we rely on the representation in (S.42) to obtainZ
Z
Kh (t − (y − x0 β))d(λ̂n − λ0 )(y − x0 β, β) = φh,t [λ̂n (y − x0 β, β) − λ0 (y − x0 β, β)]dPZ= φ̄h,t (λ̂n (y − x0 β, β) − λ0 (y − x0 β, β))dPZ+ (φh,t − φ̄h,t )(λ̂n (y − x0 β, β) − λ0 (y − x0 β, β))dP ≡ I1n + I2n .
The second term I2n in the last equality is of a negligible order Op (log n × n−2/3 ) as
shown in Lemma S11.
Given the fact that λ̂n (y − x0 β, β) is the greatest convex majorant of the cumulative
sum diagram and the function φ̄h,t is constant on the same intervals where the NPMLE
remains constant, we get"#Z
1{β ≤ v − x0 β, δ = 1}0
φ̄h,t λ̂n (v − x β, β) −dPn (v, x, δ) = 0,
(S.45)(0)Dn (v − x0 β, β)see equation (11.81) on [page 348] of Groeneboom and Jongbloed (2014). Therefore, we
insert (S.45) into I1n :Z
I1n = φ̄h,t (λ̂n (y − x0 β, β) − λ0 (y − x0 β, β))dP#"Z
1{β ≤ y − x0 β, δ = 1}0
− φ̄h,t λ̂n (y − x β, β) −dPn
(0)Dn (y − x0 β, β)Z
1{β ≤ y − x0 β, δ = 1}
= − φ̄h,t λ̂n (y − x0 β, β) −d(Pn − P )D(0) (y − x0 β, β)
#
"
Z
1
1
0
− (0)
dPn
+ φ̄h,t 1{β ≤ y − x β, δ = 1}(0)0
Dn (y − x0 β, β) D (y − x β, β)≡ J1n + J2n .

(S.46)In the above decomposition, we have made use of the fact that
λ0 (y − x0 β, β)d(y − x0 β) =dN (y − x0 β, β),D(0) (y − x0 β, β)

which leads to
λ0 (y − x0 β, β)dP =1{β ≤ y − x0 β, δ = 1}
dP.
D(0) (y − x0 β, β)S12Liu and Yu

The analysis of J2n proceeds as follows
#
"
Z
1
1
0
J2n = φ̄h,t 1{β ≤ y − x β, δ = 1}d(Pn − P )
− (0)
(0)
0
Dn (y − x0 β, β) D (y − x β, β)
"
#
Z
1
1
+ φ̄h,t 1{β ≤ y − x0 β, δ = 1}
dP
− (0)
0
(0)
Dn (y − x0 β, β) D (y − x β, β)
(a)

(b)

≡ J2n + J2n .

(S.47)

We prove in Lemma S14 that
(a)

(b)

J2n = op (n−2/7 log2 n), and J2n = Op (n−1/2 × h−1 ) = Op (n−5/14 ).
Referring to the leading term J1n , we getZ
1{β ≤ y − x0 β, δ = 1}0−J1n = φ̄h,t λ0 (y − x β, β) −d(Pn − P )D(0) (y − x0 β, β)Z+ φ̄h,t λ̂n (y − x0 β, β) − λ0 (y − x0 β, β) d(Pn − P ).
The second term is denoted by Rn term and we prove that it is of smaller order in Lemma
S13. Finally, we utilize Lemma A.7 in Groeneboom et al. (2010) to conclude:Z
1{β ≤ y − x0 β, δ = 1}0d(Pn − P )φ̄h,t λ0 (y − x β, β) −D(0) (y − x0 β, β)Z
1{β ≤ y − x0 β, δ = 1}0
= φh,t λ0 (y − x β, β) −d(Pn − P )D(0) (y − x0 β, β)Z1{β ≤ y − x0 β, δ = 1}0d(Pn − P )+φ̄h,t − φh,t λ0 (y − x β, β) −D(0) (y − x0 β, β)= Op n−2/7 log2 n .
The proof of (S.44) is very similar with some notation-wise difference. When we take(a)(b)
h̃  n−1/5 , those smaller order terms as Rn , J2n and J2n can be shown of order2−2/5
op (log n × n). Meanwhile, the leading term as in J1n is of order Op (log2 n × n−2/5 ),following Lemma A.7 in Groeneboom et al. (2010), mutatis mutandis.
Compared with the point-wise result in Groeneboom et al. (2010), The result of Lemma
S12 contains extra log2 n terms in the rate because we also need the uniform convergenceover t ∈ T for the argument of the kernel functions. Nevertheless, the obtained rates arefast enough, i.e., of op (n−1/4 ). The following two lemmas characterize a few smaller orderterms in the proof of Lemma S12.
Lemma S13. Suppose that Assumptions 4.1 to 4.6 hold. Then we have the following
characterizationZRn ≡ φ̄h,t λ̂n (y − x0 β, β) − λ0 (y − x0 β, β) d(Pn − P ) = op (log2 n × n−2/7 ),uniformly over β ∈ B and t in the support T .
Proof of Lemma S13: We first introduce some notations adapted from Groeneboom

Online supplement

S13

et al. (2010). Define for
Ct,n (u) =k(n1/7 (t − u)/c)1{t − cn−1/7 ≤ u ≤ t + cn−1/7 },cp(y, x, δ)(S.48)

the function
ξA,B,t,n (u) = Ct,n (A(u))B(u)(S.49)

where t ∈ T , A is an increasing function, and B is a function of bounded variation. And
let
G1,n ≡ {ξA,B,t,n (u) : A ∈ A, B ∈ B, t ∈ T }.

(S.50)

Given the result in (S8), for any small γ > 0 we can find finite constant term M1 suchthat for all n sufficiently large:Pr{Υ1,n,M1 } ≡ Pr{sup |λ̂n (u, β) − λ(u, β)| ≤ M1 n−1/3 log n} ≥ 1 − γ/2.
u,βNow for the vanishing sequence νn to be specified later, it is straightforward to arrive at
Pr{|n2/7 Rn | > νn } = Pr{|n2/7 Rn | > νn ∩ Υ1,n,M1 } + Pr{|n2/7 Rn | > νn ∩ Υc1,n,M1 }≤ νn−1 E|n2/7 Rn |1{Υ1,n,M1 } + γ/2.
Again by (S8), we have
E|n

2/7

Rn |1{Υ1,n,M1 } ≤ E

supn2/7−1/3Zlog nφh,t (A(u))B(u)d(Pn − P )A∈A,B∈BM1≤ n4/7−5/6 log nE supZ
ξ(u)dGn (u) .

ξ∈G1,nThe rest of the proof is to utilize Theorem 2.14.1 in Van Der Vaart and Wellner (1996) tobound the expectation in the last display. Following the construction in Groeneboom et al.
(2010), we can select a minimal n−1/7 δ/(4M1 )-net in A and a minimal δ/(2 k Cn k∞ )-net
in B. Referring to the part involving the kernel function k(·) indexed by t, we can obtaina minimal n−1/7 δ/(4M1 )-net in T given the fact that the functional class is of VC-type
Nolan and Pollard (1987). The number of functions in this net to cover G1,n is boundedby M n1/7 /δ exp(n1/7 /δ), so the entropy integral is bounded above byJ(1, G1,n ) . n1/14 log n.

(S.51)

Compared with Groeneboom et al. (2010), the uniformity over t ∈ T in the kernel
functionbrings an extra log n term. The L2 -norm of the envelope function is of order√
Op ( h) = Op (n−1/14 ) by standard arguments in kernel smoothing. Thus applying (S.17),
we haveZ
E|Rn | ≤n4/7−5/6 log n × E sup
ξ(u)dGn (u) . n4/7−5/6 log2 n,ξ∈G1,nwhich immediately leads to Rn = op (log2 n × n−2/7 ).

(a)Lemma S14. Suppose that Assumptions Assumptions 4.1 to 4.6 hold. Then J2n and(b)
J2n defined by (S.47) in Appendix B satisfy:(a)J2n = op (log2 n × n−2/7 ),

and

(b)

J2n = Op (n−1/2 × h−1 ) = op (n−5/14 ),S14

Liu and Yu

uniformly over t ∈ T , β ∈ B.
(0)

(a)

Proof of Lemma S14: After linearizing 1/Dn , the first term J2n is determined byZ(0)
D(0) (y − x0 β, β) − Dn (y − x0 β, β)d(Pn − P ),Rn(a) ≡ φ̄h,t 1{β ≤ y − x0 β, δ = 1}
D(0)2 (y − x0 β, β)(S.52)
(0)
−1 −2
(0)
modulo a second order term which is of Op (n h ). Note that both Dn and D aremonotonically decreasing functions as defined in (??), hence their difference is a function(a)of bounded variation. Now it is clear that Rn can be dealt with in a completely analogous(a)
way as the term Rn in the previous lemma, which leads to Rn = op (log2 n × n−2/7 ).
(b)
Referring to J2n , we get
#
"
Z
1
1
0
dP
− (0)
φ̄h,t 1{β ≤ y − x β, δ = 1}0(0)
Dn (y − x0 β, β) D (y − x β, β)Z
φ̄h,t 1{β ≤ y − x0 β, δ = 1}
(0)
(0)
dP
≤k Dn − D k∞(0)
Dn (y − x0 β, β)D(0) (y − x0 β, β). n−1/2 × h−1 = Op (n−5/14 ).


S2. AN ADDITIONAL EMPIRICAL APPLICATION
Methadone maintenance is a drug replacement therapy that uses methadone over a prolonged period of time as a treatment for heroin addicts. However, hospitals and clinics
often have problems to retain the addicts. Barnett et al. (2000) finds that the medianwilling to pay of heroin addicts in Baltimore was much lower than the estimated cost of
the methadone maintenance program. Meinhofer and Witman (2018) discover that the
recent expansion of health insurance in US (Medicaid expansion) substantially increased
admissions to medication-assistant treatment1 for opioid use disorder in outpatient settings while hardly changed admissions to inpatient settings. Here we use the data of anAustralian methadone program (Caplehorn and Bell, 1991) to investigate the retention
time of heroin addicts in methadone maintenance. The data was collected from 238 addicts who were assigned to one of two public clinics based on their home address between
February 1986 and August 1987. Both clinics serve adjacent areas of Sydney with similarsocioeconomic characteristics. The censoring rate is 37%. In AFT model, we regress thedays in treatment, censored at the end of the study, to the treatment variable is “Clinic”(taking value 1 or 2 to indicate two clinics that take different overall treatment polices(Kleinbaum and Klein, 2011), and two other variables “Prison” (that indicates prison
record status) and “Dose” (that records the maximum methadone dose per day in mg).
Estimates of parameters and standard errors are reported in Table S1. All methods suggest that Clinic 2 managed to retain addicts for a substantially longer period of time. The
negative association between the prison record status and retention time is significant
at 5% level for the proposed SMLE methods but insignificant for other methods. The
1 The treatment uses methadone, buprenorphine, or naltrexone.

Online supplement

S15

positive correlation between the maximum daily methadone and the retention time ishighly significant for all methods.
Table S1. AFT model estimates for retention time in methadone maintenance.

Variables
ClinicPrisonDose

SMLE-uEst.
SE.521 .121-.219 .108.033 .004

SMLE-sEst.
SE.458 .124-.174 .065
.014 .003

Gehan
Est.
SE
.510 .203-.304 .158
.032 .007

Sieve
Est.
SE
.629 .178-.277 .153
.058 .005

Note: The bandwidths for smoothing the estimated hazard rate is 1 × n−1/5 and for smoothing its
derivative is 1 × n−1/7 .

Figure S1 displays the smoothed estimates of the hazard function of the error term εand the hazard of T0 in the retention of addicts data. By the log-concavity assumption on
the error term, its estimated hazard function is increasing, as the left panel shows. The
right panel exhibits that the overall trend of the estimated hazard function is increasingin the addicts’ retention data. This suggest that the risk of patients dropping out of
the methadone treatment keeps increasing over time, given the prison record status andmaximum daily dose.

Figure S1. Estimated hazard functions of the error term ε (left) and the baseline duration
T0 (right).

S16

Liu and Yu

REFERENCES
Bagnoli, M. and T. Bergstrom (2005). Log-concave probability and its applications.
Economic theory 26 (2), 445–469.
Barnett, P. G., S. S. Hui, et al. (2000). The cost-effectiveness of methadone maintenance.
Mount Sinai Journal of Medicine 67 (5-6), 365–374.
Bickel, P. J., C. A. Klaassen, Y. Ritov, and J. A. Wellner (1993). Efficient and adaptive
estimation for semiparametric models. Johns Hopkins University Press.
Caplehorn, J. R. and J. Bell (1991). Methadone dosage and retention of patients maintenance treatment. Medical Journal of Australia 154 (3), 195–199.
Esary, J. and A. Marshall (1973). Shock models and wear processes. The Annals ofProbability 1, 627–649.
Groeneboom, P. and K. Hendrickx (2018). Current status linear regression. The Annalsof Statistics 46 (4), 1415–1444.
Groeneboom, P. and G. Jongbloed (2014). Nonparametric estimation under shape constraints. Cambridge University Press.
Groeneboom, P., G. Jongbloed, and B. I. Witte (2010). Maximum smoothed likelihood
estimation and smoothed maximum likelihood estimation in the current status model.
The Annals of Statistics 38, 352–387.
Kleinbaum, D. G. and M. Klein (2011). Survival Analysis: A Self-Learning Text. Springer.
Kuchibhotla, A. K., R. K. Patra, and B. Sen (2017). Efficient estimation in convex single
index models. working paper .
Lai, T. L. and Z. Ying (1991). Rank regression methods for left-truncated and rightcensored data. The Annals of Statistics 19, 531–556.
Marshall, A. W. and I. Olkin (2007). Life distributions. Springer.
Meinhofer, A. and A. E. Witman (2018). The role of health insurance on treatment
for opioid use disorders: Evidence from the affordable care act medicaid expansion.
Journal of health economics 60, 177–197.
Nan, B. and J. A. Wellner (2013). A general semiparametric Z-estimation approach forcase-cohort studies. Statistica Sinica 23, 1155–1180.
Nolan, D. and D. Pollard (1987). U-processes: Rates of convergence. The Annals of
Statistics 15, 780–799.
Ritov, Y. and J. A. Wellner (1988). Censoring, martingales, and the cox model. Contemporary Mathematics 80, 191–219.
Tripathi, G. (2000). Local semiparametric efficiency bounds under shape restrictions.
Econometric Theory 16, 729–739.
Tsiatis, A. A. (1990). Estimating regression parameters using linear rank tests for censored data. The Annals of Statistics 18, 354–372.
Van Der Vaart, A. and J. A. Wellner (1996). Weak convergence and empirical processes.
Springer.
Widder, D. V. (1941). The Laplace transform. Princeton University Press.
Ying, Z. (1993). A large sample study of rank estimation for censored regression data.
The Annals of Statistics 21, 76–99.