Debiased Maximum-Likelihood Estimators of Hazard Ratios Under Kernel-Based Machine-Learning Adjustment

1. Introduction

The use of hazard ratios to measure the causal effect of treatment has recently come under debate. Although it has been standard practice in epidemiological studies to examine hazard ratios using the Cox proportional hazard model and its modified versions [1,2], several studies [3,4,5] have noted that hazard ratios are uninterpretable with regard to causation (see Martinussen (2022) [6] for a review). In particular, Martinussen et al. (2020) [5] provided concrete examples of data-generating processes in which the Cox model is correctly specified, but estimated hazard ratios are difficult to interpret. The main difficulty is that time courses of different study populations with different treatment effects are described by the same Cox model with the same set of parameter values. Although a few authors have rebutted the uninterpretability of hazard ratios in the Cox model [7,8,9], the issue concerning the unidentifiability with multiple, contradictory scenarios as posed by Martinussen et al. (2020) [5] remains unresolved.

To address this issue, researchers have sought alternative measures of causal treatment effect, such as differences between counterfactual survival functions or restricted mean survival time [10,11,12]. Methods for applying machine learning to estimate these measures have also been developed [13,14,15]. However, these measures are only applicable to simple, well-controlled settings, such as randomized clinical trials or observational studies of several covariate-adjusted groups with different baseline treatments, despite such studies making few assumptions about data-generating processes in other respects. Modern epidemiology increasingly requires methods for analyzing large quantities of observational data acquired in a more uncontrolled, dynamic manner, in which many covariates are measured in real-time. Examples of such data are electronic medical or health records and data acquired by electronic devices for health promotion [16]. The marginal structural Cox model potentially used for this purpose [17] suffers from the uninterpretability problem described above, and thus, development of an alternative method is required.

The present study proposes a strategy to use a hazard model based on time-dependent treatment variables and covariates, but with no unspecified baseline hazard. Previous studies [3,4,5] attributed the uninterpretability of hazard ratios to “selection,” the fact that less frail subjects (described by unobserved factors) are more likely to remain in the risk set at later stages of the study. Although these authors did not explicitly describe the role of the baseline hazard in the description of this selection process, it is the unspecified baseline hazard that allows the Cox model to be correctly specified for selection processes of many different patterns (see Remark 3-1 below). Thus, we abandon the baseline hazard, and instead exploit the descriptive power of machine learning to capture how the risk set changes over time. This approach can now be implemented due to recent advances in incorporating machine learning into rigorous statistical analysis with effect estimation [18,19]. Prior to this development, machine learning could be used for only outcome prediction in epidemiology, because estimation with machine learning was biased. We thus develop an algorithm that debiases maximum-likelihood (ML) estimators of hazard ratios in the model, using the framework of doubly robust, debiased machine learning (DML) based on Neyman orthogonality (see Ref.[20] for an introductory review). Notably, our algorithm applies to models with latent variables, hence enabling explicit modeling of the selection in the risk set caused by unobserved factors.

The article is organized as follows. In section 2-1, we introduce our problem setting for an observational study with uncontrolled, dynamic treatment and real-time measurement of covariates, together with an exponential, machine-learning-based hazard model used for analysis. In section 2.2, we show that the ML estimator of hazard ratios can be interpreted as a measure of causal treatment effect in this setting (Proposition 1), clarifying required (mostly testable) assumptions. In this argument, we show that the unindentifiability with multiple interpretations is alleviated (Remark 3-1). In section 2.3, we construct Neyman-(near)-orthogonal scores used for debiasing the ML estimators of hazard ratios in our model (Propositions 2 and 3). For these scores, additional convergence results required for DML are also provided in Appendix H. In section 2.4, we develop an algorithm for computing the debiased estimators of hazard ratios after model selection, suitably adapting the procedure proposed in Ref.[19] to our setting. Then, in section 3, we apply the developed method to two sets of clinically plausible simulation data. In the first simulation study described in section 3.1, we show that our estimators appropriately identify the true hazard ratios with minimal bias, under the effect of complicated nonlinear confounding among treatment, comorbidities and outcome. In the second simulation study described in section 3.2, we present a case in which an unobserved factor causes selection in the risk set, resulting in less frail subjects in the later phase of the study. Our debiased estimator based on a model with a latent variable appropriately identify the true hazard ratio with minimal bias. On the basis of these results, in section 4, we discuss multiple advantages of the proposed ML approach over the conventional Cox’s maximum-partial-likelihood (MPL) approach as well as its limitations. We list mathematical notations in Table A1.

2. Theories and Methods

2.1. Exponential Parametric Hazard Model Combined with Machine Learning

We consider observational studies in which the occurrence or absence of an event of interest in subjects randomly sampled from a large population and indexed by $i(\in \mathcal{I})$ is longitudinally observed over time (indexed by $t\in \mathcal{T}\subset \mathbf{R}$) during a noninformatively right-censored period, $0\le t\le C_i$. With a set of time-dependent covariates, collectively denoted by $X_{i,t}(={\left\{X_{ij,t}\right\}}_{j\in \mathcal{J}})\in \mathcal{X}\subset {\mathbf{R}}^d$$(d\in \mathbf{N})$, we strive to identify the effect of time-dependent treatment described by a collection of binary variables $A_{i,t}(={\left\{A_{ik,t}\right\}}_{k\in \mathcal{K}}\in \mathcal{A}\subset {\{0,1\}}^{\left|\mathcal{K}\right|}$). For clarity of presentation, we assume that at most a single treatment variable can take a value of unity, and the remaining variables must be zero. Consider an exponential hazard model whose conditional hazard is given by

$$\begin{array}{c}\hfill h\left(t\right|A_{i,t},X_{i,t})\stackrel{\mathrm{def}}{=}exp\left({\theta }^{\prime }{A}_{i,t}+f\left(X_{i,t}\right)\right),\end{array}$$

(1)

where $\theta ={\left\{{\theta }_k\right\}}_{k\in \mathcal{K}}\in \Theta \subset {\mathbf{R}}^{\left|\mathcal{K}\right|}$ is a set of parameters corresponding to the natural logarithm of the hazard ratio of the untreated samples to samples treated to different extents, and the function f describes the risk variation due to the given set of covariates, which we adjust using a machine learning model (i.e., a set of functions $\mathcal{M}$ in which f lies). In addition, hereinafter, ${(·)}^{\prime }$ denotes the transposition of vectors and matrices. Note that the above model lacks the baseline hazard as a function of t (the time elapsed after enrollment). The covariate $X_{i,t}$ can include temporal information, such as date and durations of suffered disorders, but not t itself.

Suppose that a subset of subjects (indexed by i) experiences an event of interest at time $T_i$ in the observation period, and the other subjects experience no such event. The full log-likelihood for this observation is then given by

$$\begin{array}{c}\hfill ln{Q}_h\left(T\right|A,X)=\sum _{i\in \mathcal{I}}\left\{I_{[0,C_i]}\left(T_i\right)({\theta }^{\prime }{A}_{i,T_i}+f\left(X_{i,T_i}\right))-{\int }_0^{T_i\wedge C_i}exp\left({\theta }^{\prime }{A}_{i,t}+f\left(X_{i,t}\right)\right)dt\right\},\end{array}$$

(2)

(see, e.g., Ref.[21]), where $I_{[0,C_i]}\left(T_i\right)$ is the indicator function that returns unity for $T_i\in [0,C_i]$ and zero otherwise. To see the meaning of Equation (2), we consider how the hazard model describes the occurrence of an event in a small interval $(t,t+\delta t]$. Assuming for the moment that (i) $X_{i,s}$ is continuous with respect to s and (ii) $A_{i,s}=const$ over $s\in [t,t+\delta t]$, the likelihoods of occurrence and non-occurrence of an event in this period conditioned on event-free survival at time t are $h\left(t\right|A_{i,t},X_{i,t})\delta t$ and $1-h\left(t\right|A_{i,t},X_{i,t})\delta t$, respectively, up to the first order in $\delta t$. Simple computations show that the log-likelihood for this binary observation is

$$\begin{array}{ccc}\hfill ln{Q}_h\left(I_{(t,t+\delta t]}\left(T_i\right)\right|T_i>t,A_{i,t},X_{i,t})& =& I_{(t,t+\delta t]}\left(T_i\right)lnh\left(t\right|A_{i,t},X_{i,t})\delta t\hfill \\ & & +(1-I_{(t,t+\delta t]}\left(T_i\right))ln(1-h\left(t\right|A_{i,t},X_{i,t})\delta t)+O\left(\delta t^2\right),\hfill \\ & & =I_{(t,t+\delta t]}\left(T_i\right)({\theta }^{\prime }{A}_{i,t}+f\left(X_{i,t}\right))-exp\left({\theta }^{\prime }{A}_{i,t}+f\left(X_{i,t}\right)\right)\delta t\hfill \\ & & +O\left(\delta t^2\right)+const.\hfill \end{array}$$

(3)

The time integration of Equation (3) results in the summand of Equation (2) up to a constant independent of $\theta$ and f.

We assume that the two assumptions described above hold almost surely in the $\delta t\to +0$ limit, because ${\left\{A_{it}\right\}}_t$ and ${\left\{X_{it}\right\}}_t$ are continuous with respect to time except for finite numbers of time points in real applications. Here, however, we make it clear in advance that we do not go into the details about the necessary conditions on the data-generating stochastic process, ${\{A_{it},X_{it}\}}_t$. We do not specify the filtered probability space for the process, but rather assume the existence of the marginal (conditional) probability measures on the measurable spaces $(\mathcal{A}\times \mathcal{X}\times \mathcal{T},{\mathfrak{M}}_{\mathcal{A}}\times {\mathfrak{M}}_{\mathcal{X}}\times {\mathfrak{M}}_{\mathcal{T}})$, and the integrability of quantities of interest with respect to these measures. For the model $\mathcal{M}$ on the base space $(\mathcal{X},{\mathfrak{M}}_{\mathcal{X}})$ for covariates, we consider only reproducing-kernel Hilbert spaces (RKHS) ${\mathcal{H}}_k$ associated with a positive-semidefinite bounded kernel $k(·,·)$ [22]. We refer readers to Ref.[23] which provides a basic tool for our theoretical analysis of kernel-based machine learning.

2.2. Causal Inference and ML Estimation

This section clarifies the causal interpretation of the estimation results in the model described in the previous section and discusses the necessary assumptions that support this interpretation. First, we make the conventional assumptions of consistency, no unmeasured confounder (ignorability), and positivity.

Assumption 1 (consistency). 

For any counterfactual treatment schedule a, we have, $X_{(-\infty ,t]}^a=X_{(-\infty ,t]}|A_{(-\infty ,t]}=a_{(-\infty ,t]}$ and $N_{(-\infty ,t]}^a=N_{(-\infty ,t]}|A_{(-\infty ,t]}=a_{(-\infty ,t]}$, for counting processes $N_t=I_{(-\infty ,t]}\left(T\right)$ and $N_t^a=I_{(-\infty ,t]}\left(T^a\right)$.

Assumption 2 (no unmeasured confounder). 

For some $ϵ>0$, we have, for any t and any counterfactual treatment schedule a,

$$\begin{array}{ccc}\hfill T^a\perp \perp A_{[t,t+ϵ]}& |& T^a>t,X_{(-\infty ,t+ϵ]},A_{(-\infty ,t)}.\hfill \end{array}$$

(4)

Assumption 3 (positivity) [optional]. 

For any $t>0$ and any $a\in \mathcal{A}$, we have $P(A_t=a|T>t,X_t)>0$.

Here, the set subscript, such as $X_{(-\infty ,t]}$, denotes the collection of all variables with subscripts included in the set. Throughout the paper, P denotes the marginal (conditional) probability measures of the argument variables derived from the data-generating process. We have used variables for $-\infty <t<0$ for a technical reason that $A_0$ and $X_0$ must be determined on the basis of the past information. This may be replaced by some finite interval before the enrollment at $t=0$ for which the information about treatment and covariate is available.

Next, in addition to the above, we make the following assumption about the regularity of hazard.

Assumption 4 (regularity of hazard).  

$$\begin{array}{ccc}& & \underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_{(-\infty ,t+\delta t]},X_{(-\infty ,t+\delta t]})\hfill \\ & =& \underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_{(-\infty ,t]},X_{(-\infty ,t]})\hfill \\ & <& \infty .\hfill \end{array}$$

(5)

holds almost surely.

Remark 1-1. 

Assumption 1 is a natural assumption without which we have two different processes for the same treatment schedule.

Remark 1-2. 

Assumption 2 is stating that the assignment of treatment at each time point t can be considered as randomized one conditioned on the value of covariate. In the analysis of electronic medical records, for example, if the medical practitioners decide their treatment on the basis of evidence recorded as covariates and those covariates are measured in real time, this assumption is valid. It should be noted that if the decision is based on past records, a past variable should be combined with the present variable to form a vectorized covariate. Causal inference with the same assumption with discretized timesteps has been performed with marginal structural Cox models (see Ref.[17,24] for examples).

Remark 1-3. 

Assumption 3 is a technical assumption required for debiasing with inverse-treatment-probability weights (see Proposition 3 and Remark 5-1 for our model and Ref.[25] for marginal structural Cox models). If medical practioners never choose a particular treatment option for a value of covariate, this assumption is violated. Thus, if one analyzes the effect of a treatment considered as a taboo for comorbidities included in the set of covariates (for example, anticholinergics for angle-closure glaucoma), they need to take an extra care of this assumption. This assumption can be omitted if one uses Proposition 2 for estimation in our study.

Remark 1-4. 

Assumption 2 essentially implies that the treatment at time t does not affect the time evolution of covariate X for a subsequent period $(t,t+ϵ]$ as seen through the following argument. First, we have

$$\begin{array}{ccc}& & P\left(I_{(t,t+ϵ]}\left(T^a\right)\right|T^a>t,X_{(-\infty ,t+ϵ]},A_{(-\infty ,t)}=a_{(-\infty ,t)})\hfill \\ & =& P\left(I_{(t,t+ϵ]}\left(T^a\right)\right|T^a>t,X_{(-\infty ,t+ϵ]},A_{(-\infty ,t+ϵ]}=a_{(-\infty ,t+ϵ]})\hfill \\ & =& P\left(I_{(t,t+ϵ]}\left(T\right)\right|T>t,X_{(-\infty ,t+ϵ]},A_{(-\infty ,t+ϵ]}=a_{(-\infty ,t+ϵ]}),\hfill \end{array}$$

(6)

where the first and second equalities are due to no unmeasured confounder and consistency, respectively. If the time evolution of covariates is immediately affected by the treatment schedule, and if outcome is immediately affected by the covariates, then under the conditioning on $X_{[t,t+ϵ]}$, $X_{[t,t+ϵ]}^a$ depends on the value of A, and hence, $T^a$ depends on A, which contradicts the assumption. The removal of such instantaneous interactions among treatment, covariates and outcome also makes Assumption 4 reasonable.

In addition to Assumptions 1–4, we assume that the model in Equation (1) is correctly specified. More precisely, we make the following four assumptions:

Assumption 5 (hazard conditionally independent of past treatment and covariate). 

For any $t>0$, $a\in \mathcal{A}$, $x\in \mathcal{X}$ and $E_A,E_X\subset (-\infty ,t)$, we almost surely have,

$$\begin{array}{c}\hfill \underset{\delta t\to +0}{lim}{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a,X_t=x)}{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a,X_t=x,A_{E_A},X_{E_X})}}=1.\end{array}$$

(7)

Assumption 6 (homogeneous treatment effect). 

For any $a,a^{\prime }\in \mathcal{A}$ and any $x\in \mathcal{X}$, the hazard contrast,

$$\begin{array}{c}\hfill \underset{\delta t\to +0}{lim}{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a,X_t=x)}{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a^{\prime },X_t=x)}},\end{array}$$

(8)

takes a constant value regardless of the value of x.

Assumption 7 (time-homogeneous hazard): 

The hazard is independent of its timing during the observation period. Mathematically, for any $t,t^{\prime }\ge 0$, and any fixed $a\in \mathcal{A}$ and $x\in \mathcal{X}$,

$$\begin{array}{c}\hfill \underset{\delta t\to +0}{lim}{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a,X_t=x)}{P(I_{(t^{\prime },t^{\prime }+\delta t]}\left(T\right)=1|T>t^{\prime },A_{t^{\prime }}=a,X_{t^{\prime }}=x)}}=1.\end{array}$$

(9)

Assumption 8 (correctly specified machine-learning model): 

There exist unique ${\theta }^{*}\in \mathbf{R}$ and $f^{*}\in \mathcal{M}$ that satisfy

$$\begin{array}{c}\hfill \mathrm{E}\left[ln{Q}_{h_{{\theta }^{*},f^{*}}}\left(T\right|A,X)\right]=\underset{\theta \in \mathbf{R},f\in {\mathfrak{B}}_{\mathcal{X}}}{sup}\mathrm{E}\left[ln{Q}_{h_{\theta ,f}}\left(T\right|A,X)\right],\end{array}$$

(10)

where ${\mathfrak{B}}_{\mathcal{X}}$ is the set of Borel-measurable functions.

Remark 2. 

Assumptions 5–7 can be validated by extending the model to incorporate the effect of past treatment and covariates, inhomogeneous treatment effect, and time inhomogeneity, respectively. The original and extended models can then be compared in terms of, for example, Bayesian model evidence (BME) (see, e.g., Chapter 3 of Ref.[26]). If Assumption 5 is violated, the definition of the treatment variables and covariates may be extended so that the current value retains past information. This is the same strategy as employed for designing variables in a marginal structural Cox model [17]. If Assumption 6 is violated, the model may be extended by incorporating covariate-dependent heterogeneous treatment effect. Violation of Assumption 7 indicates that variations in covariate values do not account for the temporal change in the risk set. In this case, additional covariates or extending the model with latent variables (as described below) may be considered. Finally, Assumption 8 is an a priori assumption that stems from the fact in statistical learning that inference in the function space requires a regularity assumption.

Proposition 1. 

Let the maximizer of the expected likelihood be

$$\begin{array}{c}\hfill {\theta }^{*},f^{*}=\arg \underset{\theta \in \mathbf{R},f\in \mathcal{M}}{max}\mathrm{E}[ln{Q}_h\left(T\right|A,X)].\end{array}$$

(11)

Under Assumptions 1–8, suppose two possible treatment schedules $a,a^{\prime }$ such that $a_{(-\infty ,t)}=a_{(-\infty ,t)}^{\prime }$ holds, and that $a_{k,s}=1$ and $a_s^{\prime }=0$ (or $a_{k^{\prime },s}^{\prime }=1$) holds for $s\in [t,t+ϵ]$ for some $ϵ>0$. Then, we have, for any sample path $x_{(-\infty ,t+\delta t]}$ for covariate,

$$\begin{array}{ccc}& & {\theta }_k^{*}(-{\theta }_{k^{\prime }}^{*})=\hfill \\ & & \underset{\delta t\to +0}{lim}ln{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T^a\right)=1|T^a,T^{a^{\prime }}>t,A_{(-\infty ,t)}=a_{(-\infty ,t)},X_{(-\infty ,t+\delta t]}=x_{(-\infty ,t+\delta t]})}{P(I_{(t,t+\delta t]}\left(T^{a^{\prime }}\right)=1|T^a,T^{a^{\prime }}>t,A_{(-\infty ,t)}=a_{(-\infty ,t)}^{\prime },X_{(-\infty ,t+\delta t]}=x_{(-\infty ,t+\delta t]})}}.\hfill \end{array}$$

(12)

Proof. 

see Appendix C.    □

Remark 3-1. 

Suppose that the entire covariate $X_t$ is unobserved. In this case, the modified version of the above proof indicates that a hazard model expanded with respect to time:

$$\begin{array}{c}\hfill h\left(t\right|A_t)\stackrel{\mathrm{def}}{=}exp(\theta {\left(t\right)}^{\prime }{A}_t+f\left(t\right))\end{array}$$

(13)

admits the optimal solution

$$\begin{array}{ccc}\hfill {\theta }^{*}\left(t\right)& =& ln\underset{\delta t\to +0}{lim}{\displaystyle \frac{{\int }_{\mathcal{X}}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_{k,t}=1,X_t)dP\left(X_t\right|A_{k,t}=1,T>t)}{{\int }_{\mathcal{X}}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=0,X_t)dP\left(X_t\right|A_t=0,T>t)}},\hfill \\ \hfill f^{*}\left(t\right)& =& ln\underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}{\int }_{\mathcal{X}}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=0,X_t)dP\left(X_t\right|A_t=0,T>t),\hfill \end{array}$$

where one can interpret $exp\left(f^{*}\left(t\right)\right)$ as a baseline hazard. The example of the uninterpretability of hazard ratios presented by Martinussen et al. [5] was for this setting. They assumed time-independent $A_t$ and $X_t$. Their point was that multiple combinations of $dP\left(X_t\right|A_t,T>t)$ and ${lim}_{\delta t\to +0}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t,X_t)/\delta t$ yields the same solution for $\theta \left(t\right)$.

This ambiguity has been removed in our setting. We assume that ${lim}_{\delta t\to +0}P(I_{(t,t+\delta t]}$$\left(T\right)=1|T>t,A_t,X_t)/\delta t$ is independent of specific time points of the observation period, which is biologically natural. Then, time-dependence of solutions for $\theta$ and f indicates lack of observation of relevant factors and encourages us to seek for a better model. This proof mechanism does not work in the Cox model, because the Cox model does not estimate the baseline hazard. Although our framework cannot exclude a scenario in which temporal changes in the distribution of unobserved factors miraculously keep ${\theta }^{*}\left(t\right)$ and $f^{*}\left(t\right)$ constant over time, this is not a generic case.

Remark 3-2. 

The quantity in Equation (12) is a covariate-adjusted hazard contrast for two counterfactual treatments that branch at time t, so it can be interpreted as a measure of causal effect. In the setting discussed in Ref.[17], this corresponds to the treatment effect measured in the next month (for which covariates at each month can be regarded as baseline covariates [27]).

Remark 3-3. 

See Remarks 8-1 and 8-2 in Appendix I for more consideration on causal interpretation of hazard ratios.

2.3. Debiasing ML Estimators

Next, we assume that we have consistent estimators for $\theta$ and f that are the maximizers of the empirical log-likelihood in Equation (2) in the presence of a suitable regularizer. Although a large body of the machine-learning literature presents such consistent estimators, they are biased in the sense that we cannot expect $\sqrt{n}\parallel \widehat{\theta }-{\theta }^{*}{\parallel }_2,\sqrt{n}{\parallel \widehat{f}-f^{*}\parallel }_{\mathcal{M}}\stackrel{p}{\to }0$ for the number of subjects $n(\to \infty )$. This situation prohibits making a decision on the significance of estimated results. However, we can apply to this problem debiased machine learning based on Neyman orthogonality and developed by Chernozhukov et al. (2018) [19]. Chernozhukov et al. (2018) describe how to systematically debias ML estimators and other M-estimators. Applying their idea to our problem setting yields debiased estimators of ${\theta }^{*}$.

Definition 1 (Neyman near-orthogonal scores). 

Suppose that an estimator $\widehat{\theta }(\in \Theta \subset {\mathbf{R}}^d)$ of quantities of interest ${\theta }^{*}$ is given as a zero of the empirical average of score functions $\varphi :\mathcal{D}\times \Theta \times \mathcal{N}\to {\mathbf{R}}^d$ for i.i.d. data $D_i\in \mathcal{D}$:

$$\begin{array}{c}\hfill \widehat{E}\left[\varphi (D;\widehat{\theta },\widehat{\eta })\right]\stackrel{\mathrm{def}}{=}{\displaystyle \frac{1}{n}}\sum _i\varphi (D_i;\widehat{\theta },\widehat{\eta })=0,\end{array}$$

(14)

where $\widehat{E}$ denotes empirical average and the third argument of ϕ is a (possibly infinite-dimensional) nuisance parameter and an estimator $\widehat{\eta }$ of an optimal value ${\eta }^{*}\in \mathcal{N}$ is used. Also suppose that there exist convex neighbourhoods of ${\eta }^{*}$, ${\mathcal{N}}_n\subset \mathcal{N}$ to which $\widehat{\eta }$ belongs with high probabilities converging to one. Then, ϕ is said to be “Neyman near-orthogonal,” if it satisfies

$$\begin{array}{ccc}\hfill \mathrm{E}\left[\varphi (D;{\theta }^{*},{\eta }^{*})\right]& =& 0,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill {\partial }_{\eta }\mathrm{E}\left[\varphi (D;{\theta }^{*},{\eta }^{*})\right][\eta -{\eta }^{*}]& \le & {ϵ}_n\sim o\left(n^{-1/2}\right){(}^{\forall }\eta \in {\mathcal{N}}_n),\hfill \end{array}$$

(16)

where, in the second line, ${\partial }_{\eta }(·)[\eta -{\eta }^{*}]$ denotes the Gateaux derivative operator (see Table A1 and Ref.[28]. for definition). If the above condition holds for ${ϵ}_n=0$, ϕ is said to be “Neyman orthogonal”.

Given a few conditions on the regularity of score functions and the convergence of nuisance parameters, Chernozhukov et al. [19] proved $\sqrt{n}\parallel \widehat{\theta }-{\theta }^{*}\parallel \stackrel{p}{\to }0$. In the following, we construct such orthogonal scores from the log-likelihood in Equation (2). For this purpose, we modify the following lemma proven by Chernozhukov et al. [19].

Lemma 1 (Lemma 2.5 of Chernozhukov et al. 2018: Neyman orthogonal scores derived from log-likelihood). 

Consider a ML problem,

$$\begin{array}{c}\hfill {\theta }^{*},f^{*}=arg\underset{\theta \in \Theta ,f\in \mathcal{M}}{min}\mathrm{E}[\ell (D;\theta ,f)],\end{array}$$

(17)

where ℓ is a sample-wise negative log-likelihood function and $\mathcal{M}$ is a convex set of high-dimensional vectors. Define $f_{\theta }$ with

$$\begin{array}{c}\hfill f_{\theta }=arg\underset{f\in \mathcal{M}}{min}\mathrm{E}[\ell (D;\theta ,f)],\end{array}$$

(18)

and a convex set $\mathcal{N}$ of mappings of Θ into $\mathcal{M}$ with an optimal solution in this set ${\eta }^{*}\left(\theta \right)=f_{\theta }$. Further suppose that for each $\eta \in \mathcal{N}$, the function $\theta \mapsto \ell (D;\theta ,\eta (\theta \left)\right)$ is continuously differentiable almost surely. Then, under mild regularity conditions, the score ϕ:

$$\begin{array}{c}\hfill \varphi (D;\theta ,\eta )={\displaystyle \frac{d\ell (D;\theta ,\eta (\theta \left)\right)}{d\theta }}\end{array}$$

(19)

is Neyman orthogonal at $({\theta }^{*},{\eta }^{*})$. Here, the differentiation on the right-hand side of the above is the full derivative.

For a finite-dimensional f, the above lemma with the application of the implicit function theorem to ${\partial }_{\theta }E[\ell (D;\theta ,f)]=0$ (or its empirical version) yields the following orthogonal score (section 2.2.1 of Ref.[19]):

$$\begin{array}{c}\hfill \varphi (D;\theta ,(f,\mu ))={\partial }_{\theta }\ell (D;\theta ,f)-\mu {\partial }_f\ell (D;\theta ,f),\end{array}$$

(20)

with ${\mu }^{*}=E[{\partial }_{\theta }{\partial }_f\ell {|}_{{\theta }^{*},f^{*}}]E{[{\partial }_f{\partial }_f\ell {|}_{{\theta }^{*},f^{*}}]}^{-1}$ and its estimator $\widehat{\mu }=\widehat{E}[{\partial }_{\theta }{\partial }_f\ell {|}_{\widehat{\theta },\widehat{f}}]\widehat{E}{[{\partial }_f{\partial }_f\ell {|}_{\widehat{\theta },\widehat{f}}]}^{-1}$ as long as the Hessians with respect to f are invertible.

In the following, for our high-dimensional machine-learning applications, we develop an analogous approach with an extra care.

Proposition 2 (Neyman near-orthogonal scores for ML in the present model). 

Suppose that the model $\mathcal{M}$ is given as a reproducing-kernel Hilbert space ${\mathcal{H}}_k$ associated with a bounded positive-semidefinite kernel k on $\mathcal{X}\times \mathcal{X}$. Then, for the sample-wise negative log-likelihood ℓ of the present study, the gradients ${\partial }_f\ell$, ${\partial }_f{\partial }_{{\theta }_k}\ell$ and the Hessian ${\partial }_f{\partial }_f\ell$ can be identified with elements of ${\mathcal{H}}_k$ and a linear operator from ${\mathcal{H}}_k$ into ${\mathcal{H}}_k$ in the following sense:

$$\begin{array}{ccc}\hfill \ell (D;\theta ,f+h)& =& \ell (D;\theta ,f)+{({\partial }_f\ell {|}_{\theta ,f},h)}_{\mathcal{H}}+{O(\parallel h\parallel }_{{\mathcal{H}}_k}^2),\hfill \\ \hfill {\partial }_{{\theta }_k}\ell (D;\theta ,f+h)& =& {\partial }_{{\theta }_k}\ell (D;\theta ,f)+{({\partial }_f{\partial }_{{\theta }_k}\ell {|}_{\theta ,f},h)}_{\mathcal{H}}+{O(\parallel h\parallel }_{{\mathcal{H}}_k}^2),\hfill \\ \hfill {\partial }_f{\ell |}_{\theta ,f+h}& =& {\partial }_f{\ell |}_{\theta ,f}+{\partial }_f{\partial }_f{\ell |}_{\theta ,f}h+{O(\parallel h\parallel }_{{\mathcal{H}}_k}^2).\hfill \end{array}$$

(21)

If Assumption 9 (stated below) holds, for ${\zeta }_n=c{n}^{-\alpha }$ with constants $c>0$ and $\alpha >{\displaystyle \frac{1}{2}}-\beta$, the score function,

$$\begin{array}{c}\hfill \varphi (D;\theta ,(f,H))={\partial }_{\theta }\ell (D;\theta ,f)-H_{\theta f}{(H_{ff}+{\zeta }_n)}^{-1}{\partial }_f\ell (D;\theta ,f),\end{array}$$

(22)

is Neyman near-orthogonal for the shrinking set of nuisance parameters ${\mathcal{N}}_n={\mathcal{B}}_n\times {\mathcal{S}}_n$ (defined below) with the true value $H^{*}=\mathrm{E}[{\partial }_{(\theta ,f)}{\partial }_{(\theta ,f)}\ell {|}_{{\theta }^{*},f^{*}}]$ and empirical estimators $\widehat{H}=\widehat{\mathrm{E}}[{\partial }_{(\theta ,f)}{\partial }_{(\theta ,f)}\ell {|}_{\widehat{\theta },\widehat{f}}]$. Here, $H_{\theta f}$ and $H_{ff}$ are blocks of the positive-semidefinite compact linear operator $H:{\mathbf{R}}^{\left|\mathcal{K}\right|}\times {\mathcal{H}}_k\to {\mathbf{R}}^{\left|\mathcal{K}\right|}\times {\mathcal{H}}_k$ with

$$\begin{array}{c}\hfill \left(\begin{array}{c}{H}_{\theta f}\\ H_{ff}\end{array}\right):\psi \in {\mathcal{H}}_k\mapsto H\left(\begin{array}{c}0\\ \psi \end{array}\right).\end{array}$$

(23)

The other blocks are similarly defined.

Proof. 

see Appendix D.    □

Assumption 9 

Estimator $(\widehat{\theta },\widehat{f})$ falls into ${\mathcal{B}}_n\stackrel{\mathrm{def}}{=}\left\{(\theta ,f)\right|\parallel \theta -{\theta }^{*}{\parallel }_2,\parallel f-f^{*}{\parallel }_{{\mathcal{H}}_k}\le c_1{n}^{-\beta }\}$ for some positive constants $c_1$ and $\beta (<1/2)$, with probabilities converging to one for $n\to \infty$. Estimator $\widehat{H}=\widehat{E}[{\partial }_{(\theta ,f)}{\partial }_{(\theta ,f)}\ell {|}_{\widehat{\theta },\widehat{f}}]$ falls in ${\mathcal{S}}_n\stackrel{\mathrm{def}}{=}\left\{H\right|$ positive semidefinite compact op. from ${\mathbf{R}}^{\left|\mathcal{K}\right|}\times {\mathcal{H}}_k$ into ${\mathbf{R}}^{\left|\mathcal{K}\right|}\times {\mathcal{H}}_k$s.t. $\parallel H-H^{*}\parallel \le c_2{n}^{-\beta }\}$ with probabilities converging to one for $n\to \infty$. It is further assumed that, for each $k\in \mathcal{K}$, $H_{f{\theta }_k}^{*}=H_{ff}^{*}{\rho }_k$ holds for some element ${\rho }_k$ of ${\mathcal{H}}_k$.

Remark 4-1. 

It can be argued that Assumption 9 is likely to hold (see Appendix F). However, we leave this as a conjecture, because its proof needs an explicit argument about the details of the data-generating stochastic process.

Remark 4-2. 

The assumption $H_{f{\theta }_k}^{*}=H_{ff}^{*}{\rho }_k$ can be relaxed to $H_{f{\theta }_k}^{*}=H_{ff}^{*\gamma }{\rho }_k$ ($0<\gamma <1$), which is an assumption sometimes made in the analysis of the kernel method [23]. In this case, the convergence rate should be suitably modified.

Remark 4-3. 

Apart from the Neyman near-orthogonality, Chernozhukov et al. [19] provided additional conditions on score regularity and quality of nuisance parameters (Assumptions 3.3 and 3.4 of their paper). Their proof of asymptotic unbiasedness of $\theta$ relies on these assumptions. In Appendix H, we prove that these conditions are satisfied for suitably chosen $\alpha$ and $\beta$.

The above score functions may not be intuitively understood. Inspecting the concrete representations of ${\partial }_{\theta }\ell$, ${\partial }_f\ell$, $H_{\theta f}$ and $H_{ff}$, one can derive the following more intuitively understandable formula as well.

Proposition 3 (Intuitive Neyman orthogonal scores for ML in our model). 

Assume that the positive Borel measures ${\mu }_a$ ($a\in \mathcal{A}$) defined as

$$\begin{array}{c}\hfill {\mu }_a\left(E\right)\stackrel{\mathrm{def}}{=}{\int }_0^{\infty }{\int }_0^C{\int }_{E}P(A_t=a|X_t,T>t)dP\left(X_t\right|T>t)P(T>t)dtdP\left(C\right)\end{array}$$

(24)

$(E\in \mathfrak{M}(\mathcal{X}):$ Borel measurable sets) are absolutely continuous to each other and have Radon-Nykodym derivatives

$$\begin{array}{c}\hfill {\displaystyle \frac{d{\mu }_k}{d{\mu }_0}}=e^{g_k^{*}},\end{array}$$

(25)

for which ${\mu }_k$ denotes ${\mu }_a$ such that $a_k=1$, and $g_k^{*}$ is a measurable function on $\mathcal{X}$. Then, the score functions for estimating θ,

$$\begin{array}{ccc}\hfill {\varphi }_k(D;\theta ,(f,g_k))& =& I_{[0,C]}\left(T\right)A_{k,T}{e}^{-{\theta }_k}(1+e^{-g_k\left(X_t\right)})-{\int }_0^{T\wedge C}{A}_{k,t}{e}^{f\left(X_t\right)}(1+e^{-g_k\left(X_t\right)})dt\hfill \\ & & +{\int }_0^{T\wedge C}(1-\sum _{\ell \in \mathcal{K}}{A}_{\ell ,t})e^{f\left(X_t\right)}(1+e^{g_k\left(X_t\right)})dt-I_{[0,C]}\left(T\right)(1-\sum _{\ell \in \mathcal{K}}{A}_{\ell ,T})(1+e^{g_k\left(X_T\right)}),\hfill \end{array}$$

(26)

$(k\in \mathcal{K})$, are Neyman orthogonal. Here, D denotes $(T,C,X,A)$.

Proof. 

See Appendix G.    □

Remark 5-1. 

As one might notice, the crux of the construction of the above orthogonal scores lies in balancing the second and third terms of Equation (26) by $g_k$. Roughly speaking, for each value of covariates, $g_k$ weights the integrand with the inverse-treatment-probability for $A_k=1$ and $A=0$. Here, note that this treatment probability is marginalized with respect to time.

Remark 5-2. 

Note that the estimation of nuisance parameter ${\left\{g_k\right\}}_{k\in \mathcal{K}}$ can be naturally formulated as the following logistic regression:

$$\begin{array}{c}\hfill \underset{g_k}{min}\widehat{E}\left[{\int }_0^{T\wedge C}\left\{A_{k,t}ln(1+e^{-g_k\left(X_t\right)})+(1-\sum _{\ell \in \mathcal{K}}{A}_{\ell ,t})ln(1+e^{g_k\left(X_t\right)})\right\}dt\right]+{\zeta }_{n,k}{\parallel g_k\parallel }_{{\mathcal{H}}_k}^2,\end{array}$$

(27)

for each $k\in \mathcal{K}$. Here, we have introduced a regularization hyperparameter ${\zeta }_{n,k}>0$.

2.4. Extension to Models with Latent Variables

Suppose that the Assumption 7 is violated due to the presence of unobserved factors affecting the outcome, and covariates cannot adjust the resultant temporal change in the risk set. In this case, as we see in the simulation study below, the inclusion of the time elapsed after enrollment in the set of covariates may improve the description of event occurrence. For such cases, if additional data about the unobserved factors cannot be acquired, the only alternative is to explicitly model the unobserved factors as latent variables.

To simplify the argument and notation, let us assume that the unobserved factor (denoted by $W\in \{0,1\}$) is a single time-independent binary variable (e.g., a variable indicating the presence or absence of a genetic risk factor). And suppose that we use a model with a latent variable $Z\in \{0,1\}$ defined by the following prior distribution of Z and Z-dependent conditional hazard:

$$\begin{array}{c}\hfill Q_{\beta }\left(Z\right|X_0,A_0)\stackrel{\mathrm{def}}{=}{\displaystyle \frac{1}{1+exp(-(2Z-1)({\sum }_{j\in \mathcal{J}}{\beta }_j{X}_{j0}+{\beta }_0))}},\end{array}$$

(28)

and

$$\begin{array}{c}\hfill \tilde{h}\left(t\right|A_t,X_t,Z)\stackrel{\mathrm{def}}{=}exp({\theta }^{\prime }{A}_t+f\left(X_t\right)+\kappa Z).\end{array}$$

(29)

The latter yields the following Z-dependent log-likelihood:

$$\begin{array}{c}\hfill Q_{\tilde{h}}\left(T\right|X,A,Z)=I_{[0,C]}\left(T\right)({\theta }^{\prime }{A}_T+\kappa Z+f\left(X_T\right))-{\int }_0^{T\wedge C}exp({\theta }^{\prime }{A}_t+\kappa Z+f\left(X_t\right))dt.\end{array}$$

(30)

In the above, we have introduced parameters $\beta \in {\mathbf{R}}^{\dim \left(X_0\right)+1}$ relating baseline values of covariates (at $t=0$) to the prior distribution of Z, and a parameter $\kappa \in \mathbf{R}$ representing the effect of the latent variable.

The same argument for Proposition 1 applies to this model. Suppose that the data-generating process is described by

$$\begin{array}{c}\hfill P\left(W\right|X_0,A_0)={\displaystyle \frac{1}{1+exp(-(2W-1)({\sum }_{j\in \mathcal{J}}{\beta }_j^{*}{X}_{j0}+{\beta }_0^{*}))}},\end{array}$$

(31)

and

$$\begin{array}{c}\hfill \underset{\delta t\to +0}{lim}{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a,X_t=x,W=1)}{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a,X_t=x,W=0)}}=const,\end{array}$$

(32)

regardless of the values of x, a and t. Then, replacing X with the combination $(X,W)$ or $(X,Z)$ in Assumptions 1–8 and Proposition 1 and assuming that the model is identifiable, one obtains the same interpretation of hazard ratios.

Then, one can apply Proposition 2 to derive Neyman-orthogonal scores from the following marginal log-likelihood for this model:

$$\begin{array}{c}\hfill \sum _{i\in \mathcal{I}}ln{Q}_{\tilde{h}}\left(T_i\right|A_i,X_i)=\sum _{i\in \mathcal{I}}\underset{r_i}{max}\sum _{Z_i\in \{0,1\}}{r}_i\left(Z_i\right)ln{\displaystyle \frac{Q_{\tilde{h}}\left(T_i\right|X_i,A_i,Z_i)Q_{\beta }\left(Z_i\right|X_i,A_i)}{r_i\left(Z_i\right)}},\end{array}$$

(33)

with the aid of variational distributions ${\left\{r_i\right\}}_{i\in \mathcal{I}}$ constrained by $r_i\left(0\right)+r_i\left(1\right)=1$. It should be noted that, in the application of Proposition 2, f is replaced by $\mathrm{vec}(f,\kappa ,\beta )$. For the above latent-variable model, a more intuitively understandable formula corresponding to Proposition 3 is not available.

2.5. Design of Models with Multiple RKHSs

Although we have discussed a single RKHS for model $\mathcal{M}$, from a practical point of view, it is much more convenient to define f in the following form [29,30]:

$$\begin{array}{c}\hfill f\left(X\right)=\sum _{k_{\ell }\in \mathfrak{K}}{f}_{\ell }\left(X\right)+b,f_{\ell }\in {\mathcal{H}}_{k_{\ell }},b\in \mathbf{R}.\end{array}$$

(34)

In the above, we have used a collection of RKHSs each of which is associated with a kernel function $k_{\ell }\in \mathfrak{K}$. With this formulation, one can separately model the effects of different factors, such as demographics, comorbidities and drug use. We provide examples in simulation studies below. Use of multiple kernels also allow us to construct f by trial and error. One can observe whether incorporation or removal of a component function associated with a covariate improves the data description by calculating BME. On can also validate Assumptions 5–7 by checking whether incorporation of $f_{\ell }$ representing, for example, time inhomogeneity improves BME. Here, note that criteria for model selection must be BME (or Bayesian information criteria applicable to only parametric models), not cross-validation error, because the former has consistency in terms of model selection in the large sample limit, but the latter not.

The theory developed in the previous sections apply to the above multiple-kernel model as well, because the sum of functions in different RKHSs is known to belong to an RKHS associated with a single composite kernel [30,31]. The squared norm of f defined by this composite kernel is dominated by the square sum of the norms of the component functions [31]. Thus, the entire theory is not affected by the use of a multiple-kernel model. The theoretical properties of ML estimation with the model in Equation (34) under a 1-norm, 2-norm or mixed-norm regularization have been extensively studied [30,32,33,34]. We employ the 2-norm-regularized version for which calculation of BME is tractable.

2.6. Estimation Algorithm

In this section, we describe the concrete procedure for computing the debiased estimators of the hazard ratios for treatment. The entire algorithm is summarized in the pseudocode in 1 and is largely divided into three steps. How the algorithm is numerically implemented is detailed in Appendix B.

Algorithm 1: DML of hazard ratios with cross fitting

Input: log-likelihood Equation () or (), partitioned data $D=\{D^{\left(1\right)}\dots D^{\left(M\right)}\}$ Output: optimal model ${\mathcal{M}}^{*}$, debiased estimate $\overline{\theta }$ and its squared standard error $\overline{\Sigma }$ ${\mathcal{M}}^{*}\leftarrow ⌀$ while${\mathcal{M}}^{*}=⌀$do    Choose a model: ${\mathcal{M}}_0\leftarrow {\mathcal{H}}_{k_0}$.    Prepare alternative models: ${\mathcal{M}}_{\ell }\leftarrow {\mathcal{H}}_{k_{\ell }}$$(\ell =1\dots L)$    for $\lambda$, $\sigma \in$ SearchGrid do      Perform ML estimation with D and Compute ${\mathrm{BME}}_{\lambda ,\sigma ,\ell }$ for $(\ell =0\dots L)$    end for    Find ${\lambda }_{\ell }^{*},{\sigma }_{\ell }^{*}\leftarrow argmax{\mathrm{BME}}_{\lambda ,\sigma ,\ell }$ for $(\ell =0\dots L)$    if ${\mathrm{BME}}_{{\lambda }_0^{*},{\sigma }_0^{*},0}\ge {\mathrm{BME}}_{{\lambda }_{\ell }^{*},{\sigma }_{\ell }^{*},\ell }$$(1\le \ell \le L)$ then      ${\mathcal{M}}^{*},{\lambda }^{*},{\sigma }^{*}\leftarrow {\mathcal{H}}_{k_0},{\lambda }_0^{*},{\sigma }_0^{*}$    end if end while for$m=1\dots M$do    $D^{(\setminus m)}\leftarrow D\setminus D^{\left(m\right)}\cup D^{(m+1(\mathrm{mod}M\left)\right)}$    Perform ML estimation with $D^{(\setminus m)}$ to obtain ${\widehat{\theta }}^{(\setminus m)},{\widehat{\theta }}^{(\setminus m)}$ with ${\lambda }^{*}$ and ${\sigma }^{*}$. end for for$\zeta \in$ SearchGrid do    for $m=1\dots M$ do      $D^{(\setminus m)}\leftarrow D\setminus D^{\left(m\right)}\cup D^{(m+1(\mathrm{mod}M\left)\right)}$      Compute ${\widehat{H}}^{(\setminus m)}$ (or $\left\{{\widehat{g}}_k^{(\setminus m)}\right\}$) based on ${\widehat{\theta }}^{(\setminus m)},{\widehat{f}}^{(\setminus m)}$ and $D^{(\setminus m)}$.      Compute ${\mathrm{CVErr}}_{H,\zeta }^{\left(m\right)}$ (or ${\left\{{\mathrm{CVErr}}_{g_k,\zeta }^{\left(m\right)}\right\}}_k$) based on ${\widehat{H}}^{(\setminus m)}$ (or $\left\{{\widehat{g}}_k^{(\setminus m)}\right\}$) and $D^{(m+1\mathrm{mod}M)}$.    end for end for ${\zeta }^{*}\leftarrow arg{max}_{\zeta }{\sum }_{m=1}^M{\mathrm{CVErr}}_{H,\zeta }^{\left(m\right)}$ (or ${\zeta }_k^{*}\leftarrow arg{max}_{\zeta }{\sum }_{m=1}^M{\mathrm{CVErr}}_{g_k}^{\left(m\right)}$). for$m=1\dots M$do    Compute $\varphi (D_i;\theta ,{\widehat{f}}^{(\setminus m)},{\widehat{H}}^{(\setminus m)})$ (or $\varphi (D_i;\theta ,{\widehat{f}}^{(\setminus m)},{\widehat{g}}^{(\setminus m)})$) for $i\in D^{\left(m\right)}$. end for Obtain $\overline{\theta }$ by finding the zero of ${\sum }_{m=1}^M{\sum }_{i\in D^{\left(m\right)}}\varphi (D_i;\theta ,{\widehat{f}}^{(\setminus m)},{\widehat{H}}^{(\setminus m)})$ (or ${\sum }_{m=1}^M{\sum }_{i\in D^{\left(m\right)}}\varphi (D_i;\theta ,{\widehat{f}}^{(\setminus m)},{\widehat{g}}^{(\setminus m)})$). Compute $\overline{\Sigma }$ based on $\overline{\theta }$.

Step 1: Model selection and determination of hyperparameters

We first choose the best possible combination of RKHSs by trial and error by comparing BME in the manner we have described in the previous section. We also examine the validity of Assumptions 5–7 in this step. To compare BME for different models, we first perform a grid search of the optimal hyperparameter values for regularization and the bandwidth of kernels for each model, examining values with intervals of approximately $ln\left(1.5\right)$ on the logarithmic scale (e.g. $2.0,3.0,5.0,7.0,10.0\dots$) for each hyperparameter. Then, the BME values of different models calculated with the optimal hyperparameter values are compared. When we think that we have obtained the best possible model, we proceed to the next step. There is, however, a possibility that the validity of Assumptions 5–7 is denied for the best possible model at hand. In this case, one needs to review the study design, obtain new covariates or incorporate latent variables to the model.

Step 2: Estimation of nuisance parameters

For the model chosen in Step 1, we calculate nuisance parameters. Following Chernozhukov et al. [19], we perform “cross fitting” as follows. To debias $\theta$, the estimation errors in f and H (or $g_k$) must be independent of the data $D_i$ in the argument of score functions. To achieve this, the entire dataset D are first partitioned into M groups of approximately equal size ($D^{\left(1\right)}\dots D^{\left(M\right)})$. We perform the estimation of nuisance parameters M times, using $D_{\mathrm{train}}^{\left(m\right)}\stackrel{\mathrm{def}}{=}D\setminus (D^{\left(m\right)}\cup D^{(m+1(\mathrm{mod}M\left)\right)})$ as a training set, validating the estimation result with $D_{\mathrm{val}}^{\left(m\right)}\stackrel{\mathrm{def}}{=}{D}^{(m+1(\mathrm{mod}M\left)\right)}$ and then yielding the values used for score functions of $D_{\mathrm{hout}}^{\left(m\right)}\stackrel{\mathrm{def}}{=}{D}^{\left(m\right)}$ for $1\le m\le M$.

Concretely, we first perform ML estimation to obtain $\widehat{{\theta }^{(\setminus m)}}$ and $\widehat{f^{(\setminus m)}}$ as the minimizer of the negative log-likelihood of $D_{\mathrm{train}}^{\left(m\right)}$ using the optimal hyperparameter values obtained in Step 1. Then, we proceed to the computation of additional nuisance parameters $\widehat{H}$ (or ${\left\{{\widehat{g}}_k\right\}}_k$) based on $D_{\mathrm{train}}^{\left(m\right)}$. In this step, we also need to identify the value of the regularization hyperparameter ${\zeta }_n$ for these nuisance parameter. For the determination of ${\zeta }_n$ for $\widehat{H}$, we examine the following cross-validation error

$$\begin{array}{c}\hfill {\mathrm{CVErr}}_H\stackrel{\mathrm{def}}{=}\sum _{k\in \mathcal{K}}\sum _{1\le m\le M}{∥{\widehat{H}}_{{\theta }_{k}f}^{(m,\mathrm{val})}-{\widehat{H}}_{{\theta }_{k}f}^{(\setminus m)}{({\widehat{H}}_{ff}^{(\setminus m)}+{\zeta }_n)}^{-1}{\widehat{H}}_{ff}^{(m,\mathrm{val})}∥}_{{\mathcal{H}}_k}^2\end{array}$$

(35)

Here, the ${\widehat{H}}_{\theta f}^{(m,\mathrm{val})}$ and ${\widehat{H}}_{ff}^{(m,\mathrm{val})}$ are the block components of the Hessian of the negative log-likelihood of $D_{\mathrm{val}}^{\left(m\right)}$.

For the determination of ${\left\{{\zeta }_{n,k}\right\}}_k$ for ${\left\{g_k\right\}}_k$, the logarithmic loss in Equation (27) may underestimate the error due to large $e^{g_k}$ and $e^{-g_k}$ in the score functions. We therefore suggest computing the following cross-validation error:

$$\begin{array}{ccc}& & {\mathrm{CVErr}}_{g_k}\stackrel{\mathrm{def}}{=}\hfill \\ & & \sum _{1\le m\le M}{\left[\sum _{i\in D_{\mathrm{val}}^{\left(m\right)}}{\int }_0^{T_i\wedge C_i}\left\{A_{ik,t}(e^{-\widehat{g_k^{(\setminus m)}}\left(X_{i,t}\right)}-1)+(1-\sum _{\ell \in \mathcal{K}}{A}_{i\ell ,t})(e^{\widehat{g_k^{(\setminus m)}}\left(X_{i,t}\right)}-1)\right\}dt\right]}^2.\hfill \end{array}$$

(36)

The above error design exploits the fact that, for $\widehat{P}(A_k=1|X)\to P(A_k=1|X)$,

$$\begin{array}{c}\hfill P(A_k=1|X){\displaystyle \frac{1}{\widehat{P}(A_k=1|X)}}+(1-P(A_k=1|X)){\displaystyle \frac{1}{1-\widehat{P}(A_k=1|X)}}\to 2.\end{array}$$

(37)

Here, the trivial solution $\widehat{P}(A_k=1|X)={\displaystyle \frac{1}{2}}$ also minimizes ${\mathrm{CVErr}}_{g_k}$ and should be avoided by checking BME of the logistic regression (Equation (27)) as well.

Step 3: Debiased estimation of $\theta$ and its standard error.

The debiased estimator $\overline{\theta }$ of ${\theta }^{*}$ is obtained as the zero of the following:

$$\begin{array}{c}\hfill \sum _{1\le m\le M}\sum _{i\in D_{\mathrm{hout}}^{\left(m\right)}}\varphi (D_i;\theta ,{\widehat{\eta }}^{(\setminus m)}),\end{array}$$

(38)

with ${\widehat{\eta }}^{(\setminus m)}=({\widehat{f}}^{(\setminus m)},{\widehat{H}}^{(\setminus m)})$ (or $({\widehat{f}}^{(\setminus m)},{\widehat{g}}_k^{(\setminus m)}$)). According to Chernozhukov et al. [19], the asymptotic standard error of the estimator is given by ${\Sigma }^{-1/2}(\overline{\theta }-{\theta }^{*})\stackrel{d}{\to }\mathrm{Normal}(0,1)$ with

$$\begin{array}{ccc}\hfill \widehat{\Sigma }& \stackrel{\mathrm{def}}{=}& {\widehat{J}}^{-1}\left({\displaystyle \frac{1}{n^2}}\sum _{1\le m\le M}\sum _{i\in D^{\left(m\right)}}\varphi (D_i;\overline{\theta },({\widehat{f}}^{(\setminus m)},{\widehat{\eta }}^{(\setminus m)}))\varphi {(D_i;\overline{\theta },({\widehat{f}}^{(\setminus m)},{\widehat{\eta }}^{(\setminus m)}))}^{\prime }\right){\widehat{J}}^{-1\prime },\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill \widehat{J}& \stackrel{\mathrm{def}}{=}& {\displaystyle \frac{1}{n}}\sum _{1\le m\le M}\sum _{i\in D^{\left(m\right)}}{\displaystyle \frac{\partial }{\partial \theta }}\varphi (D_i;\theta ,({\widehat{f}}^{(\setminus m)},{\widehat{\eta }}^{(\setminus m)})){|}_{\theta =\overline{\theta }}.\hfill \end{array}$$

(40)

Thus, in the simulation study described below, the square roots of the diagonal elements of $\widehat{\Sigma }$ serve to scale the errors from the true value to obtain the t-statistics.

Remark 6. 

The orthogonal scores as functions of $\theta$ are given in Proposition 3 for the case with $g_k$. For the case with H, a procedure for obtaining functions of $\theta$ is given in Appendix B.

3. Results of Numerical Simulations

3.1. Simulation Result 1: Adjustment for Observed Confounders

This section presents a simulation study designed to demonstrate that, in a clinically plausible setting, the proposed method estimates hazard ratios for treatment with minimal bias. Suppose we randomly sample subjects from an electronic medical record that contains the medical history of a large local population. Suppose further that 2,000 subjects aged 50–75 years enrolled at uniformly random times between 2000 and 2005. The chosen subjects are followed up for $C\stackrel{i.i.d.}{\sim }\mathrm{Uniform}[5,10]$ years, with timely recording of the onset of comorbidities, the administration of drugs and the occurrence of a disorder defined as an outcome event.

As a typical setting in which confounders bias the estimation of the treatment effect, suppose that subjects are newly diagnosed with condition 1 at a constant rate (2.5%/month), which may be a prodrome for the outcome event. A drug that is suspected to be a risk factor for an outcome event may be used to treat this condition. In the simulation, the probability per month of initiating this drug is $0.004+0.2X_{1t}exp(-X_{1t}/0.6)/{\left(0.6\right)}^2$ with a covariate $X_{1t}$ [year] indicating the duration of suffering condition 1 before time t, and the probability per month of discontinuing the drug is 0.01. Use of this drug increases the risk of developing condition 2. The probability per month of developing condition 2 at time t is $0.05+0.05{\int }_{-\infty }^t{B}_s\Theta (t-s-1/12)e^{-3(t-s-1/12)}ds$, where $B_s$ takes a value of 1 if the subject is a current user of the drug at time s, or it takes a value of 0 otherwise. The above integral including the Heaviside step function $\Theta$ describes a short-term delayed increase in the risk of developing condition 2 subsequent to the drug use.

Finally, suppose that the occurrence of the outcome event depends on both treatment and conditions 1 and 2, and its true risk per month is $exp(-7.0+{\theta }_1^{*}{A}_{1t}+{\theta }_2^{*}{A}_{2t}+0.04{\mathrm{age}}_t+0.2sin(\pi {\mathrm{date}}_t/8)-0.2cos(\pi {\mathrm{date}}_t/6)+2.0X_{1}exp(-X_{1t}/1.5)+P_2(1-exp(-X_{2t}/2.5)))/12.0$, with ${\theta }_1^{*}=1.0$ and ${\theta }_2^{*}=2.0$, where we define two treatment variables $A_{1t}$ and $A_{2t}$, which take a value of 1 if and only if the subject use the drug for a period $\le 18$ and $>18$ months before time t, respectively. $P_2$ in the equation is a parameter for which we use different values. The covariates ${\mathrm{age}}_t$, ${\mathrm{date}}_t$ and $X_{2t}$ (duration of suffering condition 2) are measured in years.

For this dataset, we performed debiased ML estimation using multiple-kernel models for f. Concretely, we used the following model:

$$\begin{array}{c}\hfill f\left(X_t\right)=f_{\mathrm{age}}\left({\mathrm{age}}_t\right)+f_{\mathrm{date}}\left({\mathrm{date}}_t\right)+f_1\left(X_{1t}\right)+f_2\left(X_{2t}\right)+b,\end{array}$$

(41)

with $f_{\mathrm{age}}\in {\mathcal{H}}_{k_{\mathrm{age}}}$, $f_{\mathrm{date}}\in {\mathcal{H}}_{k_{\mathrm{date}}}$, $f_1\in {\mathcal{H}}_{k_1}$ and $f_2\in {\mathcal{H}}_{k_2}$. To correctly specify the model, we chose $k_{\mathrm{age}}$ to be a linear kernel and $k_{\mathrm{date}}$, $k_1$ and $k_2$ to be one-dimensional Gaussian kernels [e.g. $k_1(X_{i_1,t_1},X_{i_2,t_2})=exp(-{(X_{i_{1}1,t_1}-X_{i_{2}1,t_2})}^2/{\sigma }^2)$] with bandwidth hyperparameter, $\sigma$. Although we can a priori identify the set of necessary kernels, one can find the appropriate set of kernels by trial and error, using BME as a criterion. In this step, we can also validate the model assumptions (Assumption 5–7) by examining whether the introduction of an additional kernel improves BME. We briefly illustrate the idea through an example.

First, as we have described in the above, covariate $X_{2t}$ is necessary for predicting the onset of outcome events. This can be seen by comparing BME values for the correctly specified model and a misspecified model built by removing $f_2\left(X_{2t}\right)$ and ${\mathcal{H}}_{k_2}$ from it (Figure 1(A)). The question naturally arising here is how one can see that there is a room for improvement when they have the misspecified model at hand. We propose to examine the assumption of time homogeneity (Assumption 7) as an indicator. In our example, comparison of BME values for the second, misspecified model with or without inclusion of the time elapsed after enrollment as a covariate shows the violation of the assumption, depending on the value of $P_2$ (Figure 1(B)), although the sensitivity is not high. The risk variation due to the time elapsed after enrollment indicates a temporal change in the risk set that covariates at hand do not account for, and therefore the model is missing some necessary factors (see the argument in Remark 3-1). For a correctly specified model, such violation of the assumption is not detected (Figure 1(C)). The violation of the time-homogeneity assumption is observed for models misspecified in different manners (which we omit to avoid redundancy).

Next, suppose that we successfully identified the correctly specified model and proceeded to the calculation of debiased estimators of hazard ratios. We carried out this calculation in two different ways based on either one of Propositions 2 and 3. In both applications, we followed the procedure described in Section 2.6. Additionally, in the application of Proposition 3, we use the following multiple-kernel model for the estimation of $g_{\ell }$ ($\ell \in \mathcal{K}$):

$$\begin{array}{c}\hfill g_{\ell }\left(X_t\right)=g_{\ell ,\mathrm{age}}\left({\mathrm{age}}_t\right)+g_{\ell ,\mathrm{date}}\left({\mathrm{date}}_t\right)+g_{\ell 1}\left(X_{1t}\right)+g_{\ell 2}\left(X_{2t}\right)+b,\end{array}$$

(42)

with $g_{\ell ,\mathrm{age}}\in {\mathcal{H}}_{k_{\mathrm{age}}}$, $g_{\ell \mathrm{date}}\in {\mathcal{H}}_{k_{\mathrm{date}}}$, $g_{\ell 1}\in {\mathcal{H}}_{k_1}$ and $g_{\ell 2}\in {\mathcal{H}}_{k_2}$, where all RKHSs are associated with a one-dimensional Gaussian kernel. On top of the kernels described above, we also examined whether the inclusion of two-dimensional Gaussian kernels such as $k_{12}(X_{i_1,t_1},X_{i_2,t_2})=exp(-({(X_{i_{1}1,t_1}-X_{i_{2}1,t_2})}^2+{(X_{i_{1}2,t_1}-X_{i_{2}2,t_2})}^2)/{\sigma }_2^2),$ improves BME for the logistic regression in Equation (27). In the present case, the two-dimensional kernels did not improve BME.

In the tuning of hyperparameter values for nuisance parameters H (or $g_{\ell }$), we examined how the hyperparameters affect ${\mathrm{CVErr}}_H$ (Figure 1(D)) and ${\mathrm{CVErr}}_{g_{\ell }}$ (Figure 1(F)), respectively. In the application of Proposition 3, BME for the logistic regression and ${\mathrm{CVErr}}_{g_{\ell }}$ identified different optimal values (Figure 1(E) and (F)).

After taking all of the steps in Algorithm 1, we calculated debiased estimate $\overline{\theta }$ using the score functions obtained with either $(H_{\theta f},H_{ff})$ or $g_1$. To measure the bias in the estimator, we obtained estimates with 500 datasets generated from the identical process described above with different seeds for the random-number generator. The t-statistics ($\stackrel{\mathrm{def}}{=}({\widehat{\theta }}_k-{\theta }_k^{*})/\widehat{{\mathrm{SE}}_{{\theta }_k}}$) in Figure 1(G) and (H) shows that the debiased estimator identified the true value with minimal bias, in contrast with the naive ML estimator. In Figure 1(H), we also observed that $g_1$ tuned with ${\mathrm{CVErr}}_{g_1}$ debias the estimator more effectively than $g_1$ tuned with BME.

3.2. Simulation Result 2: Estimation of Treatment Effect in Population with Heterogeneous Risk

Suppose that most of the simulation settings are the same as the previous one, but the enrolled subjects are divided into two distinct risk groups denoted by $W\in \{0,1\}$, and this is reflected in a baseline covariate value. More precisely, the number of subjects, their age, enrollment date, censoring, onset probabilities of conditions 1 and 2 and treatment assignment are the same as the previous one. Suppose, however, that blood-test results $X_{j0}^{\left(\mathrm{test}\right)}\stackrel{i.i.d.}{\sim }\mathrm{Normal}(0,{\sigma }_j^{\left(\mathrm{test}\right)2})$ at baseline ($t=0$) are available ($1\le j\le 3$), and the risk group to which the subject belongs is probabilistically related to ${\left\{X_{j0}^{\left(\mathrm{test}\right)}\right\}}_j$ as

$$\begin{array}{c}\hfill P\left(W\right|A,X)={\displaystyle \frac{1}{1+exp(-(2W-1)({\sum }_j{\beta }_j^{*\prime }{X}_{j0}^{\left(\mathrm{test}\right)}+{\beta }_0^{*}))}},\end{array}$$

(43)

with ${\left\{{\sigma }_j^{\left(\mathrm{test}\right)}\right\}}_{1\le j\le 3}=\{2.0,1.0,4.0\}$ and ${\left\{{\beta }_j^{*}\right\}}_{0\le j\le 3}=\{-0.5,1.0,0.0,0.0\}$. With variable W, the risk of outcome event per month is given by $exp(-7.0+{\theta }_1^{*}{A}_{1t}+{\theta }_2^{*}{A}_{2t}+0.04{\mathrm{age}}_t+0.2sin(\pi {\mathrm{date}}_t/8)-0.2cos(\pi {\mathrm{date}}_t/6)+2.0X_{1}exp(-X_{1t}/1.5)+P_2(1-exp(-X_{2t}/2.5))+{\kappa }^{*}W)/12.0$ with ${\theta }_1^{*}=1.0$, ${\theta }_2^{*}=2.0$, $P_2=0.5$ and ${\kappa }^{*}=1.0,2.0$ or $3.0$.

We first performed analysis with only observed covariates $(X_t,X_t^{\left(\mathrm{test}\right)})$, not using latent variables. Here, we assumed $X_t^{\left(\mathrm{test}\right)}=X_0^{\left(\mathrm{test}\right)}$. In the validation of assumptions, the inclusion of the time elapsed after enrollment in the set of covariates improved BME (Figure 2(A)), suggesting that the enrolled subjects are heterogeneous, presumably due to an unobserved factor. Then, we analyzed the dataset with a latent variable, maximizing the marginal log-likelihood with the EM algorithm (see e.g., Chapter 9 of Ref.[26]), which led to further improvement of BME (Figure 2(B)). In this calculation, we observed that the optimization becomes numerically unstable for small $\kappa (\le 1.0)$ presumably because of the singularity of the model. We examined the estimated values of $\theta$ with 500 datasets generated by the identical process with different seeds for its random-number generator, and the results confirmed that the suggested debiased estimators identify the true value with minimal bias, in contrast with the naive ML estimator (Figure 2(C)).

4. Discussion

Preceding studies [3,4,5] showed that the Cox model allows multiple, contradictory interpretations due to the selection process in the risk set. To alleviate this difficulty, we proposed a new framework that combines an exponential hazard model with machine learning models. With the observation that the baseline hazard of the Cox model is a major source of the uninterpretability, we constructed the model without a baseline hazard. For this model, with a few testable assumptions as well as conventional assumptions for causal inference in observational studies, we clarified the context in which the estimated hazard ratios are causally interpreted. In this argument we have demonstrated that the source of uninterpretability is mostly removed (Remark 3-1). Then, we developed a framework to systematically seek for the best possible model with the aid of machine learning and then compute debiased estimates of hazard ratios with the obtained model. Numerical simulations demonstrated that wrong description of the selection process in the risk set can be detected as violation of a time-homogeneity assumption in this framework. In the simulation, it was also demonstrated that our theoretically justified debiased ML estimators of hazard ratios identify the true values with minimal bias.

Epidemiologists and statisticians might find it psychologically difficult to abandon the Cox models with baseline hazard because of its established role as a default model. The popularity of the Cox model is due to its semiparametric nature. Leaving the aspects of data difficult to model as a blackbox is an attractive idea. Furthermore, the Cox’s MPL estimator [35] for log-hazard ratios of treatment groups is asymptotically efficient [36,37], which is a natural consequence of obviating the estimation of the baseline hazard. Because of this asymptotic efficiency, all other estimators, including the ML estimator, have essentially been reduced to theoretical subjects, while a relatively small number of studies have suggested the advantage of using the others (e.g. [21,38]). The present results suggest the need to reconsider this situation by reevaluating the tradeoff between the efficiency and uninterpretability due to the unspecified baseline hazard. The disadvantage of using MPL is not only the uninterpretability. First of all, the baseline hazard is impossible to interpret from the physical and biological points of view. No physical process depends on the time elapsed after the enrollment. Thus, its use is solely justified by convenience. However, from a technical point of view, comparing the estimation of $g_k$ in the present study and the inverse-probability-of-treatment weight for marginal structural Cox models used for a similar purpose [17], one can see that the former is a marginal estimand in terms of time and hence usually easier to estimate. The standard method for the latter, that is, estimating the transition probability of treatment and censoring for each value of covariates and each time t, possibly works better in a simple model setting, such as analysis involving only a few–several discrete variables. However, carrying this out under the nonlinear effect of a high-dimensional covariate is prohibitive. In fact, the recently proposed application of doubly robust estimation to marginal structural Cox models are still limited to simple model settings to which machine learning cannot be incorporated [39]. This is why we could not directly compare the results with our method and those with the marginal structural Cox model. Another negative factor for MPL approach is that the combination of Bayes theorem and partial likelihood is only approximately justifiable, so extending Cox regression to Bayesian settings is not straightforward, although efforts in this direction have been made for Cox models as well [40]. Considering these negative factors for the MPL approach, we conclude that the debiased ML approach should be reconsidered as an option in causal inference for observational studies with uncontrolled dynamic treatment and real-time measurement of a large number of covariates.

It is, however, fair to note that the applicability of our approach at the current stage is largely limited because Assumptions 1–8 must be satisfied for the application. Among these assumptions, Assumptions 1–4 are technical assumptions that are explicitly or implicitly made in the conventional approach with the marginal structural Cox model [17]. In particular, in causal inference with observational data, Assumption 2 stating conditionally randomized assignment of treatment for each value of covariate is thought to be crucial. If this does not hold, an experiment (intervention) is needed. See Remarks 1-1–1-4 for the reasoning of the other assumptions. Assumptions 5–8 essentially concern specification of a model. Assumption 6 is potentially removable. Estimation of heterogeneous treatment effect (i.e., covariate-dependent treatment effect) is a subject of active research [41] and an approach that adjusts the heterogeneity in the treatment effect in an unbiased manner may be incorporated into our approach. Then, the remaining problem is whether one can design a suitable set of variables for X and A (Assumption 5) and a correctly specified machine-learning model (Assumption 8) without sacrificing the assumption of time-homogeneity (Assumption 7). At the current stage, the relatively simple model we have used in this study is likely to suffer misspecification in real-world epidemiological problems. Therefore, our proposals stated in this study should be considered as a suggestion of future efforts to overcome this misspecification. Taking the current success of deep learning in real-world applications into account, we would like to suggest an optimistic view that models and data of a large scale in combination with Bayesian inference may remove any such misspecification problem. Before such a goal is achieved, no one knows a priori which of our approach and the Cox’s MPL approach yields a value closer to the true hazard ratio for a given dataset.

In technical aspects, we indeed restricted our study to a simple setting for the sake of clear exposition. In this study, only non-informative right censoring is considered, only a single binary treatment variable can take the value of one, and only a single time-independent latent variable is allowed. However, the advantage of ML estimation has been demonstrated in the context of informative censoring [42]. Therefore, extending our framework in this direction seems promising. The use of multiple independently changing treatment variables is used for the sake of an exposition, and our theories apply to general cases. Application of our framework to models with continuous treatment variables is also straightforward. Proposition 2 applies to models with linear functions of continuous treatment variables, while Proposition 3 cannot be used for this case. Since our theory and the framework of DML rely on the low dimensionality of treatment variables, extension to models with nonlinear treatment effect of continuous variables is beyond the fundamental limit of the current technique. However, use of multiple binary variables can approximately describe nonlinear dependence on treatment variables. This is the reason why we chose the binary setting. The extension of our framework to models with multiple (possibly continuous, time-dependent) latent variables would face the identifiability problem. Although progress has been made in both theoretical analysis of identifiability [43,44,45,46] and numerical examination of identifiability in the context of biological modeling [47], there is still a technical challenge. Related to this topic, using a singular model can lead to effective breakdown of asymptotic normality of estimators [48]. In this case, the framework of debiased machine learning would itself need modifications, and the Laplace approximation of BME also fails and needs to be replaced by Monte-Carlo integration [49] (or information criteria for singular parametric models [50,51]). The above difficulties being mentioned, use of rich classes of latent-variable models for which inference methods have been extensively developed in machine learning (especially those for time-series data [52,53]) is yet an attractive, worthy challenge. In this challenge, our technique for debiasing will be useful.