Boosted nonparametric hazards with time-dependent covariates
Donald K. K. Lee Ningyuan Chen Hemant Ishwaran
















textfont=bf
[algorithm]labelfont=sc,font=large



























ruled
algorithmtbploa

algorithm






=30pt =5pt
byby
=12pt  plus1pt minus1pt

	<
=




#1<#1. 

#1<#1
#1section
	.#1.
#1subsection
	.#1.

	



claim






.5em
#1 #2. (#3)
claim


theoremTheorem
propositionProposition
lemmaLemma
corollaryCorollary
remarkRemark
exampleExample
definitionDefinition

assumptionAssumption


assumptionxassumption



	

	

	

	

	

	








	

    #1(<ref>)#1    #1
            
            
            αβ̱δ̣ΔεγΓκ̨łλŁΛμϕΦρ̊σθΘτζ
            
            
            å A B C D E F
             H I J K
             M
            𝒩 O P Q Rß S T U V W X Y Z
            
            
            
            
            
             argmax
             argmin
            
            
            𝔼ℙℝℚ
            
            
            
            AGEboostCENSUSESINUISANCEFailF̃F
            λ̂_boostλ_nm̂μ_n
            #1‖ #1 ‖_Φ#1#2‖ #1 ‖_Φ,#2#1#2‖ #1 ‖^2_Φ,#2R
            φα_
            
            
            =cmcsc10
            =cmcsc10 scaled  1
            =cmcsc9 scaled  2
            =cmsl12
            =cmbx9 scaled  1
            
            
            .gif,.pdf,.png,.jpg,.tiff
            
            
            empty
             Boosted nonparametric hazards with time-dependent covariates10pt
            	 Donald K.K. Lee[Correspondence: donald.lee@emory.edu. Supported by a hyperplane], Ningyuan Chen[Supported by the HKUST start-up fund R9382], Hemant Ishwaran[Supported by the NIH grant R01 GM125072]5pt
            	 Emory University, University of Toronto, University of Miami10pt
            	
            
            
            
            	Preprint of https://doi.org/10.1214/20-AOS2028Annals of Statistics 49:4:2101-2128 (2021)
            
            
            Given functional data from a survival process with time-dependent covariates, we derive a smooth convex representation for its nonparametric log-likelihood functional and obtain its functional gradient. From this we devise a generic gradient boosting procedure for estimating the hazard function nonparametrically.
            
            An illustrative implementation of the procedure using regression trees is described to show how to recover the unknown hazard.  The generic estimator is consistent if the model is correctly specified; alternatively an oracle inequality can be demonstrated for tree-based models. To avoid overfitting, boosting employs several regularization devices. One of them is step-size restriction, but the rationale for this is somewhat mysterious from the viewpoint of consistency. Our work brings some clarity to this issue by revealing that step-size restriction is a mechanism for preventing the curvature of the risk from derailing convergence.
            MSC 2010 subject classifications. Primary 62N02; Secondary
            	62G05, 90B22.
            Keywords. survival analysis, gradient boosting,
            	functional data, step-size shrinkage, regression trees, likelihood
            	functional.1pt
            
            
            
            
            § INTRODUCTION
            
            Flexible hazard models involving time-dependent covariates are
            indispensable tools for studying systems that track covariates over
            time.  In medicine, electronic health records systems make it possible
            to log patient vitals throughout the day, and these measurements can
            be used to build real-time warning systems for adverse outcomes such
            as cancer mortality <cit.>. In financial technology,
            lenders track obligors' behaviours over time to assess and revise
            default rate estimates.  Such models are also used in many other
            fields of scientific inquiry since they form the building blocks for
            transitions within a Markovian state model.  Indeed, this work was
            partly motivated by our study of patient transitions in emergency
            department queues and in organ transplant waitlist
            queues <cit.>.  For example, allocation for a donor heart
            in the U.S. is defined in terms of coarse tiers
            <cit.>, and transplant candidates are assigned to
            tiers based on their health status at the time of listing.  However, a
            patient's condition may change rapidly while awaiting a heart, and
            this time-dependent information may be the most predictive
            of mortality and not the static covariates collected far in the past.
            
            The main contribution of this paper is to introduce a fully
            nonparametric boosting procedure for hazard estimation with
            time-dependent covariates.  We describe a generic gradient boosting
            procedure for boosting arbitrary base learners for this
            setting. Generally speaking, gradient boosting adopts the view of
            boosting as an iterative gradient descent algorithm for minimizing a
            loss functional over a target function space.  Early work includes
            Breiman <cit.> and Mason et al. <cit.>.  A unified
            treatment was provided by Friedman <cit.>, who coined the term
            “gradient boosting” which is now generally taken to be the modern
            interpretation of boosting.
            
            Most of the existing boosting approaches for survival data focus on
            time-static covariates and involve boosting the Cox proportional hazards model. Examples
            include the popular R-packages  mboost (Bühlmann and Hothorn <cit.>)
            and  gbm (Ridgeway <cit.>) which apply gradient boosting to the Cox partial likelihood loss.
            Related work includes the penalized Cox partial likelihood approach of Binder and Schumacher <cit.>.  Other important approaches, but not based on
            the Cox model, include L_2Boosting <cit.> with inverse
            probability of censoring weighting
            (IPCW) <cit.>, boosted
            transformation models of parametric families <cit.>, and
            boosted accelerated failure time models <cit.>.
            
            While there are many boosting methods for dealing with time-static
            covariates, the literature is far more sparse for the case of
            time-dependent covariates.  In fact, to our knowledge there is no
            general nonparametric approach for dealing with this setting.  This is
            because in order to implement a fully nonparametric estimator, one has
            to contend with the issue of identifying the gradient, which turns out
            to be a non-trivial problem due to the functional nature of the data.
            This is unlike most standard applications of gradient boosting where
            the gradient can easily be identified and calculated.
            
            Time-dependent covariate framework
            To explain why this is so challenging, we start by formally defining the
            survival problem with time-dependent covariates.  Our description
            follows the framework of Aalen <cit.>.  Let T denote the
            potentially unobserved failure time. Conditional on the
            history up to time t- the probability of failing at T ∈ [t,t+dt)
            equals
            λ(t,X(t))Y(t)dt.
            
            Here λ(t,x) denotes the unknown hazard function,
            X(t)∈⊆^p is a predictable covariate process, and
            Y(t)∈{0,1} is a predictable indicator of whether the subject is
            at risk at time t.[The filtration of interest is σ{X(s),Y(s),I(T≤ s):s≤ t}. If X(t) is only observable when
            	Y(t)=1, we can set X(t)=x^c∉ whenever Y(t)=0.] To
            simplify notation, without loss of generality we normalize the units
            of time so that Y(t)=0 for t>1.[Since the data is always observed up to some finite time, there is no information loss from censoring at that point. For example, if T' is
            	the failure time in minutes and the longest duration in the data is
            	τ'=60 minutes, the failure time in hours, T, is at most
            	τ=1 hour. The hazard function on the minute timescale,
            	λ_T'(t',X(t')), can be recovered from the hazard function
            	on the hourly timescale, λ_T(t,X(t)), via
            	λ_T'(t',X(t')) =
            	1/τ'λ_T(t'/τ',X(t'/τ')).]
            In other words, the subject is not at risk after time t=1, so we can
            restrict attention to the time interval (0,1].
            
            If failure is observed at T∈(0,1] then the indicator Δ=Y(T)
            equals 1; otherwise Δ=0 and we set T to an arbitrary number larger than 1, e.g. T=∞. Throughout we
            assume we observe n independent and identically distributed
            functional data samples
            {(X_i(·),Y_i(·),T_i)}_i=1^n. The evolution of
            observation i's failure status can then be thought of as a
            sequence of coin flips at time increments t = 0,dt,2dt,⋯,
            with the probability of “heads” at each time point given
            by (<ref>). Therefore, observation i's contribution to
            the likelihood is
            
                [ {1-λ(0,X_i(0))Y_i(0)dt}×{1-λ(dt,X_i(dt))Y_i(dt)dt}×⋯×λ(T_i,X_i(T_i))^Δ_i;                         e^-∫_0^1 Y_i(t)λ(t,X_i(t))dtλ(T_i,X_i(T_i))^Δ_i, ]
            
            where the limit can be understood as a product integral. Hence, if the
            log-hazard function is
            
                F(t,x)=logλ(t,x),
            
            then the (scaled) negative log-likelihood functional is
            _n(F)
            =1/n∑_i=1^n∫_0^1 Y_i(t)e^F(t,X_i(t))dt
            -1/n∑_i=1^nΔ_iF(T_i,X_i(T_i)) ,
            
            which we shall refer to as the likelihood risk. The goal is to
            estimate the hazard function ł(t,x)=e^F(t,x) nonparametrically by minimizing _n(F).
            
            
            The likelihood does not have a gradient in generic function spaces
            As mentioned, our approach is to boost F using functional gradient
            descent.  However, the chief difficulty is that the canonical
            representation of the likelihood risk functional does not have a
            gradient.  To see this, observe that the directional derivative
            of eq:loglik equals
                  -25pt
            d/dθ_n(F+θ f)|_θ=0
                   = 1/n∑_i=1^n∫_0^1 Y_i(t)e^F(t,X_i(t)) f(t,X_i(t))dt
            - 1/n∑_i=1^nΔ_if(T_i,X_i(T_i)),
            
            which is the difference of two different inner products ⟨
            e^F,f⟩ _†-⟨ 1,f⟩
            _ where
            
                ⟨ g,f ⟩_†    = 1/n∑_i=1^n ∫_0^1 Y_i(t) g(t,X_i(t)) f(t,X_i(t))dt, 
                ⟨ g,f ⟩_    = 1/n∑_i=1^n Δ_i g(T_i,X_i(T_i))f(T_i,X_i(T_i)).
            
            Hence, (<ref>) cannot be expressed as a single inner
            product of the form ⟨ g_F, f ⟩ for some function
            g_F(t,x).  Were it possible to do so, g_F would then be the gradient function.
            
            In simpler non-functional data settings like regression or
            classification, the loss can be written as L(Y,(x)), where
             is the non-functional statistical target and Y is the
            outcome, so the gradient is simply ∂ L(Y,(x))/∂(x).  The negative gradient is then approximated using a base
            learner f∈ from a predefined class of functions  (this
            being either parametric; for example linear learners, or
            nonparametric; for example tree learners).  Typically, the optimal
            base learner f̂ is chosen to minimize the L^2-approximation
            error and then scaled by a regularization parameter 0<ν≤ 1 to
            obtain the updated estimate of :
            
                ← - νf̂, 15pt f̂=_f∈‖∂ L/∂-f ‖_2.
            
            Importantly, in the simpler non-functional data setting the gradient does
            not depend on the space that  belongs to.  By contrast, a key
            insight of this paper is that the gradient of _n(F) can only
            be defined after carefully specifying an appropriate sample-dependent
            domain for _n(F). The likelihood risk can then be
            re-expressed as a smooth convex functional, and an analogous
            representation also exists for the population risk. These
            representations resolve the difficulty above, allow us to describe and
            implement a gradient boosting procedure, and are also crucial to
            establishing guarantees for our estimator.
            
            
            Contributions of the paper
            A key discovery that unlocks the boosted hazard estimator is Proposition <ref> of Section <ref>. It provides an integral representation
            for the likelihood risk from which several results follow, including,
            importantly, an explicit representation for the gradient.
            Proposition <ref> relies on defining a suitable space of
            log-hazard functions defined on the time-covariate domain
            [0,1]×.  Identifying this space is the key insight that
            allows us to rescue the likelihood approach and to derive the gradient
            needed to implement gradient boosting. Arriving at this framework is
            not conceptually trivial, and may explain the absence of boosted
            nonparametric hazard estimators until now.
            
            Algorithm <ref> of Section <ref> describes our estimator. The algorithm minimizes the likelihood risk (<ref>) over the defined space of log-hazard functions.
            In the special case of regression tree learners, expressions for the likelihood risk and its gradient are obtained from Proposition <ref>, which are then used to describe a tree-based implementation of our estimator in Section <ref>. In Section <ref> we apply it to a high-dimensional dataset generated from a naturalistic simulation of patient service times in an emergency
            department.
            
            Section <ref> establishes the consistency of the
            procedure. We show that the hazard estimator is consistent if the
            space is correctly specified.  In particular, if the space is the span
            of regression trees, then the hazard estimator satisfies an oracle
            inequality and recovers ł up to some error tolerance (Propositions <ref> and <ref>).
            
            Another contribution of our work is to clarify the mechanisms used by gradient boosting to avoid overfitting.  Gradient boosting typically applies two types of regularization to invoke slow learning: (i) A small
            step-size is used for the update; and (ii) The number of boosting
            iterations is capped.  The number of iterations used in our algorithm
            is set using the framework of Zhang and Yu <cit.>, whose work shows how stopping early ensures consistency. On the other hand, the role of step-size restriction is more mysterious. While <cit.> demonstrates that small step-sizes are needed to prove consistency,
            unrestricted greedy step-sizes are already small enough for
            classification problems <cit.> and also for commonly used
            regression losses (see the Appendix of <cit.>). We show in
            Section <ref> that shrinkage acts as a counterweight to the curvature of the risk (see Lemma <ref>). Hence if the curvature is unbounded, as is the case for hazard regression, then the step-sizes may need to be explicitly controlled to ensure convergence.  This important result adds to our understanding of statistical convergence in gradient boosting.  As noted by Biau and Cadre <cit.> the literature for this is relatively sparse, which motivated them to propose another regularization mechanism that also prevents overfitting. 
            
            
            Concluding remarks can be found in Section <ref>. Proofs not appearing in the body of the paper can be found in the Appendix.
            
            
            
            The boosted hazard estimator
            In this section, we describe our boosted hazard estimator. To provide readers with concrete examples for the ideas introduced here, we will show how the quantities defined in this section specialize in the case of regression trees, which is one of a few possible ways to implement boosting.
            
            We begin by defining in Section <ref> an appropriate
            sample-dependent domain for the likelihood risk _n(F).  As
            explained, this key insight allows us to re-express the likelihood risk and its population analogue as smooth convex functionals, thereby enabling us to compute their gradients in closed form in
            Propositions <ref> and <ref> of
            Section <ref>. Following this, the boosting
            algorithm is formally stated in Section <ref>.
            
            
            Specifying a domain for _n(F)
            We will make use of two identifiability conditions <ref> and <ref> to define the domain for _n(F). Condition <ref> below is the same as Condition 1(iv)
            of Huang and Stone <cit.>.
            
            
            	The true hazard function
            	λ(t,x) is bounded between some interval
            	[Λ_L,Λ_U]⊂(0,∞) on the time-covariate
            	domain [0,1]×.
            
            
            Recall that we defined X(·) and Y(·) to be predictable processes, and so it
            can be shown that the integrals and expectations appearing in this
            paper are all well defined. Denoting the indicator function as
            I(·), define the following population and empirical
            sub-probability measures on [0,1]×:
            μ(B)
               =   ( ∫_0^1 Y(t)· I[{t,X(t)}∈ B] dt ),
            (B)
               =   1/n∑_i=1^n∫_0^1 Y_i(t)· I[{t,X_i(t)}∈ B] dt,
            
            and note that (B)=μ(B) because the data is
            i.i.d. by assumption.
            Intuitively,  measures the denseness of the observed sample
            time-covariate paths on [0,1]×. For any integrable f,
            ∫ f dμ   =   ( ∫_0^1 Y(t)· f(t,X(t)) dt ),
            ∫ f d   =   1/n∑_i=1^n∫_0^1 Y_i(t) · f(t,X_i(t)) dt.
            This allows us to define the following (random) norms and inner products
            f_,1    =    ∫|f| d
            f_,2    =    (∫ f^2 d)^1/2
            f_∞    =    sup{ |f(t,x)|:(t,x)∈[0,1]×}
            ⟨ f_1,f_2⟩ _    =    ∫ f_1f_2 d,
            
            and note that
            ·_,1≤·_,2≤·_∞
            because ([0,1]×)≤1.
            
            By careful design,  allows us to specify a natural domain for
            _n(F). Let {ϕ_j(t,x)}_j=1^d be a set of bounded
            functions [0,1]×↦[-1,1] that are linearly independent,
            in the sense that
            ∫_[0,1]×(∑_jc_jϕ_j)^2dtdx=0 if and only
            if c_1=⋯=c_d=0 (when some of the covariates are
            discrete-valued, dx should be interpreted as the product of
            a counting measure and the Lebesgue measure). The span of the functions is
            
                ={∑_j=1^dc_jϕ_j:c_j∈} .
            
            For example, the span of all regression tree functions that can be defined on [0,1]× is ={∑_j c_jI_B_j(t,x) : c_j ∈},[It is clear that said span is contained in . For the converse, it suffices to show that  is also contained in the span of trees of some depth. This is easy to show for trees with p+1 splits, because they can generate partitions of the form (-∞,t]×(-∞,x^(1)]×⋯×(-∞,x^(p)] in [0,1]× (Section 3 of <cit.>).] which are linear combinations of indicator functions over disjoint time-covariate cubes indexed[With a slight abuse of notation, the index j is only considered multi-dimensional when describing the geometry of B_j, such as in (<ref>). In all other situations j should be interpreted as a scalar index.] by j=(j_0,j_1,⋯,j_p):
            
                B_j={[                                                                                                             (t,x) ∈ [0,1]×                                                                                                                          : [         t^(j_0)<t≤ t^(j_0+1); x^(1,j_1)<x^(1)≤ x^(1,j_1+1);                            ⋮; x^(p,j_p)<x^(p)≤ x^(p,j_p+1) ] ]}.
            
            The regions B_j are formed using all possible split points
            {x^(k,j_k)}_j_k for the k-th coordinate x^(k), with the spacing determined by the precision of the measurements. For example, if weight is measured to the closest kilogram, then the set of all possible split points will be {0.5, 1.5, 2.5,⋯} kilograms. Note that these split points are the finest possible for any realization of weight that is measured to the nearest kilogram. While abstract treatments of trees assume that there is a continuum of split points, in reality they fall on a discrete (but fine) grid that is pre-determined by the precision of the data.
            
            
            
            When  is equipped with ⟨·,·⟩
            _, we obtain the following sample-dependent subspace of
            L^2(), which is the appropriate domain for
            _n(F):
            
                (,⟨·,·⟩_).
            
            Note that the elements in (,⟨·,·⟩_) are equivalence classes rather than actual functions
            that have well defined values at each (t,x). This is a problem
            because the likelihood risk (<ref>) requires evaluating F(t,x)
            at the points (T_i,X_i(T_i)) where Δ_i=1. We resolve this by fixing an orthonormal basis {_nj(t,x)}_j
            for (,⟨·,·⟩ _), and represent each member of (,⟨·,·⟩ _) uniquely in the form
            ∑_jc_j_nj(t,x). For example in the case of regression trees, applying the Gram-Schmidt procedure to {ϕ_j(t,x) = I_B_j(t,x)}_j gives
            
                {_nj(t,x) }_j =
                {I_B_j(t,x)/(B_j)^1/2: (B_j)>0},
            
            which by design have disjoint support.
            
            The second condition we impose is for {ϕ_j}_j=1^d to be
            linearly independent in L^2(μ), that is
            ∑_jc_jϕ_j_μ,2^2 = ∑_ijc_i(∫ϕ_iϕ_jdμ)c_j=0
            if and only if c_1=⋯=c_d=0. Since by construction
            {ϕ_j}_j=1^d are already linearly independent in
            [0,1]×, the condition intuitively requires the set of all
            possible time-covariate trajectories to be adequately dense in
            [0,1]× to intersect a sufficient amount of the support of
            every ϕ_j.  This is weaker than the identifiability conditions 1(ii)-1(iii) in
            <cit.> which require X(t) to have a positive joint probability
            density on [0,1]×.
            
            
            	The Gram matrix Σ_ij=∫ϕ_iϕ_jdμ is
            	positive definite.
            Integral representations for the likelihood risk
            Having deduced the appropriate domain for _n(F), we can now
            recast the risk as a smooth convex functional on (,⟨·,·⟩ _). Proposition <ref> below provides closed form expressions for this and its gradient. We note that if the risk is actually of a certain simpler form, it might be possible to estimate its gradient empirically from our risk expression using <cit.>.
            
            
            For functions F(t,x), f(t,x) of the form ∑_j c_j _nj(t,x), the likelihood risk (<ref>) can be written as
            _n(F)=∫(e^F- F)d,
            
            where ∈(,⟨·,·⟩_)
            is the function
            
                (t,x)
                =1/n∑_j{∑_i=1^nΔ_i_nj(T_i,X_i(T_i))}_nj(t,x).
            
            Thus there exists ρ∈(0,1) (depending on F and f)
            for which the Taylor representation
            
                _n(F+f)=_n(F)+
                	⟨ g_F,f⟩_+1/2∫
                	e^F+ρ ff^2d
            
            holds, where the gradient
            
                g_F(t,x)
                	=∑_j⟨ e^F,_nj⟩
                	__nj(t,x)
                	-(t,x)
            
            of _n(F) is the projection of e^F- onto
            (,⟨·,·⟩ _).
            Hence if g_F=0 then the infimum of _n(F) over
            the span of {_nj(t,x)}_j is uniquely attained
            at F.
            
            
            For regression trees the expressions (<ref>) and
            (<ref>) simplify further because  is closed under
            pointwise exponentiation, i.e. e^F∈ for F∈. This is
            because the B_j's are disjoint so F = ∑_j c_j I_B_j and hence
            e^F = ∑_j e^c_j I_B_j. Thus
            (t,x)    =   ∑_j:(B_j)>0_j/n(B_j) I_B_j(t,x),
            
            _n(F)    =   ∑_j:(B_j)>0(e^c_j(B_j)
            -c_j_j/n),
            
            
            g_F(t,x)    =   ∑_j:(B_j)>0( e^c_j - _j/n(B_j)) I_B_j(t,x),
            
            where
            
                _j=∑_iΔ_iI[{T_i,X_i(T_i)}∈ B_j]
            
            is the number of observed failures in the time-covariate region
            B_j.
            
             Fix a realization of {(X_i(·),Y_i(·),T_i)}_i=1^n.
            	Using (<ref>) we can rewrite (<ref>) as
            	
                _n(F)
                	=∫ e^Fd-1/n∑_i=1^nΔ_iF(T_i,X_i(T_i)).
            
            	We can express F in terms of the basis {_nk}_k as
            	F(t,x)=∑_kc_k_nk(t,x).  Hence
            	∫ F d   =   ∫1/n∑_j{∑_i=1^nΔ_i_nj(T_i,X_i(T_i))}_nj(t,x) F(t,x) d
               =   1/n∑_j{∑_i=1^nΔ_i_nj(T_i,X_i(T_i))}∫_nj(t,x) F(t,x) d
               =   1/n∑_j{∑_i=1^nΔ_i_nj(T_i,X_i(T_i))}∫_nj(t,x) ∑_kc_k_nk(t,x) d
               =   1/n∑_j{∑_i=1^nΔ_i_nj(T_i,X_i(T_i))}
            	c_j
               =   1/n∑_i=1^nΔ_i∑_jc_j_nj(T_i,X_i(T_i))
               =   1/n∑_i=1^nΔ_iF(T_i,X_i(T_i)),
            	
            	where the fourth equality follows from the orthonormality
            	of the basis. This completes the derivation of (<ref>).
            	
            	By an interchange argument we obtain
            	d/dθ_n(F+θ f)
            	   =   ∫(e^F+θ f-)fd,
            d^2/dθ^2_n(F+θ f)
            	   =   ∫ e^F+θ ff^2d,
            	
            	the latter being positive whenever f≠0; i.e., _n(F) is
            	convex. The Taylor representation (<ref>) then follows
            	from noting that g_F is the orthogonal projection of
            	e^F-∈ L^2() onto (,⟨·,·⟩ _).
            
            
            The expectation of the likelihood risk also has an integral
            representation. A special case of the
            representation (<ref>) below is proved in Proposition 3.2
            of <cit.> for right-censored data only, under assumptions that do not allow for internal covariates. In the statement of the proposition below recall that
            Λ_L and Λ_U are defined in <ref>.  The constant
             is defined later in (<ref>).
            
            
            For F∈∪{logλ},
            
                R(F)={_n(F)}=∫(e^F-λ F)dμ.
            
            Furthermore the restriction of R(F) to 
            is coercive:
            
                1/2R(F)≥Λ_L/F_∞
                	+Λ_Umin{0,1-log(2Λ_U)},
            
            and it attains its minimum at a unique point F^*∈(,⟨·,·⟩_μ). If 
            contains the underlying log-hazard function then F^*=logλ.
            
            
            
            
            Coerciveness (<ref>) implies that any F with expected
            risk R(F) less than R(0)≤1<3 is uniformly bounded:
            F_∞
            </Λ_L[3/2+Λ_Umax{0,log(2Λ_U)-1}]
            ≤β_Λ
            where the constant
            
                β_Λ=3/2+Λ_Umax{0,log(2Λ_U)-1}/min{1,Λ_L}
            
            is by design no smaller than 1 in order to simplify subsequent analyses.
            The boosting procedure
            In gradient boosting the key idea is to update an iterate in a direction that is approximately aligned to the negative gradient. To model this direction formally, we introduce the concept of an -gradient.
            
            
            
            Suppose g_F 0. We say that a unit vector g_F^∈(,⟨·,·⟩ _) is an -gradient at F if
            for some 0<≤1,
            ⟨g_F/g_F_,2,g_F^⟩_≥.
            
            Call -g_F^ a negative -gradient if g_F^ is an -gradient.
            
            Our boosting procedure seeks approximations g_F^ that
            satisfy eq:epsgradient for some pre-specified alignment value .
            The larger  is, the closer the alignment is between the negative gradient
            and the negative -gradient, and the greater the risk reduction. In
            particular, -g_F is the unique negative 1-gradient with maximal risk reduction.  In practice, however, we find that using a smaller value of  leads to simpler approximations that prevent overfitting in finite samples. This is consistent with other implementations of boosting: It is well known that the statistical performance of gradient descent generally improves when simpler base learners are used.
            
            
            Algorithm <ref> describes the proposed boosting procedure for
            estimating λ. For a given level of alignment , Line 3
            finds an -gradient g__m^ at _m satisfying (<ref>) at the m-th iteration, and uses its negation for the boosting update in Line 4. If the  -gradients are tree learners, as is the case with the implementation in Section <ref>, then the trees cannot be grown in the same way as the standard boosting algorithm in Friedman <cit.>. This is because the standard approach grows all regression trees to a fixed depth, which may or may not ensure -alignment at each boosting iteration.
            
            To ensure -alignment, the depth of the trees are not fixed in the implementation in Section <ref>. Instead, at each boosting iteration a tree is grown to whatever depth is needed to satisfy eq:epsgradient.alg. This can always be done because the alignment  is non-decreasing in the number of tree splits, and with enough splits we can recover the gradient g__m itself up to -almost everywhere.[Split the tree until each leaf node contains just one of the regions B_j in (<ref>) with (B_j)>0. Then set the value of the node equal to the value of the gradient function (<ref>) inside B_j.]  As mentioned earlier, we recommend using small values of , which can be determined in practice using cross-validation. This differs from the standard approach where cross-validation is used to select a common tree depth to use for all boosting iterations.
            
            
            10pt
            
            
            
            In addition to the gradient alignment , Algorithm <ref> makes use of two other regularization parameters, Ψ_n and ν_n. The first defines the early stopping criterion (how many boosting iterations to use), while the second controls the step-sizes of the boosting updates. These are two common regularization techniques used in boosting:
            
              * Early stopping. The number of boosting iterations
             is controlled by stopping the algorithm before the uniform
            norm of the estimator F__∞ reaches or
            exceeds
            Ψ_n=W(n^1/4)→∞,
            
            where W(y) is the branch of the Lambert function that returns the
            real root of the equation ze^z=y for y>0.
            
              * Step-sizes.  The step-size ν_n ≪ 1 used in gradient
            boosting is typically held constant across iterations. While we can also do this in our procedure,[The term ν_n^2 e^Ψ_n in condition (<ref>) would need to be replaced by ν_n^2 e^Ψ_n if a constant step-size is used.] the role of step-size shrinkage becomes more salient if we use ν_n/(m+1) instead as the step-size for the m-th iteration in Algorithm <ref>. This step-size is controlled in two ways. First, it is made to decrease with each iteration according to the Robbins-Monro condition that the sum of the steps diverges while the sum of squared steps converges. Second, the shrinkage factor ν_n is selected to make the step-sizes decay with n at rate
            ν_n^2e^Ψ_n<1,
            15pt ν_n^2e^Ψ_n→ 0.
            
            This acts as a counterbalance to _n(F)'s unbounded curvature:
            d^2/dθ^2_n(F+θ
            f)|_θ=0=∫ e^F f^2d,
            
            which is upper bounded by e^Ψ_n when F_∞<Ψ_n
            and f_,2=1.
            Consistency
            Under <ref> and <ref>, guarantees for our hazard estimator  in Algorithm <ref> can be derived for two scenarios of interest. The guarantees rely on the regularizations described in Section <ref> to avoid overfitting. In the following development, recall from Proposition <ref> that F^* is the unique minimizer of R(F), so it satisfies the first order condition
            ⟨ e^F^*-λ,F⟩ _μ=0
            
            for all F∈. Recall that the span of all trees is closed under pointwise exponentiation (e^F∈), in which case (<ref>) implies that ł^* = e^F^* is the orthogonal projection of λ onto (,⟨·,·⟩_μ).
            
            
              *  Consistency when  is correctly specified.
            If the true log-hazard function logλ is in ,
            then Proposition <ref> asserts that
            F^*=logλ. It will be shown in this case
            that  is consistent:
            
                ‖-λ‖_μ,2^2 = o_p(1).
            
              *  Oracle inequality for regression trees. If  is closed under pointwise exponentiation, it follows from (<ref>) that λ^* is the best L^2(μ)-approximation to λ among all candidate hazard estimators {e^F:F∈}. It can then be shown that  converges to this best approximation:
            
                ‖ -λ‖ _μ,2^2
                =‖ł^* -λ‖ _μ,2^2
                + o_p(1).
            
            This oracle result is in the spirit of the type of guarantees available for tree-based boosting in the non-functional data setting. For example, if tree stumps are used for L_2-regression, then the regression function estimate will converge to the best approximation to the true regression function in the span of tree stumps <cit.>. Similar results also exist for boosted classifiers <cit.>.
            
            
            Propositions <ref> and <ref> below
            formalize these guarantees by providing bounds on the error terms above. While sharper bounds may exist, the purpose of this paper is to introduce our generic estimator for the first time and to provide guarantees that apply across different implementations. More refined convergence rates may exist for a specific implementation, just like the analysis in Bühlmann and Yu <cit.> for L_2Boosting when componentwise spline learners are specifically used.
            
            En route to establishing the guarantees, Lemma <ref> below clarifies the role played by step-size restriction in ensuring convergence of the estimator. As explained in the Introduction, explicit shrinkage is not necessary for classification and regression problems where the risk has bounded curvature. Lemma <ref> suggests that it may, however, be needed when the risk has unbounded curvature, as is the case with _n(F). Seen in this light, shrinkage is really a mechanism for controlling the growth of the risk curvature.
            
            Strategy for establishing guarantees
            The representations for _n(F) and its population analogue R(F)
            from Section <ref> are the key ingredients for formalizing the guarantees. We use them to first show that _∈(,⟨·,·⟩_) converges to F^*∈(,⟨·,·⟩_μ): Applying Taylor's theorem to the representation for R(F) in Proposition <ref> yields
            
                ‖_-F^*‖_μ,2^2≤
                	2R(_)-R(F^*)/min_t,x(ł^*∧).
            
            The problem is thus transformed into one of risk minimization
            R(_)→ R(F^*), for which <cit.>
            suggests analyzing separately the terms of the decomposition
            
            0
               ≤   
            R(_)-R(F^*)
            
               ≤   |_n(_)-R(_)|
            40pt(I)     complexity argument
            
               +   |_n(F^*)-R(F^*)|
            40pt(II)     standard argument
            
               +   {_n(_)-_n(F^*)}.
            25pt(III)     curvature argument
            
            The authors argue that in boosting, the point of limiting the number
            of iterations  (enforced by lines 5-10 in Algorithm <ref>) is to prevent _ from growing too fast, so that (I) converges to zero as n→∞. At the same time,  is allowed to grow with n in a controlled manner so that the empirical risk _n(_) in (III) is eventually minimized as n→∞. Lemmas <ref> and <ref> below show that our procedure achieves both
            goals.  Lemma <ref> makes use of complexity theory
            via empirical processes, while Lemma <ref> deals with
            the curvature of the likelihood risk.  The term (II) will be bounded
            using standard concentration results.
            
            
            Bounding (I) using complexity
            To capture the effect of using a simple negative -gradient (<ref>) as the descent direction, we bound (I) in terms of the complexity of[
            	For technical convenience, _ has been enlarged from
            	_,boost to include the unit ball.]      _ = _,∪{F∈:F_∞=1}⊆,
            
                      where _,={_m
            =-∑_k=0^m-1ν_n/k+1g__k^:m = 0,1,…}.
            
            Depending on the choice of weak learners for the -gradients,
            _ may be much smaller than .  For example, coordinate
            descent might only ever select a small subset of basis functions
            {ϕ_j}_j because of sparsity. As another example if ł(t,x) is additively separable in time and also in each covariate, then regression trees might only ever select simple tree stumps (one tree split).
            
            The measure of complexity we use below comes from empirical process
            theory.  Define _^Ψ={F∈_:F_∞<Ψ}
            for Ψ>0 and suppose that Q is a sub-probability measure on
            [0,1]×.  Then the L^2(Q)-ball of radius δ>0
            centred at some F∈ L^2(Q) is
            {F'∈_^Ψ:F'-F_Q,2<δ}.  The covering
            number (δ,_^Ψ,Q) is the minimum number of
            such balls needed to cover _^Ψ (Definitions 2.1.5 and
            2.2.3 of van der Vaart and Wellner <cit.>), so (δ,_^Ψ,Q)=1 for
            δ≥Ψ. A complexity measure for _ is
            
                J__=sup_Ψ,Q{∫_0^1{log(uΨ,_^Ψ,Q)}^1/2 du},
            
            where the supremum is taken over Ψ>0 and over all non-zero
            sub-probability measures. As discussed, J__ is never greater than, and potentially much smaller than J_, the complexity of , which is fixed and finite.
            
            
            Before stating Lemma <ref>, we note that the result
            also shows an empirical analogue to the norm equivalences
            
                F_μ,1≤F_μ,2≤F_∞≤/2F_μ,1 for all  F∈
            
            exists, where
            
                =2sup_F∈:F_∞=1(F_∞/F_μ,1)
                	=2/inf_F∈:F_∞=1F_μ,1>1.
            
            The factor of 2 serves to simplify the presentation, and can be
            replaced with anything greater than 1.
            
            6pt
            
            There exists a universal constant
            κ such that for any 0<η<1, with probability at least
            
                1-4exp{ -(η n^1/4/κ J__)^2}
            
            an empirical analogue to (<ref>) holds for all F∈:
            
                F_,1≤F_,2≤F_∞≤F_,1,
            
            and for all F∈_^Ψ_n,
            
                |{_n(F)-_n(0)}-{R(F)-R(0)}|<η.
            6pt
            
            
            The equivalences (<ref>) imply that
            (,⟨·,·⟩ _)
            equals its upper bound =d. That is, if
            ∑_jc_jϕ_j_,2=0, then
            ∑_jc_jϕ_j_∞=0, so c_1=⋯=c_d=0
            because {ϕ_j}_j=1^d are linearly independent on
            [0,1]×.
            Bounding (III) using curvature
            We use the representation in Proposition <ref> to study the
            minimization of the empirical risk _n(F) by boosting. Standard
            results for exact gradient descent like Theorem 2.1.15
            of Nesterov <cit.> are in terms of the norm of the minimizer, which
            may not exist for _n(F).[The infimum of _n(F) is not always attainable: If f is non-positive and vanishes on the set {{T_i,X_i(T_i)}:Δ_i=1}, then _n(F+θ f)=∫(e^F+θ f- F)d is decreasing in θ so f is a direction of recession. This is however not an issue for boosting because of early stopping.]  If coordinate descent is used instead, Section 4.1 of <cit.> can be applied to convex functions whose infimum may not be attainable, but its curvature is
            required to be uniformly bounded above. Since the second derivative of
            _n(F) is unbounded (<ref>),
            Lemma <ref> below provides two remedies: (i) Use the
            shrinkage decay (<ref>) of ν_n to counterbalance
            the curvature; (ii) Use coercivity (<ref>) to show that
            with increasing probability, {_m}_m=0^ are uniformly
            bounded, so the curvatures at those points are also uniformly
            bounded. Lemma <ref> combines both to derive a result
            that is simpler than what can be achieved from either one alone. In
            doing so, the role played by step-size restriction becomes clear. The
            lemma relies in part on adapting the analysis in Lemma 4.1
            of <cit.> for coordinate descent to the case for generic
            -gradients. The conditions required below will be shown to hold
            with high probability.
            
            
            Suppose (<ref>) holds and that
            
                |_n(F^*)-R(F^*)|<1,
                15pt
                sup_F∈_^Ψ_n|_n(F)-R(F)|<1.
            
            Then the largest gap between F^* and {_m}_m=0^,
            
                γ̂=max_m≤m̂_m-F^*_∞∨1,
            
            is bounded by a constant no greater than 2β_Λ,
            and for n≥55,
            
                _n(_)-_n(F^*)
                	<2eβ_Λ(log
                		n/4n^1/4)^/(γ̂)
                	+ν_n^2e^Ψ_n.
            
            The last term in (<ref>) suggests that the role of the
            step-size shrinkage ν_n is to keep the curvature of the
            risk in check, to prevent it from derailing convergence. Recall
            from (<ref>) that e^Ψ_n describes the
            curvature of _n(_m). Thus our result clarifies the
            role of step-size restriction in boosting functional data.
            
            
            
            
            
            
            
            Regardless of whether the risk curvature is bounded or not, smaller
            step-sizes always improve the convergence bound. This can be seen from
            the parsimonious relationship between ν_n
            and (<ref>).  Fixing n, pushing the value of
            ν_n down towards zero yields the lower limit
            
                2eβ_Λ(log n/4n^1/4)^/(γ̂).
            
            However, this limit is unattainable as ν_n must be positive in
            order to decrease the risk.
            This effect has been observed in practical applications of boosting.
            Friedman <cit.> noted improved performance for gradient boosting with
            the use of a small shrinkage factor ν.  At the same time, it was
            also noted there was diminishing performance gain as ν became
            very small, and this came at the expense of an increased number of
            boosting iterations.  This same
            phenomenon has also been observed for
            L_2Boosting <cit.> with componentwise linear learners.
            It is known that the solution path for L_2Boosting closely matches
            that of lasso as ν→ 0. However, the algorithm exhibits
            cycling behaviour for small ν, which greatly increases the number
            of iterations and offsets the performance gain in trying to
            approximate the lasso (see Ehrlinger and Ishwaran <cit.>).
            Formal statements of guarantees
            As a reminder, we have defined the following quantities:
                =    e^_, the boosted hazard estimator in Algorithm <ref>
            ł^*            =     e^F^*,  where F^* is the unique minimizer of R(F) in Proposition <ref>
            Λ_L,Λ_U     =    lower and upper bounds on λ(t,x) as defined in <ref>
            γ̂    =    maximum gap between F^* and
            		{_m}_m=0^ defined in (<ref>)
             =    a universal constant
                =    constant defined in (<ref>)
            β_Ł    =    constant defined in (<ref>)
            
            	J__    =    complexity measure (<ref>), bounded above by  J_
            To simplify the results, we will assume that n≥55 and also set
            the shrinkage to satisfy ν_n^2e^Ψ_n=log
            n/(64n^1/4). Our first guarantee shows that our hazard estimator is consistent if the model is correctly specified.
            
             (Consistency under correct model specification). Suppose  contains the true log-hazard function logλ. Then with probability
            
                1-8exp{ -(log n/κ (Λ_L^-1∨Λ_U)
                	J__)^2}
            
            we have that _∞ is bounded and
            
                ‖-ł‖ _μ,2^2
                <13β_Łmax_t,x(ł∨)^2/min_t,x(ł∧)(log n/4n^1/4)^/(γ̂).
            
            Thus  is consistent.
            
            
            Via the tension between  and J__, Proposition <ref> captures the trade-off in statistical performance between weak and strong learners in gradient boosting.  The advantage of low complexity (weak learners) is reflected in the increased probability of the L^2(μ)-bound holding, with this probability being maximized when J__→0, which generally occurs as → 0. However, diametrically opposed to this, we find that the L^2(μ)-bound is minimized by → 1, which occurs with the use of stronger learners that are more aligned with the gradient.  This same trade-off is also captured by our second guarantee which establishes an oracle inequality for tree learners.
            
             (Oracle inequality for tree learners). Suppose e^F∈ for F∈. Then among {e^F:F∈}, ł^* is the best L^2(μ)-approximation to ł, that is
            
                ł^* = min_e^F:F∈‖ e^F-λ‖ _μ,2.
            
            Moreover,  converges to this best approximation ł^*: With probability
            
                1-8exp{-(log n/κ(Λ_L^-1∨Λ_U) J__)^2}
            
            we have that _∞ is bounded and
            
                ‖-λ‖ _μ,2^2 <
                ρ_^2 +
                13β_Λmax_t,x(Λ_U∨)^2/min_t,x(Λ_L∧)(log n/4n^1/4)^/(γ̂),
            
            where ρ_^2 = ł^*-ł_μ,2^2 is the smallest error one can achieve from using functions in {e^F:F∈} to approximate λ.
            
            
            For tree learners, ł^*(t,x) is constant over each region B_j in (<ref>), and its value equals the local average of ł over B_j,
            
                ł^*(t,x)|_B_j = 1/μ(B_j)∫_B_jλ dμ.
            
            Hence if the B_j's are small, ł^* should closely approximate ł (recall from Remark <ref> that the size of the B_j's is fixed by the data). To estimate the approximation error ρ_ in terms of B_j, suppose that λ is sufficiently smooth, e.g. Hölder continuous |λ(t,x)-λ(t',x')| ≾(t-t',x-x')^b for some b>0. Then since inf_B_jλ≤ł^*|_B_j≤sup_B_jλ,
            
                ρ_≤λ^*-λ_∞≾max_j(B_j)^b.
            A tree-based implementation
            Here we describe an implementation of Algorithm <ref> using
            regression trees, whereby the -gradient g__m^ is obtained by growing a tree to satisfy eq:epsgradient.alg for a pre-specified .
            
            To explain the tree growing process, first observe that the m-th step log-hazard estimator is an additive expansion of CART basis functions. Thus it can be written as
            _m(t,x)    =   ∑_b=0^m-1∑_l=1^L_bγ_b,l I_A_b,l(t,x) 
            
               =   ∑_j c_m,j I_B_j(t,x), 
            where A_b,l is the l-th leaf region of the b-th tree. Recall from Section <ref> that each tree is grown until eq:epsgradient.alg is satisfied, so the number of leaf nodes L_b can vary from tree to tree. The leaf regions are typically large subsets of the time-covariate space [0,1]× adaptively determined by the tree growing process (to be discussed shortly). Since each leaf region can be further decomposed into the finer disjoint regions B_j in eq:regions, _m(t,x) can be rewritten as (<ref>). However, many of these regions will share the same coefficient value, so (<ref>) can be written more compactly as
            
                _m(t,x) = ∑_j c_m,j I_B_m,j'(t,x),
            
            where B_m,j' is the union of contiguous regions whose coefficient equals c_m,j. This smooths the hazard estimator (t,x) over [0,1]×, thanks to the regularization imposed by limiting the number of trees (early stopping) and also by the use of weak tree learners.  This is unlike the unconstrained hazard MLE (t,x) defined in eq:mle.hazard, which can take on a different value in each region B_j, making it prone to overfit the data.
            
            To construct an -gradient g__m^ with -alignment to g__m defined by (<ref>),
            
                g__m(t,x) =
                ∑_j:(B_j)>0( e^c_m,j -
                _j/n(B_j)) I_B_j(t,x),
            
            the tree splits are adaptively chosen to reduce the L^2()-approximation error between g__m^ and g__m. We implement tree splits for both time and covariates. Specifically, suppose we wish to split a leaf region A ⊆ [0,1]× into left and right daughter subregions A_1 and A_2, and assign values γ_1 and γ_2 to them.  For example, a split on the k-th covariate could propose left and right daughters such as
            
            A_1={(t,x)∈ A: x^(k)≤ s},
            A_2={(t,x)∈ A: x^(k)> s},
            
            or a split on time t could propose regions
            
            A_1={(t,x)∈ A: t≤ s},
            A_2={(t,x)∈ A: t> s}.
            
            Now note that g__m is constant within each region B_j. We denote its value by g__m(t_B_j,x_B_j) where (t_B_j,x_B_j) is the centre of B_j. Hence the best split of A into A_1 and A_2 is the one that minimizes
            
                -35pt
                	min_γ_1∫_A_1{ g__m(t,x)
                	- γ_1 }^2 d + min_γ_2∫_A_2{g__m(t,x)
                	- γ_2 }^2 d
                   =
                	min_γ_1∑_j:B_j⊆ A_1(B_j) ·{ g__m(t_B_j,x_B_j) - γ_1 }^2 
                       +
                	min_γ_2∑_k:B_k⊆ A_2(B_k)·{ g__m(t_B_k,x_B_k) - γ_2 }^2 
                   =
                	min_γ_1∑_j:z_j∈ A_1,
                w_j>0 w_j· (ỹ_j - γ_1)^2
                	+ min_γ_2∑_k:z_k∈ A_2,
                w_k>0 w_k· (ỹ_k - γ_2)^2,
            
            where
            
                ỹ_j = g__m(t_B_j,x_B_j)
                = e^c_m,j - _j/n(B_j)
            
            represents the j-th pseudo-response, z_j = (t_B_j,x_B_j)
            its covariate and w_j = (B_j) its weight. Thus the splits use a weighted least squares criterion, which can be efficiently computed as usual.
            
            
            
            
            
            
            
            
            
            
            
            We split the tree until eq:epsgradient.alg is satisfied, resulting in L_m leaf nodes (L_m-1 splits).  As discussed in Section <ref>, we can always find a deep enough tree that is an -gradient because with enough splits we can recover the gradient g__m itself. Recall also that a small value of  performs best in practice, and this can be chosen by cross-validating on a set of small-sized candidates: For each one we implement Algorithm <ref>, and we select the one that minimizes the cross-validated risk _n(F) defined in (<ref>). By contrast, the standard boosting algorithm <cit.> uses cross-validation to select a common number of splits to use for all trees, which does not ensure that each tree is an -gradient.
            
            Regarding the possible split points for the covariates covariate.split.point, note that the k-th covariate x^(k) = x^(k)(t) is a time series that is sampled periodically. This yields a set of unique values equal to the union of all of the sampled values for the n observations. In direct analogy to non-functional data boosting, we place candidate split points in-between the sorted values in this set.  In other words, splits for covariates only occur at values corresponding to the observed data just as in non-functional boosting.
            
            The resolution for the grid of candidate time
            splits time.split.point is set equal to the temporal
            resolution. For example, the covariate trajectories in the simulation
            in Section <ref> are piecewise constant and may change every
            0.002 days. Placing the candidate split points at {0.002,0.004,…} days simplifies the exact computation of (B_j)
            because every covariate trajectory is constant between these points.
            Again, notice that the splits for time only occur at values informed by
            the observed data.
            
            Putting it together, the setup above leverages our insight in
            (<ref>) by transforming the survival functional data
            into the data values {w_j, ỹ_j, z_j}_j:w_j>0, which
            enables the implementation to proceed like standard gradient boosting
            for non-functional data. Only the pseudo-response ỹ_j in
            {w_j, ỹ_j, z_j} needs to be updated at each boosting
            iteration, while the other two do not change. In terms of storage it
            costs 𝒪(np|𝒯|) to store {w_j, ỹ_j,
            z_j}_j:w_j>0, where |𝒯| is the cardinality of the set
            of candidate time splits.[Each {w_j, ỹ_j, z_j} is
            	of dimension p+3 and the number of time-covariate regions B_j
            	with w_j>0 is at most n(|𝒯|+1). To show the latter,
            	observe that B_j will only have w_j=(B_j)>0 if it is
            	traversed by at least one sample covariate trajectory. Then note
            	that each of the n sample covariate trajectories can traverse at
            	most |𝒯|+1 unique regions.] Computationally, choosing a
            new tree split requires testing 𝒪(np|𝒯|)
            candidate splits.[A sample covariate trajectory can have at
            	most |𝒯| unique observed values for the k-th covariate
            	x^(k), so there are at most n|𝒯| candidate splits
            	for x^(k). Thus there are 𝒪(np|𝒯|)
            	candidate splits for p covariates. The number of candidate splits
            	on time is obviously |𝒯|.] The space and time
            complexities of the implementation are reasonable given that they are
            𝒪(np) for non-functional data boosting: In the functional
            data setting, each sample can have up to |𝒯| observations,
            so n functional data samples is akin to
            𝒪(n|𝒯|) samples in a non-functional data
            setting.
            
            
            
            Numerical experiment
            We now apply the boosting procedure of Section <ref> to a
            high-dimensional dataset generated from a naturalistic
            simulation. This allows us to compare the performance of our estimator
            to existing boosting methods. The simulation is of patient service
            times in an emergency department (ED), and the hazard function of
            interest is patient service rate in the ED. The study of patient
            transitions in an ED queue is an important one in healthcare
            operations, because without a high resolution model of patient flow
            dynamics, the ED may be suboptimally utilized which would deny patients
            of timely critical care.
            
            
            Service rate The service rate model used in the simulation is based upon a service time dataset from the ED of an academic hospital in the United States. The dataset contains information on 86,983 treatment encounters from 2014 to early 2015. Recorded for each encounter was: Age, gender, Emergency Severity Index (ESI)[Level 1 is the most severe
            	(e.g., cardiac arrest) and level 5 is the least (e.g., rash). We
            	removed level 1 patients from the dataset because they were treated
            	in a separate trauma bay.], time of day when treatment in the ED
            ward began, day of week of ED visit, and ward census. The last one
            represents the total number of occupied beds in the ED ward, which
            varies over the course of the patient's stay. Hence it is a
            time-dependent variable. Lastly, we also have the duration of the
            patient's stay (service time).
            
            The service rate function is developed from the data in the following
            way. First, we apply our nonparametric estimator to the data to
            perform exploratory analysis.  We find that:
            
              *  The key variables affecting the service rate (based on relative variable importance <cit.>) are ESI, age, and ward census. In addition, two of the most pronounced interaction terms identified by the tree splits are (≥34, =5) and (≥34, ≤4).
            
              *  Holding all the variables fixed, the shapes of the estimated service rate function resemble the hazard functions of log-normal distributions. This agrees with the queuing literature that find log-normality to be a reasonable parametric fit for service durations.
            
            Guided by these findings, we specify the service rate
            λ(t,X(t)) for the simulation as a log-normal accelerated
            failure time (AFT) model, and estimate its parameters from data.
            
            This yields the service rate
            
                λ(t,x) = θ(x) ·ϕ_l(θ(x)t;m,σ)/1-Φ_l(θ(x)t;m,σ),
            
            where ϕ_l(·;m,σ) and Φ_l(·;m,σ) are the
            PDF and CDF of the log-normal distribution with log-mean m=-1.8 and
            log-standard deviation σ=0.74. The function θ(x) captures
            the dependence of the service rate on the covariates:
            
                logθ(X(t))    = -0.0071· + 0.022· - min{ a·_t/70,2 }
                        + 0.10· I(≥ 34,=5) - 0.10· I(≥ 34,≤ 4) 
                        + 0·_1 + ⋯ + 0·_43.
            
            The specification for θ(X(t)) above is a slight modification of
            the original estimate, with the free parameter a allowing us to
            study the effect of time-dependent covariates on hazard
            estimation. When a=0, the service rate does not depend on
            time-varying covariates, but as a increases, the dependency becomes
            more and more significant. In the data, the ward census never exceeds
            70, so we set the capacity of the simulated ED to 70 as well. The
            min operator caps the impact that census can have on the simulated
            service rate as a grows. The irrelevant covariates _1,⋯,_43 are added to the data in order to assess how boosting performs in high dimensions. We explicitly include them in (<ref>) to remind ourselves that the simulated data is high-dimensional. Forty of the irrelevant variables are generated synthetically as described in the next subsection, while the rest are variables from the original dataset not used in the simulation.
            
            
            Simulation model
            Using (<ref>) and (<ref>), we simulate a naturalistic
            dataset of 10,000 patient visit histories. The value of a will be
            varied from 0 to 3 in order to study the impact of time-dependent
            covariates on hazard estimation. Each patient is associated with a
            46-dimensional covariate vector consisting of:
            
              *  The time-varying ward census. The initial value is sampled from its marginal empirical distribution in the original dataset. To simulate its trajectory over a patient's stay, for every timestep advance of 0.002 days (≈3 minutes), a Bernoulli(0.02) random variable is generated. If it is one, then the census is incremented by a normal random variable with zero mean and standard deviation 10. The result is truncated if it lies outside the range [1,70], the upper end being the capacity of the ED.
            
              *  The other five time-static covariates in the original dataset. These are sampled from their marginal empirical distributions in the original dataset. Two of the variables (age and ESI) influence the service rate, while the other three are irrelevant.
            
              *  An additional forty time-static covariates that do not affect the service rate (irrelevant covariates). Their values are drawn uniformly from [0,1].
            
            
            We also generate independent censoring times (rounded to the nearest 0.002 days) for each visit from an exponential distribution. For each simulation, the rate of the exponential distribution is set to achieve an approximate target of 25% censoring.
            
            
            Comparison benchmarks
            When the covariates are static in time, a few software packages are available for performing hazard estimation with tree ensembles. Given that the data is simulated from a log-normal hazard, we compare our nonparametric method to two correctly specified parametric estimators:
            
              *  The  blackboost estimator in the R package  mboost <cit.> provides a tree boosting procedure for fitting the log-normal hazard function. In order to apply this to the simulated data, we make ward census a time-static covariate by fixing it at its initial value.
            
              *  Transformation forests <cit.> in the R package  trtf can also fit log-normal hazards. Moreover, it allows for left-truncated and right-censored data. Since the ward census variable is simulated to be piecewise constant over time, we can treat each segment as a left-truncated and right-censored observation. Thus for this simulation, transformation forests are able to handle time-dependent covariates with time-static effects. This falls in between the static covariate/static effect  blackboost estimator and our fully nonparametric one.
            
            
            Since the service rate model used in the simulations is in fact log-normal, the benchmark methods above enjoy a significant advantage over our nonparametric one, which is not privy to the true distribution. In fact, when a=0 the log-normal hazard (<ref>) depends only on time-static covariates, so the benchmarks should outperform our nonparametric estimator. However, as a grows, we would expect a reversal in relative performance.
            
            To compare the performances of the estimators, we use Monte Carlo integration to evaluate the relative mean squared error
            
                %MSE = 𝔼_X[∫_0^1 {λ(t,X)-λ̂(t,X)}^2dt]/𝔼_X [∫_0^1 λ(t,X)^2dt].
            
            The Monte Carlo integrations are conducted using an independent test set of 10,000 uncensored patient visit histories. For the test set, ward census is held fixed over time at the initial value, and we use the grid {0,0.02,0.04,⋯,1} for the time integral. The nominator above is then estimated by the average of {λ(t,x)-λ̂(t,x)}^2 evaluated at the 51×10,000 points of (t,x). The denominator is estimated in the same manner.
            
            
            Results
            For the implementation of our estimator in Section <ref>, the value of  and the number of trees m̂ are jointly determined using ten-fold cross validation.  The candidate values we tried for  are {0.003, 0.004, 0.005, 0.006, 0.007}, and we limit m̂ to no more than 1,000 trees.  A wider range of values can be of course be explored for better performance  (at the cost of more computations). As comparison, we run an ad-hoc version of our algorithm in which all trees use the same number of splits, as is the case in standard boosting <cit.>.  This approach does not explicitly ensure that the trees will be -gradients for a pre-specified . The number of splits and the number of trees used in the ad-hoc method are jointly determined using ten-fold cross-validation.
            
            In order to speed up convergence at the m-th iteration for both approaches, instead of using the step-size ν_n/(m+1) of Algorithm <ref>, we performed line-search within the interval (0,ν_n/(m+1)].  While Lemma <ref> shows that a smaller shrinkage ν_n is always better, this comes at the expense of a larger m̂ and hence computation time. For simplicity we set ν_n=1 for all the experiments here.
            
            For fitting the  blackboost estimator, we use the default setting of  nu=0.1 for the step-size taken at each iteration. The other hyperparameters,  mstop (the number of trees) and  maxdepth (maximum depth of trees), are chosen to directly minimize the relative MSE on the test set. This of course gives the  blackboost estimator an unfair advantage over our estimator, which is on top of the fact that it is based on the same distribution as the true model. Transformation forest (using also the true distribution) is fit using code kindly provided by Professor T. Hothorn.[In the code 100 trees are used in the forest, which takes about 700 megabytes to store the fitted object when applied to our simulated data.]Variable selection.  The relative importance of
            variables <cit.> for our estimator are given in
            Table <ref> for all four cases a=0,1,2,3. The four
            factors that influence the service rate (<ref>) are
            explicitly listed, while the irrelevant covariates are grouped
            together in the last column. When a=0, the service rate does not
            depend on census, and we see that the importance of census and the
            other irrelevant covariates are at least an order of magnitude smaller
            than the relevant ones. As a increases, census becomes more and more
            important as correctly reflected in the table. Across all the cases
            the importance of the relevant covariates are at least an order of
            magnitude larger than the others, suggesting that our estimator is
            able to pick out the influential covariates and largely avoid the
            irrelevant ones.
            Presence of time-dependent covariates. Table <ref> presents the relative MSEs for the estimators as the service rate function (<ref>) becomes increasingly
            dependent on the time-varying census variable. When a=0 the service rate depends only on time-static covariates, so as expected, the parametric log-normal benchmarks perform the best when applied to data simulated from a log-normal AFT model.
            
            However, as a increases, the service rate becomes increasingly
            dependent on census. The corresponding performances of both benchmarks deteriorate dramatically, and is handily outperformed by the proposed estimator. We note that the inclusion of just one time-dependent covariate is enough to degrade the
            performances of the benchmarks, despite the fact that they have the
            exact same parametric form as the true model.
            
            Finally we find comparable performance among the ad-hoc boosted
            estimator and our proposed one, although a slight edge goes to the latter especially in the more difficult simulations with larger a. The results here demonstrate that there is a place in the survival
            boosting literature for fully nonparametric methods like this one that
            can flexibly handle time-dependent covariates.
            
            
            
            Discussion
            Our estimator can also potentially be used to evaluate the
            goodness-of-fit of simpler parametric hazard models. Since our
            approach is likelihood-based, future work might examine whether model
            selection frameworks like those in <cit.> can be extended to
            cover likelihood functionals. For this, <cit.> provides
            some guidance for determining the effective degrees of freedom for the
            boosting estimator. The ideas in <cit.> may also be germane.
            
            The implementation presented in Section <ref> is one of many possible ways to implement our estimator. We defer the design of a more refined implementation to future research, along with open-source code.
            10pt Acknowledgements. The review team provided many insightful comments that significantly improved our paper. We are grateful to Brian Clarke, Jack Hall, Sahand Negahban, and Hongyu Zhao for helpful discussions. Special thanks to Trevor Hastie for early formative discussions. The dataset used in Section <ref> was kindly provided by Dr. Kito Lord.
            
            
            
            40pt
             APPENDIX: PROOFS
            
             §.§ Proof of Proposition <ref>
            
            	Writing
            	
                R(F)=( ∫_0^1Y(t)· e^F(t,X(t)) dt-Δ F(T,X(T)) ),
            
            	we can apply (<ref>) to establish the first part of the
            	integral in (<ref>) when F∈∪{logλ}.  To
            	complete the representation, it suffices to show that the point
            	process
            	
                M(B)=Δ· I[{T,X(T)}∈ B]
            
            	has mean ∫_Bλ dμ, and then apply Campbell's formula. To
            	this end, write N(t)=I(T≤ t) and consider the filtration σ{X(s),Y(s),N(s):s≤ t}. Then N(t) has the Doob-Meyer form
            	dN(t)=λ(t,X(t))Y(t)dt+dM(t) where M(t) is a
            	martingale. Hence
            	{ M(B) }   =   ( ∫_0^1I[{t,X(t)}∈ B] dN(t))
               =   ( ∫_0^1 Y(t)· I[{t,X(t)}∈ B]·λ(t,X(t)) dt )
                  + (∫_0^1I[{t,X(t)}∈ B] dM(t))
               =   ∫_Bλ dμ
            	+(∫_0^1I[{t,X(t)}∈ B] dM(t)),
            	
            	where the last equality follows from (<ref>). Since
            	I[{t,X(t)}∈ B] is predictable because X(t) is, the desired
            	result follows if the stochastic integral ∫_0^1I[{t,X(t)}∈ B]dM(t) is a martingale. By Section 2 of Aalen <cit.>, this is true if M(t) is square-integrable. In fact, M(t)=N(t)-∫_0^tλ(t,X(t))dt is bounded because λ(t,x) is bounded above by <ref>. This establishes
            	(<ref>).
            	
            	Now note that for a positive constant Λ the function e^y-Λ y
            	is bounded below by both -Λ y and Λ y+2Λ{1-log2Λ},
            	hence e^y-Λ y≥Λ|y|+2Λmin{0,1-log2Λ}.
            	Since Λmin{0,1-log2Λ} is non-increasing in Λ, <ref> implies that
            	
                e^F(t,x)-λ(t,x)F(t,x)    ≥min{ e^F(t,x)-Λ_LF(t,x),e^F(t,x)-Λ_UF(t,x)}
                   ≥Λ_L|F(t,x)|+2Λ_Umin{0,1-log(2Λ_U)}.
            
            	Integrating both sides and using the norm equivalence relation (<ref>)
            	shows that
            	
                R(F)    ≥Λ_LF_μ,1+2Λ_Umin{0,1-log(2Λ_U)}
                   ≥2Λ_L/F_∞+2Λ_Umin{0,1-log(2Λ_U)}
                   ≥2Λ_L/F_μ,2+2Λ_Umin{0,1-log(2Λ_U)}.
            
            	The lower bound (<ref>) then follows from the second
            	inequality. The last inequality shows that R(F) is coercive on
            	(,⟨·,·⟩ _μ). Moreover the
            	same argument used to derive (<ref>) shows that R(F) is
            	smooth and convex on (,⟨·,·⟩
            	_μ).  Therefore a unique minimizer F^* of R(F) exists in
            	(,⟨·,·⟩ _μ).  Since <ref>
            	implies there is a bijection between the equivalent classes of
            	(,⟨·,·⟩ _μ) and the functions
            	in , F^* is also the unique minimizer of R(F) in
            	. Finally, since e^F(t,x)-λ(t,x)F(t,x) is pointwise
            	bounded below by λ(t,x){1-logλ(t,x)},
            	R(F)≥∫(λ-λlogλ)dμ=R(logλ) for all
            	F∈.
            
            
             §.§ Proof of Lemma <ref>
            
            	By a pointwise-measurable argument (Example 2.3.4 of <cit.>) it
            	can be shown that all suprema quantities appearing below are
            	sufficiently well behaved, so outer integration is not
            	required. Define the Orlicz norm
            	X=inf{C>0:Φ(|X|/C)≤ 1} where
            	Φ(x)=e^x^2-1.  Suppose the following holds:
            	sup_F∈_^Ψ_n|{_n(F)-_n(0)}-{R(F)-R(0)}|    ≤   κ' J__/n^1/4, 
            sup_G∈_: G_∞≤ 1|G_,1-G_μ,1|    ≤   κ”J__/n^1/2,
            	where J__ is the complexity measure (<ref>),
            	and κ',κ” are universal constants. Then by Markov's
            	inequality, (<ref>) holds with probability
            	at least 1-2exp[-{η n^1/4/(κ'J__)}^2], and
            	sup_G∈_: G_∞≤ 1{G_μ,1-G_,1} <1/
            holds with probability at least
            1-2exp[-{n^1/2/(κ”J__)}^2].  Since
            >1 and η<1, (<ref>) and
            (<ref>) jointly hold with probability at least
            1-4exp[-{η n^1/4/(κ J__)}^2].  The
            lemma then follows if (<ref>)
            implies (<ref>). Indeed, for any non-zero F∈,
            its normalization G=F/F_∞ is in _ by
            construction (<ref>). Then (<ref>) implies that
            
                F_∞/F_,1=1/G_,1≤
            
            because
            
                1/>G_μ,1-G_,1≥2/-G_,1,
            
            where the last inequality follows from the definition of  (<ref>).
            
            Thus it remains to establish (<ref>) and (<ref>),
            which can be done by applying the symmetrization and maximal
            inequality results in Sections 2.2 and 2.3.2 of <cit.>. Write
            _n(F)=(1/n)∑_i=1^nl_i(F) where
            l_i(F)=∫_0^1Y_i(t)e^F(t,X_i(t))dt-Δ_iF(T_i,X_i(T_i))
            are independent copies of the loss
            
                l(F)=∫_0^1 Y(t)· e^F(t,X(t))dt-Δ· F(T,X(T)),
            
            which is a stochastic process indexed by F∈. As was shown in
            Proposition <ref>, {l(F)}=R(F). Let
            ζ_1,⋯,ζ_N be independent Rademacher random
            variables that are independent of Z={(X_i(·),Y_i(·),T_i)}_i=1^n. It
            follows from the symmetrization Lemma 2.3.6 of <cit.> for
            stochastic processes that the left hand side of (<ref>) is
            bounded by twice the Orlicz norm of
                  -30pt
            sup_F∈_^Ψ_n|1/n∑_i=1^nζ_i{l_i(F)-l_i(0)}|
            
                  ≤1/nsup_F∈_^Ψ_n|∑_i=1^nζ_i∫_0^1
            Y_i(t) { e^F(t,X_i(t))-1} dt| 
                     + 1/nsup_F∈_^Ψ_n|∑_i=1^nζ_iΔ_iF(T_i,X_i(T_i))|. 
            
            Now hold Z fixed so that only ζ_1,⋯,ζ_n are
            stochastic, in which case the sum in the second line
            of (<ref>) becomes a separable subgaussian
            process. Since the Orlicz norm of ∑_i=1^nζ_ia_i is
            bounded by (6∑_i=1^na_i^2)^1/2 for any constant
            a_i, we obtain the following the Lipschitz property for any
            F_1,F_2∈_^Ψ_n:
                  -30pt
            ∑_i=1^nζ_i∫_0^1Y_i(t){ e^F_1(t,X_i(t))-e^F_2(t,X_i(t))} dtζ|Z
                  ≤ 6∑_i=1^n[∫_0^1 Y_i(t){ e^F_1(t,X_i(t))-e^F_2(t,X_i(t))} dt]^2
                  ≤ 6e^2Ψ_n∑_i=1^n(∫_0^1 Y_i(t)· |F_1(t,X_i(t))-F_2(t,X_i(t))|dt)^2
                  ≤ 6e^2Ψ_n∑_i=1^n∫_0^1Y_i(t){F_1(t,X_i(t))-F_2(t,X_i(t))}^2dt
                  =    6ne^2Ψ_nF_1-F_2_,2^2,
            
            where the second inequality follows from |e^x-e^y|≤ e^max(x,y)|x-y|
            and the last from the Cauchy-Schwarz inequality. Putting the Lipschitz
            constant (6n)^1/2e^Ψ_n obtained above into Theorem 2.2.4
            of <cit.> yields the following maximal inequality: There is
            a universal constant κ' such that
                  -30pt
            sup_F∈_^Ψ_n|∑_i=1^nζ_i∫_0^1Y_i(t)
            	{ e^F(t,X_i(t))-1} dt| ζ|Z
                  ≤κ' n^1/2e^Ψ_n∫_0^Ψ_n{log𝒩(u,_^Ψ_n,)} ^1/2du
                  ≤κ' n^1/2e^Ψ_nΨ_nJ__,
            
            where the last line follows from (<ref>). Likewise the
            conditional Orlicz norm for the supremum of
            |∑_i=1^nζ_iΔ_iF(T_i,X_i(T_i))|
            is bounded by κ' J__n^1/2Ψ_n.  Since neither
            bounds depend on Z, plugging back into (<ref>)
            establishes (<ref>):
                  -30pt
            sup_F∈_^Ψ_n|{_n(F)-_n(0)}-{R(F)-R(0)}| 
                  ≤ 2κ' J__Ψ_ne^Ψ_n/n^1/2{1+e^-Ψ_n}
                  ≤4κ' J__/n^1/4,
            
            where Ψ_ne^Ψ_n=n^1/4 by (<ref>).  On noting that
            
                G_,1 =
                1/n∑_i=1^n∫_0^1Y_i(t)|G(t,X_i(t))|dt,10pt
                G_μ,1={∫_0^1Y(t)|G(t,X(t))| dt},
            (<ref>) can be established using the same approach.
            
            
             §.§ Proof of Lemma <ref>
            
            	For m<, applying (<ref>) to
            	_n(_m+1)=_n(_m-ν_n/m+1g__m^)
            	yields
            	
                _n(_m+1)
                			   =
                			_n(_m)-ν_n/m+1⟨
                			g__m,g__m^⟩_
                      +ν_n^2/2(m+1)^2∫(g__m^)^2exp{_m-ρ (_m+1 - _m)   } d
                   <_n(_m)-ν_n/m+1g__m_,2+ν_n^2e^Ψ_n/2(m+1)^2,
            
            	where the bound for the second term is due to (<ref>)
            	and the bound for the integral follows from
            	∫(g__m^)^2d=1
            	(Definition <ref> of an -gradient) and
            	_m_∞, _m+1_∞<Ψ_n for
            	m< (lines 5-6 of Algorithm <ref>). Hence for
            	m≤, (<ref>) implies that
            	
                _n(_m)<_n(0)
                	+∑_m=0^∞ν_n^2e^Ψ_n/2(m+1)^2<_n(0)+1≤2
            
            	because ν_n^2 e^Ψ_n < 1 under (<ref>). Since
            	max_m≤m̂_m_∞<Ψ_n, and using our
            	assumption sup_F∈_^Ψ_n|_n(F)-R(F)|<1 in
            	the statement of the lemma, we have
            	
                R(_m)
                	≤_n(_m)+|_n(_m)-R(_m)|
                	<3.
            
            	Clearly the minimizer F^* also satisfies R(F^*)≤ R(0)<3.
            	Thus coercivity (<ref>) implies that
            	
                _m_∞,F^*_∞<α_β_Λ,
            
            	so the gap γ̂ defined in (<ref>) is bounded
            	as claimed.
            	
            	It remains to establish (<ref>), for
            	which we need only consider the case
            	_n(_m̂)-_n(F^*)>0.  The termination
            	criterion g__m=0 in Algorithm <ref> is never
            	triggered under this scenario, because by Proposition <ref> this
            	would imply that _ minimizes _n(F)
            	over the span of {_nj(t,x)}_j, which also contains F^* (Remark <ref>). Thus either
            	=∞, or the termination criterion _-ν_n/+1ĝ__^_∞≥Ψ_n in line 5 of Algorithm <ref> is met. In the latter case
            	
                Ψ_n   ≤‖_-ν_n/+1g__^‖_∞≤‖_-ν_n/+1g__^‖ _,2
                   ≤α_(∑_m=0^-1ν_n/m+1+1)
            
            	where the inequalities follow from (<ref>) and from
            	g__m^_,2=1. Since the sum is diverging, the inequality also holds for  sufficiently large (e.g. =∞).
            	
            	Given that F^* lies in the span
            	of {_nj(t,x)}_j, the Taylor expansion (<ref>) is valid for _n(F^*).
            	Since the remainder term in the expansion is non-negative, we have
            	_n(F^*)
            	   =   _n(_m + F^* - _m)
            
            	
            	
            	   ≥   _n(_m)  + ⟨
            	g__m,F^*-_m⟩_.
            	
            	Furthermore for m ≤m̂,
            	⟨ g__m,_m-F^*⟩_   ≤   _m-F^*_,2·g__m_,2
               ≤   _m-F^*_∞·g__m_,2
               ≤   γ̂·g__m_,2.
            	
            	Putting both into (<ref>) gives
            	_n(_m+1)
            	   <   _n(_m)+ν_n/γ̂(m+1)⟨ g__m,F^*-_m⟩_
            	+ν_n^2e^Ψ_n/2(m+1)^2
               ≤   _n(_m)+ν_n/γ̂(m+1){_n(F^*)-_n(_m)}+ν_n^2e^Ψ_n/2(m+1)^2.
            	
            	Subtracting _n(F^*) from both sides above and denoting
            	δ_m=_n(_m)-_n(F^*),
            	we obtain
            	
                δ_m+1
                	<(1-ν_n/γ̂(m+1))δ_m
                	+ν_n^2e^Ψ_n/2(m+1)^2.
            
            	Since the term inside the first parenthesis is between 0 and 1, solving
            	the recurrence yields
            	δ_   <   δ_0∏_m=0^-1(1-ν_n/γ̂(m+1))
            	+ν_n^2e^Ψ_n∑_m=0^∞1/2(m+1)^2
               ≤   max{0,δ_0}exp(-/γ̂∑_m=0^-1ν_n/m+1)+ν_n^2e^Ψ_n
               ≤   
            	emax{0,δ_0}exp(-/γ̂Ψ_n)+ν_n^2e^Ψ_n,
            	
            	where in the second inequality we used the fact that 0≤1+y≤ e^y
            	for |y|<1, and the last line follows from (<ref>).
            	
            	The Lambert function (<ref>) in Ψ_n=W(n^1/4) is
            	asymptotically log y-loglog y, and in fact by Theorem 2.1
            	of <cit.>, W(y)≥log y-loglog y for y≥ e. Since
            	by assumption n≥55>e^4, the above becomes
            	
                δ_
                	<emax{0,δ_0}(log  n/4n^1/4)^/(γ̂)+ν_n^2e^Ψ_n.
            
            	The last step is to control δ_0, which is bounded
            	by 1-_n(F^*) because _n(_0)=_n(0)≤1.
            	Then under the hypothesis |_n(F^*)-R(F^*)|<1,
            	we have
            	
                δ_0≤ 1-R(F^*)+1<2-R(F^*).
            
            	Since (<ref>)
            	implies R(F^*)≥2Λ_Umin{0,1-log(2Λ_U)},
            	
                δ_0
                	<2-R(F^*)≤2
                	+2Λ_Umax{0,log(2Λ_U)-1}<2β_Λ.
            
            
             §.§ Proof of Proposition <ref>
            
            	
            	Let δ=log n/(4n^1/4) which is less than one for n≥ 55 >
            	e^4.  Since ,γ̂≥1 it follows that
            	
                δ<(log n/4n^1/4)^/(γ̂).
            
            	Now define the following probability sets
            	
            	S_1   =   {sup_F∈_^Ψ_n|{_n(F)-_n(0)}-{R(F)-R(0)}|<δ/3}
            
            	S_2   =   {|_n(0)-R(0)|<δ/3}
            
            	S_3   =   {|_n(F^*)-R(F^*)|<δ/3}
            
            	S_4   =   {} ,
            	
            	and fix a sample realization from ∩_k=1^4S_k. Then the
            	conditions required in Lemma <ref> are satisfied with
            	sup_F∈_^Ψ_n|_n(F)-R(F)|<2δ/3, so
            	γ̂ (and hence _∞) is bounded and (<ref>) holds. Since Algorithm <ref> ensures that __∞<Ψ_n, we have _∈_^Ψ_n and therefore it also follows that |_n(_)-R(_)|<2δ/3. Combining (<ref>) and (<ref>) gives
            	‖_-F^*‖ _μ,2^2   ≤   2/min_t,x(λ^*∧)(2δ/3+δ/3
            	+ {_n(_)-_n(F^*)})
               <   2/min_t,x(λ^*∧)(δ+2eβ_Λ(log n/4n^1/4)^/(γ̂)
            	+1/16·log n/4n^1/4)
               <   13β_Λ/min_t,x(λ^*∧)(log n/4n^1/4)^/(γ̂),
            	
            	where the second inequality follows from (<ref>) and
            	ν_n^2 e^Ψ_n=log n/(64n^1/4), and the last from (<ref>). Now, using the inequality |e^x-e^y|≤max(e^x,e^y)|x-y| yields
            	
                ‖-ł^*‖ _μ,2^2
                	<13β_Łmax_t,x(ł^*∨)^2/min_t,x(ł^*∧)(log n/4n^1/4)^/(γ̂),
            
            	and the stated bound follows from F^*=logł since  is correctly specified (Proposition <ref>).
            	
            	
            	The next task is to lower bound (∩_k=1^4S_k).
            	It follows from Lemma <ref> that
            	
                (S_1∩ S_4)
                	≥
                	1-4exp{ -(log  n/12κ J__)^2} .
            
            	Bounds on (S_2) and (S_3) can be obtained using
            	Hoeffding's inequality. Note from (<ref>) that
            	_n(0)=∑_i=1^n∫_0^1 Y_i(t)dt/n and
            	_n(F^*)=∑_i=1^nl_i(F^*)/n for the loss
            	l(·) defined in (<ref>). Since 0≤∫_0^1
            	Y_i(t)dt≤1 and
            	-F^*_∞<l(F^*)≤e^F^*_∞+F^*_∞,
            	
                (S_2)
                	≥
                	1-2exp{ -2n^1/2(log n/12)^2},(S_3)
                	≥
                	1-2exp{ -2n^1/2(log n/36e^F^*_∞)^2} .
            
            	By increasing the value of κ and/or replacing J__
            	with max(1,J__) if necessary, we can combine the
            	inequalities to get a crude but compact bound:
            	
                {∩_k=1^4S_k}≥
                		1-8exp{ -(log  n/κ J__e^F^*_∞)^2}.
            
            	Finally, since F^*_∞ = logλ_∞ < max{|logΛ_L|,|logΛ_U|}, we can replace e^F^*_∞ in the probability bound above by Λ_L^-1∨Λ_U.
            
            
             §.§ Proof of Proposition <ref>
            
            	It follows from (<ref>) that ł^* is the orthogonal projection of
            	λ onto (,⟨·,·⟩ _μ).
            	Hence
            	‖-λ‖ _μ,2^2    =    ‖ e^F^*-λ‖ _μ,2^2
            	+‖ e^_-e^F^*‖ _μ,2^2
                =    min_F∈‖ e^F-λ‖ _μ,2^2
            	+‖ e^_-e^F^*‖ _μ,2^2
               ≤   min_F∈‖ e^F-λ‖ _μ,2^2
            	+max_t,x(λ^*∨)^2‖_-F^*‖ _μ,2^2,
            	
            	where the inequality follows from |e^x-e^y|≤max(e^x,e^y)|x-y|.
            	Bounding the last term in the same way as Proposition <ref>
            	completes the proof. To replace e^F^*_∞ in (<ref>) by Λ_L^-1∨Λ_U, it suffices to show that Λ_L ≤λ^*(t,x) ≤Λ_U. Since the value of λ^* over one of its piecewise constant regions B is ∫_Bλ dμ/μ(B), the desired bound follows from <ref>. We can also replace max_t,x(ł^*∨) and min_t,x(ł^*∧) with max_t,x(Λ_U∨) and min_t,x(Λ_L∧) respectively.
             REFERENCES
            abbrvnat

