Debiased Maximum-Likelihood Estimators of Hazard Ratios Under 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 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, as suggested by semiparametric theories for causal inference (see, e.g., Ref.[16]). 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 [17]. The marginal structural Cox model potentially used for this purpose [18] suffers from the uninterpretability problem described above, and thus, development of an alternative method is required.

The present study proposes a strategy to model heterogeneous risk groups in a study population by using an exponential hazard model based on time-dependent treatment variables and covariates, but with no indefinite baseline hazard. Although Martinussen et al. (2020) did not explicitly describe the role of the baseline hazard in the problem about interpretability, it was the indefinite baseline hazard that allows the Cox model to be correctly specified for arbitrary temporal change in the risk set, namely, what they call “selection”. Thus, we abandon the semiparametric nature of the Cox models, 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 [19,20]. 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.[21] for an introductory review). Under mostly testable assumptions, the estimated hazard ratios can be interpreted as a measure of causal treatment effect. We show that this approach can be systematically extended to settings with latent variables. In simulation studies with clinically plausible settings, the results confirm that the proposed method identifies the ground truth with minimal bias. In these results, we also show that the proposed ML approach has multiple advantages over Cox’s maximum-partial-likelihood (MPL) approach.

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) 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}$, 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})=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 {\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^{min(T_i,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.[22]), 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 Eq.(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 Eq.(3) results in the summand of Eq.(2) up to a constant independent of $\theta$ and f. Here, it is reasonable to consider that the two assumptions hold almost surely in the $\delta t\to +0$ limit.

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 $T^a=T|T\le t,A_{(-\infty ,t]}=a_{(-\infty ,t]}$.

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): 

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 ground-truth (conditional) probability of the argument variables. 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: 

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)

holds, 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. A similar assumption is conventionally made in the analysis of observational data with a marginal structural Cox model [18].

In addition to Assumptions 1–4, we assume that the model in Eq.(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{M}\left(\mathcal{X}\right)}{sup}\mathrm{E}\left[ln{Q}_{h_{\theta ,f}}\left(T\right|A,X)\right],\end{array}$$

(10)

where $\mathfrak{M}\left(\mathcal{X}\right)$ 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.[23]). 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 [18]. 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. 

We have,

$$\begin{array}{ccc}& \underset{\delta t\to +0}{lim}& {\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 \\ & =& \underset{\delta t\to +0}{lim}{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T^a\right)=1|T^a,T^{a^{\prime }}>t,A_{(-\infty ,t+\delta t]}=a_{(-\infty ,t+\delta 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+\delta t]}=a_{(-\infty ,t+\delta t]}^{\prime },X_{(-\infty ,t+\delta t]}=x_{(-\infty ,t+\delta t]})}}\hfill \\ & =& \underset{\delta t\to +0}{lim}{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_{(-\infty ,t+\delta t]}=a_{(-\infty ,t+\delta t]},X_{(-\infty ,t+\delta t]}=x_{(-\infty ,t+\delta t]})}{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_{(-\infty ,t+\delta t]}=a_{(-\infty ,t+\delta t]}^{\prime },X_{(-\infty ,t+\delta t]}=x_{(-\infty ,t+\delta t]})}}\hfill \\ & =& \underset{\delta t\to +0}{lim}{\displaystyle \frac{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a_t,X_t=x_t)}{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=a_t^{\prime },X_t=x_t)}}.\hfill \end{array}$$

(13)

The first and second equalities are due to no unmeasured confounder and consistency (Assumptions 2 and 1), respectively. The last equality is due to Assumptions 4 and 5. Identifying the last line with the ML estimator is straightforward as follows. Using Assumptions 4 and 5 and Eq.(3), the expected likelihood is represented as

$$\begin{array}{ccc}& & \mathrm{E}\left[Q_h\left(T\right|X,A)\right]\hfill \\ & =& {\int }_0^{\infty }{\int }_0^C{\int }_{\mathcal{X}\times \mathcal{A}}\underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}\left\{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t,X_t)({\theta }^{\prime }{A}_t+f\left(X_t\right))\right.\hfill \\ & & \left.-exp({\theta }^{\prime }{A}_t+f\left(X_t\right))\delta t\right\}dP(X_t,A_t|T>t)P(T>t)dtdP\left(C\right)+const\hfill \\ & =& \sum _{k\in \mathcal{K}}{\int }_0^{\infty }{\int }_0^C{\int }_{\mathcal{X}\times \mathcal{A}}\underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}\left\{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_{k,t}=1,X_t)({\theta }_k+f\left(X_t\right))\right.\hfill \\ & & \left.-exp({\theta }_k+f\left(X_t\right))\delta t\right\}dP\left(X_t\right|T>t,A_{k,t}=1)P(T>t,A_{k,t}=1)dtdP\left(C\right)\hfill \\ & & +{\int }_0^{\infty }{\int }_0^C{\int }_{\mathcal{X}\times \mathcal{A}}\underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}\left\{P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t=0,X_t)f\left(X_t\right)\right.\hfill \\ & & \left.-exp\left(f\left(X_t\right)\right)\delta t\right\}dP\left(X_t\right|T>t,A_t=0)P(T>t,A_t=0)dtdP\left(C\right)+const.\hfill \end{array}$$

(14)

Here, terms independent of $\theta$ and f have been regarded as constants, and the integration with respect to the ground-truth probability measure has been denoted by $dP$. Differentiating this representation gives

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial }{\partial {\theta }_k}}\mathrm{E}\left[Q_h\left(T\right|X,A)\right]& =& {\int }_0^{\infty }{\int }_0^C{\int }_{\mathcal{X}\times \mathcal{A}}\left\{\underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_{k,t}=1,X_t)\right.\hfill \\ & & \left.-exp({\theta }_k+f\left(X_t\right))\right\}dP\left(X_t\right|T>t,A_{k,t}=1)P(T>t,A_{k,t}=1)dtdP\left(C\right),\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & \underset{ϵ\to +0}{lim}{\displaystyle \frac{1}{ϵ}}\mathrm{E}[Q_{h_{\theta ,f+ϵ\psi }}\left(T\right|X,A)-Q_{h_{\theta ,f}}\left(T\right|X,A)]\hfill \\ & =& {\int }_0^{\infty }{\int }_0^C{\int }_{\mathcal{X}\times \mathcal{A}}\left\{\underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}P(I_{(t,t+\delta t]}\left(T\right)=1|T>t,A_t,X_t)\psi \left(X_t\right)\right.\hfill \\ & & \left.-\psi \left(X_t\right)exp({\theta }^{\prime }{A}_t+f\left(X_t\right))\right\}dP(X_t,A_t|T>t)P(T>t)dtdP\left(C\right).\hfill \end{array}$$

(16)

With Assumptions 6 and 7, it is seen that these derivatives vanish if and only if

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

(17)

hold almost everywhere in $\mathcal{X}\times [0,\infty )$ with respect to $dP(X_t,A_t|T>t)P(T>t)P(C>t)$, and the concavity of the expected log-likelihood ensures that this solution is the maximizer. □

Remark 3-1. 

Suppose that Assumption 2 holds for arbitrary $ϵ>0$. Then, for arbitrary $t>t^{\prime }$, we obtain

$$\begin{array}{c}\hfill \underset{\delta t\to +0}{lim}{\displaystyle \frac{1}{\delta t}}P(I_{(t,t+\delta t]}\left(T^a\right)=1|T^a>t,X_{(-\infty ,t+\delta t]},A_{(-\infty ,t^{\prime })}=a_{(-\infty ,t^{\prime })})\\ \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_t=a_t,X_t),\end{array}$$

(18)

with the aid of Assumptions 1–5 in the same manner as we have done for Eqs.(6) and (13). If we can assume $A_{(-\infty ,t^{\prime })}=a_{(-\infty ,t^{\prime })}=0$ with probability one at some $t^{\prime }$, the above equation implies that the counterfactual survival function is given by

$$\begin{array}{c}\hfill S^a\left(t\right)=\mathrm{E}\left[exp\left(-{\int }_0^{t}exp({\theta }^{*\prime }{a}_s+f^{*}\left(X_s\right))ds\right)\right].\end{array}$$

(19)

Although measures of causal effects can be calculated from this relation, doing so requires an unbiased estimation of $f^{*}$, which is possible, only if a simple model can be used to estimate f without modeling error.

Remark 3-2. 

The quantity in Eq.(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.[18], this corresponds to the treatment effect measured in the next month (for which covariates at each month can be regarded as baseline covariates [24]).

Remark 3-3. 

As discussed in Martinussen et al. (2022) [6], of more interest is the quantity

$$\begin{array}{c}\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,X_{(-\infty ,t+\delta t]}^a=X_{(-\infty ,t+\delta t]}^{a^{\prime }})}{P(I_{(t,t+\delta t]}\left(T^{a^{\prime }}\right)=1|T^a,T^{a^{\prime }}>t,X_{(-\infty ,t+\delta t]}^a=X_{(-\infty ,t+\delta t]}^{a^{\prime }})}},\end{array}$$

(20)

for two arbitrary counterfactual treatment schedules that branch before t. A stronger untestable assumption is required for this to be related to the hazard ratios. For example,

$$\begin{array}{ccc}\hfill T\perp \perp T^a,X^a& |& T>t,{\left\{A_s\right\}}_{s\ge t},{\left\{X_s\right\}}_{s\ge t}\hfill \\ \hfill T^a\perp \perp A,T,X& |& T^a>t,{\left\{X_s^a\right\}}_{s\ge t},\hfill \\ \hfill T^a\perp \perp T^{a^{\prime }},X^{a^{\prime }}& |& T^a>t,{\left\{X_s^a\right\}}_{s\ge t},\hfill \end{array}$$

(21)

for any t and any counterfactual treatment schedules $a,a^{\prime }$. This condition is satisfied if no unmeasured factors affect the outcome, as described by the causal graph in Figure 1(A). If unmeasured factors affect the outcome, as described by the causal graph in Figure 1(B), the above conditions are not satisfied. Although the assumption of time-homogeneity (Assumption 7) is strong and does not allow temporal change in the heterogeneity of risk to be left unmodeled, it does not completely exclude the existence of unmeasured factors affecting outcome. For example, suppose that the survival of subjects for an arbitrary counterfactual treatment schedule a is described by

$$\begin{array}{c}\hfill P(T_i^a>t)=exp\left(-{\int }_0^{t}exp({\theta }^{\prime *}{a}_s+f^{*}\left(X_{i,s}\right))(1+W_{i,s})ds\right),X_{i,s}=X_{i,s}^a,W_{i,s}=W_{i,s}^a,\end{array}$$

(22)

where $W_{i,t}={\sum }_j\delta (t-t_{ij})$ is the so called “shot noise” described by Dirac delta functions centered at time points determined by a homogeneous Poisson point process ${\left\{t_{ij}\right\}}_j$. The subject-wise random effect of $W_i$ shared between $T,T^a$ and $T^{a^{\prime }}$ then violates Eq.(21), whereas $W_i$ at different time points are independent and do not influence the time-homogeneity.

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 Eq.(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 ,\sqrt{n}\parallel \widehat{f}-f^{*}\parallel \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) [20]. 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 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{H}\to {\mathbf{R}}^d$ for i.i.d. data $D_i\in \mathcal{D}$:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{n}}\sum _i\varphi (D_i;\widehat{\theta },\widehat{\eta })=0,\end{array}$$

(23)

where the third argument in the summands is a (possibly infinite-dimensional) nuisance parameter and an estimator $\widehat{\eta }$ of an optimal value ${\eta }^{*}\in \mathcal{H}$ is used. Then, ϕ is said to be “Neyman orthogonal,” if it satisfies

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

(24)

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

(25)

where, in the second line, ${\partial }_{\eta }(·)[\eta -{\eta }^{*}]$ denotes the Gateaux derivative operator [25].

Given a few conditions on the regularity of score functions and the convergence of nuisance parameters, Chernozhukov et al. [20] 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 Eq.(2).

Proposition 2 (Neyman orthogonal scores derived from log-likelihood). 

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

$$\begin{array}{c}\hfill {\mu }_a\left(E\right)={\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}$$

(26)

($E\in \mathfrak{M}\left(\mathcal{X}\right)$) 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}$$

(27)

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^{\min (T,C)}{A}_{k,t}{e}^{f\left(X_t\right)}(1+e^{-g_k\left(X_t\right)})dt\hfill \\ & & +{\int }_0^{\min (T,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}$$

(28)

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

Proof. 

Consider the following hazard model that includes the heterogeneous treatment effect:

$$\begin{array}{c}\hfill {\tilde{h}}_{\theta ,f,f_{\mathrm{het}}}\left(t\right|A_t,X_t)=exp({\theta }^{\prime }{A}_t+f\left(X_t\right)+f_{\mathrm{het}}\left(X_t\right)A_{k,t}).\end{array}$$

(29)

Because we assumed the homogeneity of the treatment effect (Assumption 6), we have, for arbitrary $\psi$,

$$\begin{array}{ccc}& \underset{ϵ\to +0}{lim}& {\displaystyle \frac{1}{ϵ}}\mathrm{E}\left[ln{Q}_{{\tilde{h}}_{{\theta }^{*},f^{*},ϵ\psi }}\left(T\right|A,X)-ln{Q}_{{\tilde{h}}_{{\theta }^{*},f^{*},0}}\left(T\right|A,X)\right]\hfill \\ & & =\mathrm{E}\left[I_{[0,C]}\left(T\right)A_{k,T}\psi \left(X_T\right)-{\int }_0^{min(T,C)}{A}_{k,t}\psi \left(X_t\right)exp({\theta }_k^{*}{A}_{k,t}+f^{*}\left(X_t\right))dt\right]\hfill \\ & & =0.\hfill \end{array}$$

(30)

Let the first line on the right-hand side of Eq.(28) be ${\varphi }_{k,1}(D;\theta ,(f,g_k))$ and the second line be ${\varphi }_{k,2}(D;\theta ,(f,g_k))$. Note that $\mathrm{E}\left[{\varphi }_{k,1}\right]$ with ${\theta }_k={\theta }_k^{*}$ and $f=f^{*}$ coincides with the above equation with $\psi \left(X_t\right)=e^{-{\theta }_k^{*}}(1+e^{-g_k\left(X_t\right)})$. Similarly, for arbitrary $\psi$, we have

$$\begin{array}{ccc}& \underset{ϵ\to +0}{lim}& {\displaystyle \frac{1}{ϵ}}\mathrm{E}\left[ln{Q}_{h_{{\theta }^{*},f^{*}+ϵ\psi }}\left(T\right|A,X)-ln{Q}_{h_{{\theta }^{*},f^{*}}}\left(T\right|A,X)\right]\hfill \\ & & =\mathrm{E}\left[I_{[0,C]}\left(T\right)\psi \left(X_T\right)-{\int }_0^{min(T,C)}\psi \left(X_t\right)exp({\theta }^{*\prime }{A}_t+f^{*}\left(X_t\right))dt\right]\hfill \\ & & =0.\hfill \end{array}$$

(31)

Subtracting the sum of Eq.(30) for all k from Eq.(31) for $\psi \left(X_t\right)=1+e^{g_k\left(X_t\right)}$, we obtain $\mathrm{E}\left[{\varphi }_{k,2}(D;{\theta }^{*},(f^{*},g_k))\right]=0$. Hence, $E\left[{\varphi }_k(D;{\theta }^{*},(f^{*},g_k))\right]=0$ has been proven.

Next, we examine the Gateaux derivatives of the expected scores with respect to the nuisance parameters $(f,g_k)$. First, the derivative of $\mathrm{E}\left[{\varphi }_k\right]$ with respect to $g_k$ at ${\theta }_k={\theta }_k^{*}$ and $f=f^{*}$ is

$$\begin{array}{ccc}& \underset{ϵ\to +0}{lim}& {\displaystyle \frac{1}{ϵ}}\mathrm{E}\left[{\varphi }_k(D;{\theta }^{*},(f^{*},g_k+ϵ\psi ))-{\varphi }_k(D;{\theta }^{*},(f^{*},g_k))\right]\hfill \\ & =& \mathrm{E}\left[-I_{[0,C]}\left(T\right)A_{k,T}\psi \left(X_T\right)e^{-{\theta }_k^{*}-g_k\left(X_T\right)}+{\int }_0^{\min (T,C)}{A}_{k,t}\psi \left(X_t\right)e^{f^{*}\left(X_t\right)-g_k\left(X_t\right)}dt\right]\hfill \\ & & +\mathrm{E}\left[{\int }_0^{\min (T,C)}(1-\sum _{\ell \in \mathcal{K}}{A}_{\ell ,t})\psi \left(X_t\right)e^{f^{*}\left(X_t\right)+g_k\left(X_t\right)}dt-(1-\sum _{\ell \in \mathcal{K}}{A}_{\ell ,T})\psi \left(X_T\right)e^{g_k\left(X_T\right)}dt\right]\hfill \\ & =& 0.\hfill \end{array}$$

(32)

The combination of Eqs.(30) and (31) again proves the last equality. The derivative with respect to f at $\theta ={\theta }^{*}$, $f=f^{*}$ and $g_k=g_k^{*}$ is given by

$$\begin{array}{ccc}& \underset{ϵ\to +0}{lim}& {\left.{\displaystyle \frac{1}{ϵ}}\mathrm{E}\left[{\varphi }_k(D;{\theta }^{*},(f^{*}+ϵ\psi ,g_k^{*})-{\varphi }_k(D;{\theta }^{*},(f^{*},g_k^{*}))\right]\right|}_{f=f^{*}}\hfill \\ & =& \mathrm{E}\left[{\int }_0^{\min (T,C)}(1-\sum _{\ell \in \mathcal{K}}{A}_{\ell ,t})\psi \left(X_t\right)e^{f^{*}\left(X_t\right)}(1+e^{g_k^{*}\left(X_t\right)})-A_{k,t}\psi \left(X_t\right)e^{f^{*}\left(X_t\right)}(1+e^{-g_k^{*}\left(X_t\right)})dt\right]\hfill \\ & =& 0.\hfill \end{array}$$

(33)

Eq.(27) shows that the expected values of the two integrands in the second line cancel out, which proves the last equality. □

Remark 4-1. 

In the above proof, Assumption 8 assures that the argument of the Gateaux derivative $\psi$ does not need to belong to $\mathcal{M}$. As one might notice, the crux of the construction lies in balancing the second and third terms of Eq.(28) 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 4-2. 

Chernozhukov et al. (2018) [20] present a systematic method to derive orthogonal scores from likelihood [see Eq.(2.29) of Ref.[20]. The straightforward application of their method results in more complex scores, each of which contains all of ${\left\{{\theta }_k\right\}}_{k\in \mathcal{K}}$. As we see in this fact, multiple methods for orthogonalization are available. Prioritizing the simple explanation, we propose using the above.

2.4. Estimation Algorithm

• In this section, we describe a concrete procedure for computing the debiased estimators for the hazard ratios for treatment.
• Step 1: Estimation of f and $g_k$.
• Plug-in estimators for f and $g_k$ are used to debias the estimation of $\theta$, following a procedure called “cross-fitting”. To debias $\theta$, the estimation errors in f and $g_k$ must be independent of the data $D_i$ in the argument of score functions. Thus, the entire dataset D are first partitioned into M groups of approximately equal size ($D^{\left(1\right)}\dots D^{\left(M\right)})$. Using $D_{\mathrm{val}}^{\left(m\right)}=D^{\left(m\right)}$, $D_{\mathrm{hout}}^{\left(m\right)}=D^{(m+1(modM\left)\right)}$ and $D_{\mathrm{train}}^{\left(m\right)}=D\setminus (D_{\mathrm{val}}^{\left(m\right)}\cup D_{\mathrm{hout}}^{\left(m\right)})$, $\widehat{f^{(\setminus m)}}$ is estimated as the maximizer of the log-likelihood of $D_{\mathrm{train}}^{\left(m\right)}$ with a suitable regularizer. Similarly, $\widehat{g_k^{(\setminus m)}}$ is estimated as a regularized solution for the following logistic regression:

  $$\begin{array}{c}\hfill \underset{g_k^{(\setminus m)}}{min}\sum _{i\in D_{\mathrm{train}}^{\left(m\right)}}{\int }_0^{min(T_i,C_i)}\left\{A_{ik,t}ln(1+e^{-g_k^{(\setminus m)}\left(X_{i,t}\right)})+(1-\sum _{\ell \in \mathcal{K}}{A}_{i\ell ,t})ln(1+e^{g_k^{(\setminus m)}\left(X_{i,t}\right)})\right\}dt\\ \hfill +\left(\mathrm{regularization}\mathrm{term}\right).\end{array}$$

  (34)

  To determine the regularization parameters for estimating f, we use BME averaged over the M folds. Because the logarithmic loss in Eq.(34) may underestimate the error due to large $e^{g_k}$ and $e^{-g_k}$ in the score functions, we suggest computing the following cross-validation error in Eq.(33) which $g_k$ is targeting:

  $$\begin{array}{ccc}& & {\mathrm{CVErr}}_{g_k}=\hfill \\ & & \sum _{1\le m\le M}{\left[\sum _{i\in D_{\mathrm{val}}^{\left(m\right)}}\left\{{\int }_0^{min(T_i,C_i)}{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}$$

  (35)

  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}$$

  (36)

  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 as well.
• Step 2: Validation of hypotheses. The assumptions for Proposition 1 should be validated. Validating Assumptions 5–7 can be done by carrying out the ML estimation with replacement of $A_t$ by $(A_t,A_{t-\tau })$, $X_t$ by $(X_t,X_{t-\tau })$, $f\left(X\right)$ by $f_1\left(X\right)+f_2(X,t)$ and ${\theta }_{\ell }{A}_{\ell }$ by $A_{\ell }({\theta }_{\ell }+f_{\ell }\left(X\right))$, respectively, and comparing the resultant BME (i.e., the posterior probabilities of the null and alternative hypotheses).
• Step 3: Debiased estimation of $\theta$ and its standard error.
• The estimator $\widehat{\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^{\left(m\right)}}\varphi (D_i;\theta ,({\widehat{f}}^{(\setminus m)},{\widehat{g}}^{(\setminus m)})).\end{array}$$

  (37)

  The asymptotic standard error of the estimator is given by ${\Sigma }^{-1/2}(\widehat{\theta }-{\theta }^{*})\stackrel{d}{\to }\mathcal{N}(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;\widehat{\theta },({\widehat{f}}^{(\setminus m)},{\widehat{g}}^{(\setminus m)}))\varphi {(D_i;\widehat{\theta },({\widehat{f}}^{(\setminus m)},{\widehat{g}}^{(\setminus m)}))}^{\prime }\right){\widehat{J}}^{-1\prime },\hfill \end{array}$$

  (38)

  $$\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{g}}^{(\setminus m)})){|}_{\theta =\widehat{\theta }}.\hfill \end{array}$$

  (39)

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

2.5. 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 in terms of BME. 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). Suppose that we attempt to explicitly model W with a latent variable $Z\in \{0,1\}$ and define the distribution of Z and its influence on the outcome as $Q_{\beta }\left(Z\right|X_{[0,min(T,C\left)\right]},A_{[0,min(T,C\left)\right]})$ and $Q_{\kappa }\left(T\right|Z,X_{[0,min(T,C\left)\right]},A_{[0,min(T,C\left)\right]})$ with parameters $\beta$ and $\kappa$, respectively. For example, in the simulation study described below, we use a logistic model that determines the distribution of the latent variable based on baseline covariate values:

$$\begin{array}{c}\hfill Q_{\beta }\left(Z\right|X_{[0,min(T,C\left)\right]},A_{[0,min(T,C\left)\right]})={\displaystyle \frac{1}{1+exp(-(2Z-1)({\beta }^{\prime }{X}_0+{\beta }_0))}},\end{array}$$

(40)

and we extend the hazard model as

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

(41)

The same argument for Proposition 1 applies to this model. Suppose that the ground truth is described by

$$\begin{array}{c}\hfill P\left(W\right|X_{[0,min(T,C\left)\right]},A_{[0,min(T,C\left)\right]})={\displaystyle \frac{1}{1+exp(-(2W-1)({\beta }^{*T}{X}_0+{\beta }_0^{*}))}},\end{array}$$

(42)

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}$$

(43)

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 assuming that the model is identifiable, one obtains the same result as in Proposition 1.

Then, one can also derive Neyman-orthogonal scores for this model as well. The logarithm of marginal likelihood of outcome is given by

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

(44)

with the aid of variational distributions ${\left\{r_i\right\}}_{i\in \mathcal{I}}$. Hereafter, we use the variable without time subscripts, such as $A_i,X_i$, which denotes the collection of variables throughout the observation period and can be understood from the context. Although optimal values for ${\left\{r_i\right\}}_i$ depend on $\theta ,\kappa ,\beta$ and f, variable $r_i$ can also be regarded as an independent nuisance parameter. Naive non-orthogonal scores for ML estimation are obtained by simply differentiating the marginal likelihood with respect to $\theta$ as

$$\begin{array}{ccc}& & {\varphi }_{\mathrm{naive},k}(D;\theta ,f,\kappa ,r)=\hfill \\ & & \sum _{Z\in \{0,1\}}r\left(Z\right)\left\{I_{[0,C]}\left(T\right)A_{k,T}exp(-{\theta }_k)-{\int }_0^{min(T,C)}{A}_{k,t}exp(f\left(X_t\right)+\kappa Z)dt\right\}.\hfill \end{array}$$

(45)

These score functions can be orthogonalized in the same manner as for the model without latent variables.

Proposition 3 (Neyman orthogonal scores with a binary latent variable) 

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

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

(46)

($E\subset \mathcal{X}\times \{0,1\}$) are absolutely continuous to each other with Radon-Nykodym derivatives

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

(47)

for which ${\mu }_k$ denotes ${\mu }_a$ such that $a_k=1$, and $g_k^{*}$ is a measurable function on $\mathcal{X}\times \{0,1\}$. Here, we have used random variable $r^{*}\left(Z\right)$

$$\begin{array}{c}\hfill r^{*}\left(Z\right)={\displaystyle \frac{Q_{\tilde{h}}\left(T\right|Z,X,A)Q_{\beta }\left(Z\right|X,A)}{{\sum }_{\tilde{Z}\in \{0,1\}}{Q}_{\tilde{h}}\left(T\right|\tilde{Z},X,A)Q_{\beta }\left(\tilde{Z}\right|X,A)}}.\end{array}$$

(48)

Then, score functions for estimating ${\left\{{\theta }_k\right\}}_{k\in \mathcal{K}}$

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

(49)

are Neyman orthogonal scores that can be used to debias the maximizer of the marginal log-likelihood in Eq.(44).

Proof. 

The outline of the proof is the same as that for Proposition 2, and so is omitted here. □

2.6. Inference with Multiple-Kernel Models

In the simulation study described below, we use multiple-kernel models for the estimation of f and $g_k$[26]. Statistical properties of estimation with multiple-kernel models under a 1-norm, 2-norm or mixed-norm regularization have been extensively studied [27,28,29,30]. We employ the 2-norm-regularized version for which calculation of BME is tractable. Regularised ML estimation of $(\theta ,f)$ (or g) with discretized timesteps under a 2-norm regularization is formulated as

$$\begin{array}{c}\hfill \underset{\theta ,f}{max}\underset{0\le t=m\Delta t<min(C_i,T_i)+\Delta t}{{\sum }_{i\in \mathcal{I}}}\mathcal{L}(D_{i,t};\theta ,f)-\sum _k{\displaystyle \frac{{\lambda }_k}{2}}{\alpha }_k^{\prime }{G}_k{\alpha }_k,\end{array}$$

(50)

where the Gram matrix, ${\left(G_k\right)}_{it,i^{\prime }{t}^{\prime }}=k(X_{i,t},X_{i^{\prime },t^{\prime }})$ relates f and ${\alpha }_k$ via

$$\begin{array}{c}\hfill f\left(X_{i,t}\right)=\sum _k\sum _{i^{\prime }{t}^{\prime }}{\left(G_k\right)}_{it,i^{\prime }{t}^{\prime }}{\alpha }_{k,i^{\prime }{t}^{\prime }}\end{array}$$

(51)

, and

$$\begin{array}{c}\hfill \mathcal{L}(D_{i,t};\theta ,f)=I_{[t,t+\Delta t)}\left(T_i\right)({\theta }^{\prime }{A}_{it}+f\left(X_{it}\right))-exp({\theta }^{\prime }{A}_{it}+f\left(X_{it}\right))\Delta t\end{array}$$

(52)

denotes the timestep-wise likelihood function. Concretely, we use linear or 1–3 dimensional Gaussian kernels, such as $k_j(X_{i,t},X_{i^{\prime },t^{\prime }})=X_{ij,t}^{\prime }{X}_{i^{\prime }j,t^{\prime }}$ and $k_j(X_{i,t},X_{i^{\prime },t^{\prime }})=exp(-{\sum }_{1\le \ell \le d}{(X_{i{j}_{\ell },t}^{\prime }-X_{i^{\prime }{j}_{\ell },t^{\prime }})}^2/2{\sigma }_d^2)$. See the Results section for concrete designs of kernels. Since ${\left\{{\alpha }_k\right\}}_k$ has large dimensions, we use incomplete Cholesky decomposition [31] to approximate $N\times N$ Gram matrix with $N\times N^{\prime }$ matrix $L_k$ as $G_k\approx L_k{L}_k^{\prime }$ for $N^{\prime }\ll N$. We then have $f={\sum }_k{L}_k{u}_k$ with $u_k=L_k^{\prime }{\alpha }_k$ and ${\displaystyle \frac{{\lambda }_k}{2}}{\alpha }_k^{\prime }{G}_k{\alpha }_k={\displaystyle \frac{{\lambda }_k}{2}}{\parallel u_k\parallel }_2^2$. We solved this maximization problem, by applying limited-memory BFGS algorithm [32] with backtracking and line search stopped by the Armijo condition. The optimization was stopped, if the ${\ell }_2$-norm of the gradient gets smaller than ${ϵ}_{\mathrm{stop}}=1.0\times {10}^{-2}$. For ML estimation, one can approximate the BME by using the Laplace approximation. Assigning the Hessian of the log-likelihood with the regularization term to $\widehat{H}$ and regarding the regularization term as the log-likelihood of a Gaussian process prior [33], we calculate BME (see, e.g., Chapter 3 of Ref.[23]) using Gaussian integrals as follows:

$$\begin{array}{c}\hfill ln{P}_{\mathrm{BME}}\left(\mathcal{M}\right|D)\approx \underset{0\le t=m\Delta t<min(C_i,T_i)+\Delta t}{{\sum }_{i\in \mathcal{I}}}\mathcal{L}\left(D_{i,t};\widehat{\theta },\sum _k{L}_k\widehat{u_k}\right)+\sum _k{\displaystyle \frac{{\lambda }_k}{2}}\dim \left(u_k\right)-{\displaystyle \frac{1}{2}}ln\left|\widehat{H}\right|.\end{array}$$

(53)

In the above approximation, the space of functions perpendicular to the range of $G_k$ does not contribute to $\mathcal{L}$ (the representer theorem [34]). Therefore, the contributions of the prior and posterior cancel out. Similarly, the space perpendicular to the range of $L_k$, but not to that of $G_k$, makes negligible contributions to $\mathcal{L}$ relative to their variations in the prior process.

The above procedure can also be used to obtain $g_k$, where

$$\begin{array}{c}\hfill {\mathcal{L}}_{g_k}(D_{i,t};g_k)=A_{ik,t}ln(1+e^{-g_k\left(X_{i,t}\right)})+(1-\sum _{\ell }{A}_{i\ell ,t})ln(1+e^{g_k\left(X_{i,t}\right)})\end{array}$$

(54)

replaces $\mathcal{L}$ for $\theta$ and f in the above.

In the simulation with latent-variable Z, we applied the EM algorithm for optimization (see, e.g., Chapter 9 of Ref.[23]). Since the latent variable model has indefiniteness with respect to permutation of the range of Z as do most other models with discrete latent variables, we started optimization from values close to ${\theta }^{*},f^{*},{\beta }^{*},{\kappa }^{*}$. Repeating sequences of updating $r(Z_i=1|X,A,T)$ with Eq.(48), optimization of $\theta ,f$ and a gradient-ascent step in $\beta ,\kappa$, we maximized the marginal log-likelihood in Eq.(44). Nuisance parameter $g_k$ was obtained by performing logistic regression with ${\mathcal{L}}_{g_k}(D_{i,t},Z_i;g_k)={\sum }_{Z_i\in \{0,1\}}\widehat{r_i}\left(Z_i\right)(A_{ik,t}ln(1+e^{-g_k(X_{i,t},Z_i)})+(1-{\sum }_{\ell }{A}_{i\ell ,t})ln(1+e^{g_k(X_{i,t},Z_i)}))$. In this case, we estimated $g_k(X,Z)=g_{k,1}\left(X\right)+Z{g}_{k,2}\left(X\right)$ with $g_{k,1},g_{k,2}\in {\mathcal{H}}_{k_X}$.

In the numerical analysis, all covariate values are normalized to have a zero mean and a unit variance. For ${\sigma }_d$ and ${\lambda }_k$, we performed a grid search, 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$). In this grid search, we use the same values of ${\lambda }_k$ and ${\sigma }_d$ for all d-dimensional Gaussian kernels and use ${\lambda }_k=0$ for linear kernels. The goodness of fit was measured in terms of the average of BME for ${\left\{D_{\mathrm{train}}^{\left(m\right)}\right\}}_{1\le m\le 10}$ and the error defined in Eq.(35).

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.

To simulate 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 ground-truth 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 used 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 and $g_k$. 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}$$

(55)

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}$, where ${\mathcal{H}}_k$ and b are a reproducing-kernel Hilbert space (RKHS) associated with kernel k and a bias parameter, respectively. 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,X^{\prime })=exp(-{(X-X^{\prime })}^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 know from the model construction, covariate $X_{2t}$ is necessary for describing 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 (Fig.2(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$ (Fig.2(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. For a correctly specified model, such violation of the assumption is not detected (Fig.2(C)). The violation of the time-homogeneity assumption is observed for models misspecified in different manners (which we omit to avoid redundancy).

The nuisance parameter $g_{\ell }$ can also be similarly modeled as

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

(56)

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. In practice, on top of the kernels described above, we additionally included two-dimensional Gaussian kernels such as $k_{12}(X,X^{\prime })=exp(-({(X_1-X_1^{\prime })}^2+{(X_2-X_2^{\prime })}^2)/{\sigma }_2^2),$ and examined whether the inclusion of these additional kernels improves BME. In the present case, the two-dimensional kernels did not improve BME. In the tuning of hyperparameter values, we observe that BME and ${\mathrm{CVErr}}_{g_{\ell }}$ identify different optimal values (Fig.2(D) and (E)). We show debiased estimates of ${\theta }_1$ based on $g_1$ estimated with each of these optimal hyperparameter values ${\lambda }_g=70.0$ and ${\lambda }_g=1.0$ in Fig.2(F). 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 Fig.2(F) showed that the debiased estimator actually identified the ground truth with minimal bias, in contrast with the naive ML estimator. In Fig.2(F), we also observe 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 }\mathcal{N}(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)({\beta }^{*\prime }{X}_0^{\left(\mathrm{test}\right)}+{\beta }_0^{*}))}},\end{array}$$

(57)

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 }^{*}=2.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 3(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 EM algorithm (see e.g., Chapter 9 of Ref.[23]), which led to further improvement of BME (Figure 3(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 ground truth with minimal bias, in contrast with the naive ML estimator (Figure 3(C)).

4. Discussion

Preceding studies [3,4,5] showed that the Cox model allows multiple, contradictory interpretations. In the present study, we therefore proposed to estimate hazard ratios of treatment options in an exponential hazard model with or without latent variables that lacks baseline hazard as a function of time elapsed after enrollment. We assumed that conditional hazard should not depend on a specific time during the observation, and measured the degree of its violation as an indicator of the necessity to seek additional covariates or latent variables. Given this and several other 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. Then, we proposed methods to debias ML estimators of hazard ratios in cases with and without latent variables, simultaneously validating the model assumptions by calculating BME. Numerical simulations confirm that the proposed method allows us to systematically identify necessary variables and estimate the treatment effect 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 (thanks to the cancellation in the representation of partial likelihood). 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. [22,38]). The present results suggest the need to reconsider this situation by reevaluating the tradeoff between the efficiency and uninterpretability due to the indefinite baseline hazard. The disadvantage of using MPL is not only the uninterpretability. First of all, the indefinite 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 [18], 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 the conventional method based on 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 in causal inference for observational studies with uncontrolled dynamic treatment and real-time measurement of a large number of covariates.

In technical aspects, for the sake of clear exposition, we restricted our study to a relatively simple setting, in which 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. For the first condition, the advantage of ML estimation has been demonstrated in the context of informative censoring [41]. Therefore, extending our framework in this direction seems promising. Second, the use of multiple independently changing treatment variables is straightforward, but requires cumbersome modifications, where finding the root of scores amounts to solving multivariate algebraic equations of degrees greater than one. To apply our framework to models with continuous treatment variables, the conditional distribution of treatment assignment for each value of covariates must be identified in a manner similar to what is done with propensity scores for continuous treatment [42]. There is still a technical challenge in this step, because machine learning of the conditional distribution of continuous variables remains in the phase of theoretical argument (see e.g., Refs.[43,44]), not practical use. 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 [45,46,47,48] and numerical examination of identifiability in the context of biological modeling [49], there is still a technical challenge. Related to this topic, using a singular model can lead to effective breakdown of asymptotic normality of estimators [50]. 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 [51] (or information criteria for singular parametric models [52,53]). 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 [54,55]) is yet an attractive, worthy challenge. In this challenge, our technique to use variational distributions as a nuisance parameter in the orthogonalization will be useful.