text:
\section{Introduction}
Composite fermions (CF's) [\onlinecite{jain}] describe the physics of electrons in fractional quantum Hall (FQH) regime. At filling factor $\nu=1/2$, essentially they absorb the external flux, and make a metallic state [\onlinecite{hlr}] with its own Fermi surface - Fermi surface of CF's. By slightly modifying Read's dipole construction of composite (neutral) fermions in the half-filled lowest Landau level [\onlinecite{read}], an argument can be given for the accumulation of Berry phase equal to $\pi$ as a CF encircles its own Fermi surface [\onlinecite{ws}].
This has motivated a description of the CF's in this setting in terms of Dirac fermions, which has been recently introduced in Ref. [\onlinecite{son}], and have attracted some interest [\onlinecite{ws,gz,psv,msone,kach,mrf,btj,bmf,mv}].
The PH symmetric description of the half-filled LL is given in terms of a Dirac system of composite quasiparticles - Dirac CF's at a finite chemical potential [\onlinecite{ws}] and in the presence of a gauge field. However, the implied existence of singularity at zero momentum in the CF
spectrum was criticized [\onlinecite{fdmh,bj}].
We may add that, due to the requirement of gauge invariance in two dimensions (2D), a small mass
must be introduced into the Dirac theory (``parity anomaly") [\onlinecite{red}]. This may be a way
to heal and complete in the high-energy domain (``UV completion" [\onlinecite{ws2}]) Dirac
description of CF's, and avoid singularity.

Thus the description in terms of Dirac fermions
may have capacity to capture essential, at least qualitative, aspects of the CF's physics. To further examine this possibility in this work we consider BCS pairing of Dirac CF's. First, in the framework of Dirac description of a single CF, we point out that, assuming Cooper pairing between spinor components, besides so-called PH symmetric Pfaffian, also anisotropic states  can be realized. This is  analogous to the ${^3}$He system in which both  $B$ and $A$ (anisotropic) phases are possible [\onlinecite{vol}].  Next we discuss unconventional $p$-wave  pairing of two kinds of Dirac CF's, motivated by the situation in the quantum Hall bilayer (QHB) at total filling factor one, i.e.
 with each layer  with half-filled lowest LL. In this system $p$-wave pairing between two kinds of non-relativistic Halperin-Lee-Read (HLR) composite fermions at intermediate interlayer distances was proposed in Ref. [\onlinecite{msr}], and, recently, this scenario was further substantiated by a detection of the topological signatures of the $p$-wave system in the torus geometry [\onlinecite{mdp}]. Therefore, it is natural to ask how this picture may be modified if we take into account the description by two Dirac CF's of the two half-filled LL  monolayers, and consider their possible pairing.

One of the main conclusions that we can draw by applying the Dirac CF formalism in the context of BCS pairing is that
due to the Dirac two-component nature, isotropic (gapped) pairing states may coexist with anisotropic ones, and this is in the accordance with the results on PH symmetric, single and bilayer fractional quantum Hall systems obtained by employing  exact diagonalization [\onlinecite{park,wy}], as well as  with experimental findings [\onlinecite{l1,l2}], in which anisotropy is probed by in-plane magnetic field. This may be a direct consequence of the dipole nature of CF's that is captured by Dirac formalism. Anisotropic pairing states may serve as seed states for Pfaffian and anti-Pfaffian through a process in which  PH asymmetry increases by introducing a mass term, while rotational symmetry gradually sets in. Furthermore, we find that the features, in particular low-energy spectrum, of the QHB at intermediate distances between the layers are better captured if we assume Dirac rather than HLR $p$-wave paired CF's at large distances (decoupled layers).
Already at the  effective field theory level, modeling the evolution with the distance between layers by Dirac CFs, we can detect the main feature of CF-composite boson (CB) mixed states [\onlinecite{msr,srm}]:  the decrease in the number of CF's with decreasing distance.

The paper is organized as follows. In Section II , based on Dirac formalism, we discuss
the single layer case and its pairing instabilities,
including the situation when the PH symmetry is spoiled by a mass term. In Section III we discuss the pairing instabilities in the bilayer system when the PH symmetry inside each layer is intact. In Section IV we examine the evolution of low-energy properties of the QHB with distance between layers, by including a mass term with an opposite sign in the two layers. The last section, Section V, is devoted to discussion and conclusions. Mean-field analysis of the coexistence of the isotropic and anisotropic pairings is presented in the Appendices.
\section{Dirac composite fermion and Cooper pairing}
We begin by considering a single Dirac fermion which was  proposed to effectively describe
half-filled lowest Landau level of electrons [\onlinecite{son}], with $s$-wave pairing between spinor components. The $s$-wave pairing suggested in Ref. [\onlinecite{son}], can be expressed by the following Bogoliubov - de Gennes Hamiltonian in the Nambu - Gorkov notation,
\begin{eqnarray}\label{BdGsD}
H &=&\frac{1}{2} \sum_{\bf k} \left[
  \begin{array}{cc}
   \Psi^\dagger ({\bf k}) & \tilde{\Psi}(-{\bf k}) \\
  \end{array}
\right]  \\ \nonumber
&\times&\left[
  \begin{array}{cc}
   {\cal D} ({\bf k}) & {\cal P}({\bf k}) \\
   {\cal P}^\dagger ({\bf k}) & - {\cal D} (-{\bf k}) \\
  \end{array}
\right]
\left[
  \begin{array}{c}
   \Psi ({\bf k}) \\
   \tilde{\Psi}^\dagger (-{\bf k}) \\
  \end{array}
\right],
\end{eqnarray}
where $\Psi ({\bf k})$ denotes a two-component spinor with momentum ${\bf k}$,
\begin{equation}
\Psi ({\bf k}) = \left[
  \begin{array}{c}
   \Psi_{a} ({\bf k}) \\
   \Psi_{b} ({\bf k}) \\
  \end{array}
\right], \,\,\,\,
\tilde{\Psi} ({\bf k}) = \left[
  \begin{array}{c}
   \Psi_{b} ({\bf k}) \\
   \Psi_{a} ({\bf k}) \\
  \end{array}
\right],
\end{equation}
and
\begin{equation}
{\cal D}({\bf k}) = \left[
  \begin{array}{cc}
   - \mu & k_x - i k_y \\
   k_x + i k_y & - \mu \\
  \end{array}
\right]=-\mu\sigma_0+k_x \sigma_x+k_y \sigma_y, \label{dmatrix}
\end{equation}
and  $2 \times 2$ matrix ${\cal P}({\bf k})$ describes Cooper pairing between $a$ and $b$ spinor components
\begin{equation}
{\cal P}({\bf k}) = \left[
  \begin{array}{cc}
   \Delta_s & 0\\
    0 & - \Delta_s \\
  \end{array}
\right]=\Delta_s \sigma_z  , \label{right_pairing_sD}
\end{equation}
or more explicitly
\begin{equation}
\delta {\cal H} = \sum_{\bf k} \{ - \Delta_s \Psi_{a}({\bf k}) \Psi_{b}(-{\bf k}) +  h.c. \}. \label{rightH}
\end{equation}
Here, $\sigma_0$ is the $2\times2$ identity matrix, while ${\bm \sigma}$ are the standard Pauli matrices. Throughout the paper we set $\hbar = 1$, and the Fermi velocity, $v_F = 1$.  $\mu$ denotes a chemical potential equal to $\mu = \sqrt{B} = k_F$, where $B$ and $k_F$ are the external magnetic field and Fermi
vector, respectively.

Since the pairing matrix anticommutes with the free Dirac Hamiltonian at the zero chemical potential, the dispersion of  Bogoliubons has
 the rotationally symmetric form
\begin{equation}
E_k^2 = (k \pm \mu)^2 + \Delta_{s}^2,
\end{equation}
where $k\equiv |{\bf k}|$.
This construction is considered in the literature
as a basis for a PH symmetric Pfaffian system.

However,  a different type of pairing is also possible with the pairing matrix
\begin{equation}
{\cal P}({\bf k}) = \left[
  \begin{array}{cc}
   0  & \alpha k_x + \beta k_y\\
    \alpha k_x - \beta k_y & 0 \\
  \end{array}
\right],
\end{equation}
or more explicitly
\begin{eqnarray}
\delta {\cal H}'& =&\sum_{{\bf k}} \alpha k_x  \{ \Psi_{a}^\dagger ({\bf k}) \Psi_{a}^\dagger(-{\bf k}) + \Psi_{b}^\dagger({\bf k}) \Psi_{b}^\dagger(-{\bf k}) \} + h.c. \nonumber\\
 &+&\sum_{\bf k}  \beta k_y  \{ \Psi_{a}^\dagger({\bf k}) \Psi_{a}^\dagger(-{\bf k}) - \Psi_{b}^\dagger({\bf k}) \Psi_{b}^\dagger(-{\bf k}) \} + h.c.  \nonumber \\
\label{anisotropicH}
\end{eqnarray}
where $\alpha$ and $\beta$ are, in general, allowed to be complex coefficients. The overall form of $ \delta {\cal H}'$ is fixed by the requirement of the $CP$ symmetry, which, as emphasized in Ref. [\onlinecite{son}], is equivalent to the requirement of the PH symmetry in the real electron system. Namely, the $CP$ symmetry is a product of the charge conjugation, $C$,
\begin{equation}
C \Psi({\bf k}) C^{-1} = \sigma_x \Psi^* ({\bf k}),
\end{equation}
and a parity transformation, $P$,
\begin{equation}
P \Psi({\bf k}) P^{-1} =  \Psi^* ({\bf k}'),
\end{equation}
where ${\bf k} = (k_x, k_y) \rightarrow {\bf k}' = (k_x, - k_y)$ under  the parity transformation. Thus,
\begin{equation}
CP \Psi({\bf k}) (CP)^{-1} = \sigma_x \Psi({\bf k}').
\label{CP}
\end{equation}
The starting Dirac Hamiltonian\ (\ref{BdGsD}) with ${\cal P}=0$, as well  $ \delta {\cal H}'$ are both invariant under the $CP$ transformation  (\ref{CP}).  Notice that (\ref{rightH}) is invariant up to a sign change under $CP$ transformation. This is also a property of the small mass term that seems necessary to ensure the gauge invariance of the theory, and to avoid the singularity at ${\bf k} = 0$ [\onlinecite{ws2}]. The BCS pairing terms as the one in (\ref{rightH}) may accommodate the sign change by gauge transformations [\onlinecite{son}]. Thus the theory is invariant under $CP$ transformation in a more general sense, allowing for terms that are invariant up to a change of the sign. This
makes our choice for $p$ wave not unique. Indeed, other $p$ wave pairing order parameters are also possible, including one analogous to the A phase of ${^3}$He system that features two (gapless) Fermi points. This case can be analyzed analogously to the one considered here, and the main conclusions remain.  In the following, we restrict our discussion to the
$p$-wave case (\ref{anisotropicH}) invariant under $CP$ transformation in the strict sense.

We now consider the pairings given by Eq.\ (\ref{anisotropicH}), recently also discussed in Ref. [\onlinecite{wc}], in light of the possibility of introducing an anisotropy. The choice $\alpha = \Delta$ and  $\beta  = - i \Delta $ yields the pairing matrix, $ {\cal P}({\bf k})$, proportional to the Dirac Hamiltonian, ${\cal D}({\bf k})$, at chemical potential $\mu = 0$, and thus explicitly rotationally invariant. (See
 also Sec. III for further analysis of the rotational symmetry.). In that case, the dispersion relation of Bogoliubons, $E_{\bf k}^2 = k^2(1+\Delta^2)+\mu^2 \pm 2k \sqrt{\mu^2 + k^2 \Delta^2}$, implies that the the pairing just renormalizes
 chemical potential. On the other hand, by
 choosing $\alpha = \Delta$ and  $\beta  = + i \Delta $, we obtain
 \begin{equation}
E_{\bf k}^2 = k^2 (1 + \Delta^2) + \mu^2 \pm 2 \sqrt{\mu^2 k^2 + \Delta^2 (k_x^2 - k_y^2)^2}. \label{andispersion}
\end{equation}
This dispersion describes an {\em anisotropic} gapless system with four nodes at
\begin{equation}
k_x = \pm \frac{\mu}{\sqrt{1 - \Delta^2}}, \;\;\;\; {\rm and}\;\;\; k_y = 0,
\end{equation}
and
\begin{equation}
k_y = \pm \frac{\mu}{\sqrt{1 - \Delta^2}}, \;\;\;\; {\rm and}\;\;\; k_x = 0.
\end{equation}
The appearance of the four nodes related by the discrete $C_4$ symmetry is a consequence of the  $C_4$ symmetry of the pairing\ (\ref{anisotropicH}) with $\alpha = \Delta$ and  $\beta  = + i \Delta $.
In fact, Eq. (\ref{anisotropicH}) describes a whole family of gapless anisotropic solutions.

 If we consider both $s$-wave (\ref{rightH}) and $p$-wave (\ref{anisotropicH}) with $\alpha = \Delta$ and  $\beta  = + i \Delta $ pairings, the dispersion of the Bogoliubov quasiparticles is
\begin{equation}
\tilde{E}_{\bf k}^2 = \Delta_s^2 + k^2 (1 + \Delta^2) + \mu^2 \pm 2 \sqrt{\mu^2 k^2 + \Delta^2 (k_x^2 - k_y^2)^2},
\end{equation}
i.e. the dispersion (\ref{andispersion}) simply acquired a shift of
$\Delta_s^2$ in the presence of the $s$-wave pairing. This is a consequence of the anticommutation of the matrices corresponding to the two pairings, similarly to the situation in
Refs. [\onlinecite{RJ2014}, \onlinecite{roy-goswami}], which makes their coexistence likely at a finite chemical potential. Assuming a generic form of the two couplings driving the instabilities in the isotropic and anisotropic channels, in the presence of a small mass term, we show in Appendix A that the low energy description implies that the isotropic instability (\ref{rightH}) may coexist with the anisotropic one. This is consistent with experimental [\onlinecite{l1,l2}], and theoretical  [\onlinecite{park,wy}] findings pointing out that gapped states at half-filled Landau level can sustain and even harbor anisotropy.

In connection with the possible pairings given by Eq.\ (\ref{anisotropicH}) when $\alpha = \Delta$ and  $\beta  = + i \; \Delta $, we may notice that if we break $CP$ (particle-hole symmetry) by a mass term of the form $\sim \Psi^\dagger({\bf k}) \sigma_3 \Psi({\bf k})$,
one component, $a$ or $b$, of the Dirac field will remain in the low energy  sector. The remaining fermion should correspond to HLR (spinless) fermion which in turn pairs in the manner of $p$-wave. This should correspond to Pfaffian and anti-Pfaffian states (that comprise possible $( k_x \pm i k_y)$ states), in the absence of PH symmetry, but with an emergent rotational symmetry. A closely related proposal for the existence of the Pfaffian (Moore-Read) state in the presence of an excitonic instability already appeared in the context of Dirac CF physics in graphene
Ref. [\onlinecite{cai}].

 To further understand the pairings in Eqs.\ (\ref{rightH}) and (\ref{anisotropicH}), we now consider the chirality operator $\frac{{\bm\sigma} \cdot{\bf k}}{|k|}$, and its eigenstates
\begin{equation}
|+ \rangle = \frac{1}{\sqrt{2}} \left[
  \begin{array}{c}
   1 \\
   \frac{k_+}{k} \\
  \end{array}
\right], \,\,\,\,
|- \rangle = \frac{1}{\sqrt{2}} \left[
  \begin{array}{c}
   - 1 \\
   \frac{k_+}{k} \\
  \end{array}
\right].\label{chiral-eigenstates}
\end{equation}
We can  introduce Dirac operators with a definite chirality
\begin{equation}
\Psi_+ ({\bf k}) = \frac{1}{\sqrt{2}}( \Psi_a ({\bf k}) +  \frac{k_-}{k} \Psi_b ({\bf k})) ,
\end{equation}
and
\begin{equation}
\Psi_- ({\bf k}) = \frac{1}{\sqrt{2}}( - \Psi_a ({\bf k}) +  \frac{k_-}{k} \Psi_b ({\bf k})) ,
\end{equation}
to find that
\begin{equation}
\Psi_a ({\bf k})  \Psi_b ( - {\bf k})
=- \frac{1}{2}  \frac{k_+}{k} [\Psi_+ ({\bf k})  \Psi_+ ( - {\bf k}) +   \Psi_- ({\bf k}) \Psi_- (-{\bf k})],
\label{s_into_p}
\end{equation}
with $k_\pm\equiv k_x\pm i k_y$.
We can clearly see from Eq.\ (\ref{s_into_p}) that in the chirality basis i.e. the eigenbasis of the non-interacting system, the pairing (\ref{rightH}), in fact, describes a pairing in the odd ($p$-wave) channel. This
can be understood as a consequence of the non-trivial Berry phase contributions, as discussed in Ref. [\onlinecite{son}]; see also Ref. [\onlinecite{mn}]
for the influence of the singularities (topological charges) on the vorticity of Cooper pairs. On the other hand,
the anisotropic pairing (\ref{anisotropicH}) is a combination of odd channel components in the chirality basis.

We now analyze an alternative scenario for the coexistence with the $p$-wave pairing represented by the pairing matrix ${\cal P}({\bf k})=(\alpha k_x+\beta k_y)\sigma_x$ that features two Fermi points and does not require a mass for the coexistence with the isotropic state. In particular, as shown in Appendix B,  a special anisotropic pairing with
\begin{equation}\label{anisotropicHH}
{\cal P} \sim i  k_y\sigma_x
\end{equation}
can coexist with the isotropic pairing. Analogously, we can discuss pairing with ${\cal P}({\bf k})\sim(\gamma k_x+\delta k_y)\sigma_y$,
where $\gamma$ and $\delta$ are, in general, allowed to be complex coefficients.
The ensuing pairing is then given by
\begin{equation}\label{anisotropicHH2}
{\cal P}({\bf k})\sim k_x \sigma_y.
\end{equation}
Both these pairings
are invariant up to a change of sign (up to a gauge transformation) under the $CP$ transformation. Each pairing on its own features two Fermi points, and is likely energetically advantageous over the pairing in (\ref{anisotropicH}) that has four Fermi points. As we explicitly show in Appendix B, these pairings do not need a mass term to coexist with the isotropic state. Furthermore, in the presence of a mass term, they develop new components, and may thus evolve into the rotationally symmetric pairings of HLR fermions. These are the reasons that make states given by Eqs.\ (\ref{anisotropicHH}) or (\ref{anisotropicHH2}) likely present when considering pairing instabilities  in the half-filled LL, consistent with the exact diagonalization results  of Refs. [\onlinecite{park,wy}].

 Finally, we  point out that the Dirac based microscopic wave functions of pairing instabilities have not been proposed and
tested yet. The effective field theory approach seems currently to be the most efficient tool for treating the Dirac-based pairing instabilities and their properties. Once the microscopic description is provided, most importantly for the case of PH Pfaffian, anisotropic modifications may be induced in the manner described and discussed in Refs. [\onlinecite{h}, \onlinecite{bj2}].

\section{Dirac fermions and $p$-wave  pairing}

We consider the following general form of the Bogoliubov-de Gennes Hamiltonian,  motivated by the situation in a QHB system with each of the two layers at half filling,
\begin{eqnarray}\label{BdG}
H &=& \sum_{\bf k} \left[
  \begin{array}{cc}
   \Psi_\uparrow^\dagger ({\bf k}) & \Psi_\downarrow (-{\bf k}) \\
  \end{array}
\right]  \\ \nonumber
&\times&\left[
  \begin{array}{cc}
   {\cal D}_\uparrow ({\bf k}) & {\cal P}({\bf k}) \\
   {\cal P}^\dagger ({\bf k}) & - {\cal D}_\downarrow (-{\bf k}) \\
  \end{array}
\right]
\left[
  \begin{array}{c}
   \Psi_\uparrow ({\bf k}) \\
   \Psi_\downarrow^\dagger (-{\bf k}) \\
  \end{array}
\right],
\end{eqnarray}
where $\Psi_\uparrow ({\bf k})$ and $\Psi_\downarrow ({\bf k})$ are two component spinors,
\begin{equation}
\Psi_\uparrow ({\bf k}) = \left[
  \begin{array}{c}
   \Psi_{a \uparrow} ({\bf k}) \\
   \Psi_{b \uparrow} ({\bf k}) \\
  \end{array}
\right], \,\,\,\,
\Psi_\downarrow ({\bf k}) = \left[
  \begin{array}{c}
   \Psi_{b \downarrow} ({\bf k}) \\
   \Psi_{a \downarrow} ({\bf k}) \\
  \end{array}
\right].
\end{equation}
Matrices ${\cal D}_\uparrow ({\bf k})$ and ${\cal D}_\downarrow ({\bf k})$ describe two identical Dirac systems, ${\cal D}_\uparrow ({\bf k})={\cal D}_\downarrow ({\bf k})={\cal D} ({\bf k})$, with ${\cal D} ({\bf k})$ given by Eq.\ (\ref{dmatrix}), while $2 \times 2$ matrix ${\cal P}({\bf k})$ describes Cooper pairing between the two systems $\uparrow$ and $\downarrow$.

A triplet $p$-wave pairing between the same spinor components can be expressed
as the following term in the Hamiltonian,
\begin{eqnarray}
\delta {\cal H} &=& \sum_{\bf k} \{ [ \Delta^*_{\bf k} \Psi_{a \downarrow}(-{\bf k}) \Psi_{a \uparrow}({\bf k})  \\ \nonumber
&+&\Delta^*_{\bf k} \Psi_{b \downarrow}(-{\bf k}) \Psi_{b \uparrow}({\bf k}) ] \, + \, h.c. \},
\end{eqnarray}
with a pairing function $ \Delta_{\bf k} = \Delta (k_x \pm i k_y)$.
The corresponding pairing matrix in the Hamiltonian (\ref{BdG}) is
\begin{equation}
{\cal P}({\bf k}) = \left[
  \begin{array}{cc}
   0 & \Delta_{\bf k} \\
   \Delta_{\bf k} & 0 \\
  \end{array}
\right] = \Delta_{\bf k} \sigma_x.\label{anisotropic_pairing}
\end{equation}
A rotation around $z$ axis by an angle $\phi$ in both subsystems $\uparrow$ and $\downarrow$ is represented  by a matrix $ R = \exp (i \sigma_z  \phi/2)$ so that
\begin{eqnarray}
R \sigma_x R^{-1} &=& \sigma_x \cos \phi - \sigma_y \sin \phi, \\ \nonumber
R \sigma_y R^{-1} &=& \sigma_x \sin \phi + \sigma_y \cos \phi.
\end{eqnarray}
It can be readily seen that $ \tilde{R} H({\bf k}) \tilde{R}^{-1} \neq H({\bf k}')$ where
$ k_{x}^{'} = k_x\cos\phi  - k_y\sin\phi $ and $ k_{y}^{'} = k_x\sin\phi  + k_y\cos\phi $, and ${\tilde R}=\tau_0\otimes R$, with $\tau_0$ as the $2\times2$ unity matrix in the subsystem space.  Therefore, the system with the pairing matrix ${\cal P}({\bf k}) = \Delta_{\bf k}  \sigma_x$ is not rotationally invariant, and may lead to anisotropic behavior. In fact, the system is gapless, and supports two anisotropic Dirac cones at $ k_x^2 = \mu^2/(1 - \Delta^2) = k_0^2$ and $k_y = 0$. Expanding around $\pm k_0$ we obtain for $\Delta \ll 1$, $ E^2 \approx (1 - 2 \Delta^2) (\delta k_x)^2 + \Delta^2 (\delta k_y)^2$.
We find similar results if we choose,
\begin{equation}
{\cal P}({\bf k}) = \left[
  \begin{array}{cc}
   0 & \Delta_{\bf k} \\
   -\Delta_{\bf k} & 0 \\
  \end{array}
\right].\label{sign_anisotropic_pairing}
\end{equation}
Therefore the systems that we considered by now do not possess quantum spin Hall effect, and due to anisotropy are likely to be fragile under disorder, and certainly can not represent stable phases in realistic circumstances.

On the other hand, the system with the pairing matrix
\begin{equation}
{\cal P}({\bf k}) = \left[
  \begin{array}{cc}
   \Delta_{\bf k} & 0\\
    0 & - \Delta_{\bf k} \\
  \end{array}
\right] = \Delta_{\bf k} \sigma_z , \label{right_pairing}
\end{equation}
yields the dispersion relation of the Bogoliubons
\begin{equation}
E_{\pm} = \sqrt{(k \pm \mu)^2 + |\Delta_{\bf k}|^2}, \label{disper_triplet}
\end{equation}
and therefore resembles very closely $p$-wave pairing of ordinary fermions.
We now express the pairing in the chirality basis to obtain
\begin{eqnarray}
 &&\Delta^*_{\bf k} (\Psi_{b \downarrow}(-{\bf k}) \Psi_{a \uparrow}({\bf k}) -
 \Psi_{a \downarrow}(-{\bf k}) \Psi_{b \uparrow}({\bf k}))=- \Delta^*_{\bf k} \frac{1}{2}  \frac{k_+}{k} \nonumber \\
&\times&
 (\Psi_{+ \downarrow}(-{\bf k}) \Psi_{+ \uparrow}({\bf k}) -
 \Psi_{- \downarrow}(-{\bf k}) \Psi_{- \uparrow}({\bf k})). \nonumber \\
 \end{eqnarray}
 Thus depending whether  $ \Delta_{\bf k} = \Delta (k_x + i k_y)$ or
 $ \Delta_{\bf k} = \Delta (k_x - i k_y)$, we obtain $s$-wave or $d$-wave pairing,
 respectively, in the chirality basis.
 In this sense there is no surprise to find that the pairing matrix (\ref{right_pairing}) gives rise to a singlet state for $\uparrow$ and $\downarrow$ electrons. The choice for the pairing without the minus sign in Eq.\ (\ref{right_pairing}), i.e., ${\cal P}({\bf k})=\Delta_{\bf k}\sigma_0$, is not energetically favorable, since  pairing just renormalizes chemical potential in that case.

We now provide a topological characterization of pairing in Eq.\ (\ref{right_pairing}) through the  (pseudo)spin Chern number, $C_s$. In fact we find that in this case is $C_s = 1$, if we use the low-energy theory
with (\ref{right_pairing}) and $ \Delta_{\bf k} = \Delta (k_x + i k_y)$.  We calculated the Chern number by taking the eigenvectors
of the two lower Bogoliubov bands, $ |v_{-}({\bf k})\rangle$ and $ |v_{+}({\bf k})\rangle$,
corresponding to eigenvalues $ - E_{-}({\bf k})$ and $ - E_{+}({\bf k})$, respectively.
We first computed the Berry curvature of each vector,
\begin{equation}
F^{\sigma}_{xy}({\bf k}) = i  (\partial_x \langle v_{\sigma}({\bf k})| \partial_y |v_\sigma ({\bf k}) \rangle -  \partial_y \langle v_{\sigma}({\bf k})| \partial_x |v_\sigma ({\bf k}) \rangle ),
\end{equation}
and then the Chern number,
\begin{equation}
C_s = \frac{1}{2 \pi} \sum_\sigma \int d{\bf k} F^{\sigma}_{xy}({\bf k}),
\label{chern}
\end{equation}
where the sum in (\ref{chern}) is over the two lowest bands. Nevertheless, as discussed in the previous paragraph, and, also, due to the form of eigenvectors below, we expect that the real winding number is zero or two if a complete description is taken into account.

To further characterize  the pairing state, let us consider the four-component vectors of  the Bogoliubov bands with positive energy, $  E_{-}({\bf k})$ and $ E_{+}({\bf k})$,
\begin{eqnarray}
u_{-}(k) &=& \frac{1}{2 \sqrt{E_{-}}} \left\{ - \sqrt{E_{-} - (\mu - k)} \left(1, \frac{k_{+}}{k}\right),\right. \nonumber\\
 & & \left.\frac{\Delta \cdot k}{\sqrt{E_{-} - (\mu - k)}}  \left(-  \frac{k_{-}}{k},1\right) \right\},
\end{eqnarray}
and
\begin{eqnarray}
u_{+}(k) &=& \frac{1}{2 \sqrt{E_{+}}} \left\{  \sqrt{E_{+} - (\mu + k)} \left(1, - \frac{k_{+}}{k}\right),\right.\nonumber\\
 & &\left. \frac{\Delta \cdot k}{\sqrt{E_{+} - (\mu + k)}} \left(\frac{k_{-}}{k}, 1\right) \right\},
\end{eqnarray}
where we regrouped components to appear with common factors. In each  Bogoliubov eigenstate, the first two-component spinor, ( , ), is an eigenstate of the chirality
operator, given by Eq.\ (\ref{chiral-eigenstates}), while the second one is the eigenstate that is complex conjugated and with inverted components due to the ordering in the Nambu - Gorkov representation, and we fix $ \Delta_{\bf k} = \Delta (k_x + i k_y)$. From the coefficients in front of the fixed chirality states,  we find  the long-distance behavior of the pairing function ($g_{\bf k} \sim v_{\bf k}/u_{\bf k}$ in the usual BCS problem) in each band
\begin{equation}
g(z) \sim 1/|z| ,
\end{equation}
where $g(z)$ is the  pairing function in the real space, and $z = x + i y$. Thus the pairing function has the characteristic $s$-wave feature.

In this case the lowest gap
is at Fermi surface, $\Delta E \sim  \Delta \cdot k_F$, in  contrast with ordinary $p$-wave
pairing where the lowest gap is  at zero momentum, and it is equal to $\Delta E \sim k_F$ [\onlinecite{rg}].

\section{Quantum Hall bilayer and $p$-wave paired composite fermions}

 In light of recent advance in understanding  of each (isolated PH symmetric) half-filled  quantum Hall monolayer
based on Dirac CF's
it is quite natural to consider the physics of the bilayer, especially at the intermediate distances, in the same framework. It is important to take into account the $p$-wave
pairing [\onlinecite{msr}] which was initially expressed in terms of ordinary  HLR CF's. The picture based on the ordinary CF's does not have a clear answer for the lowest lying spectrum which appears nearly gapless (with small gap) or gapless when the system is put on a  torus, while the topological $p$-wave pairing of ordinary  fermions [\onlinecite{rg}] would likely produce a clear gap of the order $\mu$. However, even if we neglect possible insufficiencies with $p$-wave pairing of ordinary fermions, it is fundamentally important to address the problem of the QHB  in terms of Dirac CF's.

First we may  notice that the presence of the interlayer Coulomb interaction, which increases with decreasing distance between layers, spoils PH symmetry inside a layer.
We incorporate this breaking of the PH symmetry by introducing a mass $r$ in the Dirac matrices ${\cal D}_\uparrow ({\bf k})$ and ${\cal D}_\downarrow ({\bf k})$,  with opposite signs in each layer,
\begin{equation}
{\cal D}_\uparrow ({\bf k}) = \sigma_x k_x + \sigma_y k_y - \mu + r \sigma_z  = {\cal D}_\downarrow ({\bf k}).
\end{equation}
Second, the components of the  spinors in different layers are inverted with respect to each other, and thus the mass term  of the opposite sign in the two layers enters with the same sign in the matrices ${\cal D}_\uparrow ({\bf k})$ and ${\cal D}_\downarrow ({\bf k})$.
The dispersion relation in this case acquires the form
\begin{equation}
E_{\pm} = \sqrt{(\sqrt{k^2 + r^2} \pm \mu)^2 + |\Delta_{\bf k}|^2}. \label{disper_triplet_with_r}
\end{equation}
The masses in the two layers are of the opposite sign, due to the requirement of the PH symmetry of the whole system. Namely, under the transformation in each layer masses change sign [\onlinecite{son}], and if we, in addition, exchange layer index we recover the original Hamiltonian.

There are two important things to notice regarding the evolution of the CF state with increasing mass $r$:

(a) The minimum of the lower Bogoliubov band shifts from a finite value at $k_F^2 = \mu^2/(1 + \Delta^2)^2 - r^2$ to $k = 0$,
and this transition - without  closing of the gap - occurs at $r = \mu/(1 + \Delta^2)$;

  (b) Because   $k_F^2 = \mu^2/(1 + \Delta^2)^2 - r^2$, the Fermi momentum decreases with the mass, and therefore  the number of CF's reduces as the distance between the layers decreases.

Therefore the most important consequence of the assumed Dirac description of individual layers at large distances is that the number of CF's decreases as the distance between the layers decreases. For large distances we may assume that the pairing is weak, the order parameter
is small, and the pairing  cannot be detected then due to finite temperature effects, for instance. In any case we may choose, $ \Delta_{\bf k} = \Delta (k_x + i k_y)$, so that there is no Hall drag (pseudospin Hall effect) at large distances, but it develops gradually as the interlayer distance decreases and reaches  the  quantized value in agreement with experiments [\onlinecite{kim}]. This choice of the order parameter agrees with Refs. [\onlinecite{mdp,msr}].  For smaller distances $(r \sim \mu$ but $ r < \mu)$ we may assume that
the upper Bogoliubov band is pushed to high energies and an effective description in terms of quadratically dispersing CF's paired via weak $p$-wave pairing emerges,  implying an algebraically decaying Cooper pair wave function [\onlinecite{rg}].
The description of the system within this scenario then implies that at intermediate distances CB-CF mixture accounts for the total number of electrons [\onlinecite{srm}, \onlinecite{msr}]. As a consequence, composite bosons cannot have long range order, and likely have critical, algebraic pairwise correlations [\onlinecite{mdp}].

If at intermediate distances solely a collection of $p$-wave paired composite fermions, quadratically dispersing as in Ref. [\onlinecite{rg}], were present, signals of a topological phase with a large gap, $\Delta E \sim \mu$ would appear.
Instead, as detected on a torus in Ref. [\onlinecite{mdp}], there is an abundance of various low-energy excitations. This is in accordance with the above physical picture that implies a small portion of CF's at intermediate distances  in a topological phase with a small gap $\Delta E \sim \mu - r$, and $\mu \simeq r$.

As in the single layer case, anisotropic gapless solution\ (\ref{anisotropic_pairing}) is possible also for a bilayer. In the presence of a mass term $\sim r$ and in the case of the pairing  (\ref{anisotropic_pairing}) we obtain two anisotropic  Dirac cones at $ k_x^2 = \mu^2/(1 - \Delta^2) = k_0^2$ and $k_y = 0$. Expanding around $\pm k_0$ with $r \ll \mu$ we obtain  $ E^2 \approx (1 - 2 \Delta^2 - \frac{r^2}{\mu^2}) (\delta k_x)^2 + \Delta^2 (\delta k_y)^2$. The absence of a gap suggests a non-topological behavior. On the other hand, topological signatures were detected at intermediate distances in Ref. [\onlinecite{mdp}],
  in agreement with the characterization of isotropic weak $p$-wave pairing.  Thus the presence of the isotropic pairing, which may be accompanied by anisotropic ones,  seems crucial for the explanation of the properties at intermediate distances.

\section{Discussion and Conclusions}
The existence of anisotropic candidates for BCS paired states, in the case of monolayer (Sec. II), and bilayer (Sec. III and IV), is in agreement with the results in Ref. [\onlinecite{park}]. In that paper, the physics of the PH symmetric case of half-filled second Landau level is studied by exact diagonalization on a torus.  The main
 result of this numerical study is that the paired quantum Hall state in that case, as well closely related (by antisymmetrization) bilayer state, made of two kinds of electrons that each occupy quarter of the available single particle states in the second Landau level, are susceptible to anisotropic instabilities. By using the Dirac description of  the dipole nature of CFs, we can identify the paired quantum Hall state of Ref. [\onlinecite{park}] with PH Pfaffian, and its closeness to anisotropy as a sign of the relevance of  anisotropic solutions discussed in Sec. II. On the other hand, the relevance of the anisotropy for the bilayer state at effective $\nu=1/2 = 1/4 + 1/4$ total filling factor [\onlinecite{park}], may be again due to the composite - dipole nature of the CF's at filling factor $\nu=1/4$.
The Dirac description could be the easiest way to capture the dipole nature of a CF, despite the doubling of the fermionic degrees of freedom. In other words, we need particles and holes to describe dipoles [\onlinecite{sm}],  and the Dirac formalism could be a way to achieve that even in the cases when CF's have a Berry phase equal to $\pi/2$ (at quarter filling), with appropriate mass and chemical potential. If the Diracness is the cause of the anisotropic behavior, we can conclude that the Dirac formalism is equally applicable at $\nu=1/2$ and $\nu=1/4$. In this sense ``nothing is special at $\nu=1/2$" (Ref. [\onlinecite{fdmh}]) since only PH symmetry singles out Dirac description. The PH symmetry is sufficient but not necessary to cause the Diracness at the filling equal to one half.

If we restrict our discussion only to the case when CF's possess Berry phase equal to $\pi$, and thus Dirac formalism
seems appropriate for the bilayer case at total filling one, we demonstrated that the description by Dirac fermions is justified due to a global appearance and characterization of low-energy spectrum from the exact diagonalization on a torus [\onlinecite{mdp}]. In fact, the Dirac CF in the bilayer changes its Berry phase from value $\pi$ at large distances, to value $\sim 0$, at small distances (HLR fermion), while retaining its fermionic character. The second important consequence, due to the use of the Dirac formalism, is that the number of CF's is decreasing with the decreasing distance between layers. This is in in agreement with the necessity to use CF-CB mixed states to describe the bilayer at intermediate distances [\onlinecite{msr}].

Thus we can conclude that Dirac formalism can capture the basic phenomenology  of the bilayer at $\nu=1$, and the nature of the gapped paired states
in the single layer quantum Hall systems  with half-filled LL. We therefore expect it to become
an indispensable tool for further understanding of the paired states  in this context.

{\it Note added:} While this manuscript was in the final stage of preparation, Ref. [\onlinecite{wc}] appeared. It is a study of
possible pairings, based on Dirac formalism, and their realization in the case of a single layer with the half-filled LL.
 Ref. [\onlinecite{wc}] considered
pairings in the low-energy subspace of Dirac spectrum in
the context of a specific pairing mechanism.  In our work the low-energy projection is in place after the consideration of the pairing  instabilities within the Dirac formalism.  In this way we are able to account for the anisotropic pairings, with the consequences consistent with
theoretical and experimental findings, as we already emphasized.
\begin{acknowledgments}
We would like to thank S. Simon for a discussion.
The work was supported by
the Ministry of Education, Science, and Technological
Development of the Republic of Serbia under projects
ON171017 and  ON171031.
\end{acknowledgments}


text:
\section{Introduction}

Multistage sequential decision-making is commonly seen in many real-world applications for the sake of accuracy, efficiency, and cost-effectiveness, such as heathcare audit \cite{ekin2021integrated}, optimal investment \cite{sirbiladze2014multistage}, energy
\cite{gerking1987modeling}, optimization \cite{han2021multiple}, treatment \cite{tao2017adaptive}, and design of experiments \cite{mukkula2020robust}. Meanwhile, active learning is a special case of machine learning in which a learning algorithm can interactively query a user or some other information source to label new data points with the desired outputs \cite{settles2009active}, and has been integrated with decision-making in recent years in various contexts such as imbalanced data \cite{sundin2019active,lee2015Unbalanced} and Bayesian decision-making \cite{filstroff2021targeted}. In this paper, we aim at developing an active learning-based multistage sequential decision-making model to assist doctors making reliable diagnostic decisions and treatment recommendations in a cost-effective and convenient manner.

Our research is motivated from  the evaluation of gallstone. Gallstone are solid particles that can form from cholesterol, bilirubin, and other substances within the gallbladder. These densities are often benign when localized to the gallbladder, but can cause pain, infection, and liver damage when they become stuck in the common bile duct (CBD) and impede the flow of bile into the digestive tract. A stone that becomes impacted in the CBD can be difficult to detect definitively, but its presence requires a procedural intervention, endoscopic retrograde cholangiopancreatography (ERCP). As with most procedures, ERCP requires anesthesia and carries the risk of complications, including pancreatitis, infection, and bleeding. As a result, it is imperative to ensure that ERCP is only performed when there is a stone definitively obstructing the duct. This is particularly true in children \cite{zev2021creation}.
Thus, one important objective is how to efficiently predict the presence of a CBD stone in a child with high accuracy and specificity.

From the machine learning or statistical point of view, this might be straightforward:  one would like to collect as many informative features as possible for each patient, which will lead to higher prediction accuracy. However, from the medical point of view, this is highly non-trivial, since some informative features could be expensive and hard to obtain in many real-world scenarios.
Obtaining additional clinical data often comes with a price, including delayed time to intervention, cost, and complications. This price can be particularly high in pediatric care. For instance, while magnetic resonance imaging (MRI) is the most accurate way to diagnose a CBD stone, only 5\% of  pediatric  patients in our real dataset have such diagnostic, as pediatric MRI often involves sedation, long scheduling time, etc. Indeed a less cost-effective approach such as laboratory tests or ultrasound (US) often might be sufficient. Thus a common clinical challenge is how to reduce the cost of data collection for diagnosis decision-making, yet avoid computational issues associated with limited or incomplete data.

In this article, a multistage sequential decision-making model is developed to assist doctors making reliable diagnostic decisions and treatment recommendations with high accuracy in a cost-effective and convenient manner. In our approach, some patients are allowed to take fewer examinations if their symptoms are apparent according to early-stage features, and some patients will be recommended to take more examinations if their symptoms are ambiguous and doctors have difficulty making reliable treatment recommendations.

There are several benefits to the proposed model.
This novel platform allows the healthcare provider to actively collect only the necessary data sequentially rather than collecting all the data at inception. As a result, the diagnostic process can be faster, more convenient, and cost effective, while retaining significant accuracy.
In the research setting, our model offers a novel solution to address challenges raised by incomplete data through offering sequential decision-making in an active learning way.
The active learning-based model allows some missing data yet guarantees reliability and prediction accuracy, which makes the incomplete data missing ``at control" rather than missing at random.

In the proposed framework, one important thing is to identify a proper predictive model for each stage as the building block.
General categorization is required to utilize the model. Patients must be stratified to: certainly healthy (``0"), indeterminate (``0.5"), and certainly sick (``1").
After comparing the performance of different machine learning models, including neural networks, random forest, logistic regression, we selected the logistic regression model, which provides stability and accuracy with easy interpretability, see our previous work \cite{zev2021creation}.
In the sections to follow, a generalized logistic regression model will be utilized as a building block for each stage.

Multiclass logistic regression, also called  softmax/multinoulli/multinomial regression, is a classification method that generalizes logistic regression to multiclass problems ($K>2$). The idea behind this classification model is one-vs-all or one-vs-rest, which means $K$ binary classifier models will be trained to ensure every class has a binary classifier to identify whether a sample belongs to this specific class or not.
This model enables multiple outcomes but does not involve the ordering of the categorical dependent variable.
The ordinal logistic regression model, also called the ordered logit model or proportional odds model, overcomes this limitation by calculating the log of the odds of cumulative events \cite{mccullagh1980regression,kim2004topics,watson1986generalized}.
Unlike simple logistic regression, ordinal logistic models consider the probability of an event and all the events that are below the focal event in the ordered hierarchy.
The existing ordinal logistic regression model satisfies our needs in single stages but failed to deal with multistage scenarios simultaneously.
Therefore, in this research, we'll develop our own multistage model with the existing ordinal logistic regression model as a building block for each stage, but estimate the parameters in all stages together.

The contribution of our work is twofold. In terms of statistical analysis, an active learning model is proposed in a multistage decision-making process. This learning algorithm can interactively query a patient collect additional clinical data for collecting necessary information \cite{olsson2009literature}.
Also, the proposed model introduces a novel way of dealing with incomplete data by allowing some data to be absent and modeling with incomplete but necessary data directly, rather than guessing the incomplete data with uncertainty.
In terms of healthcare application,
a general strategy is developed to help the doctor making diagnostic decisions and treatment recommendations sequentially when some laboratory data are more expensive, invasive, or time-consuming to collect than other data.
By recommending requisite testing to those patients in whom symptoms do not support a definitive diagnosis, the proposed process is fast, cost-effective, yet reliable.
Also, the methodology applies to a large variety of disease diagnose that involve a sequence of indirect examinations and some unpleasant or inaccurate examinations.

Note that our proposed model/method is closely related to incomplete data analysis in the statistics literature, where many widely-used methods have been developed, e.g., simple deletion on feature or record, mean/median/mode/zeros substitution, regression imputation, last observation carried forward, EM, etc. \cite{kang2013prevention,ma2016missing}. However, these existing methods were developed from two primary approaches. One is deleting the records or features with missing values, which will lead to a waste of available data and some intrinsic problems associated with limited data size. The other approach is ``guessing" or ``learning" the missing data from the observed data, which can be dangerous for clinical modeling since the synthetic data brings extra uncertainties. Our approach is different from these existing methods in the sense that we actively decide which kinds of data to be observed and which kinds of data to be missing.

In Section 2, the proposed model is described.
Section 3 discusses the novel approach to parameter estimation for the proposed model.
In Section 4, the performance of the proposed methodology on synthetic data is reported. Finally, in section 5, the analytical result on real clinical data is demonstrated

\section{Proposed Active Learning Multistage Sequential Decision-Making Model}
In this section,
we first present a high-level framework of our proposed active learning-based multistage sequential decision-making model, and then introduce necessary notation to represent observed data in our framework. Finally, our model is rigorously developed based on the existing ordinal logistic regression model.

\subsection{Model Framework}

\begin{figure}[t]
   \centering
   \includegraphics[width = \textwidth]{figures/stagemodel.png}
    \caption{The framework of the proposed multistage sequential decision-making model, where $K$ is the number of stages, ${\pi}_i^{(j)}$ is the probability of patient $i$ having disease estimated in stage $j$, $U_j$ and $L_j$ is the cutoff points on the estimated probability in stage $j$, and $C_{K}$ is the single cutoff point in the last stage, $K$.}
   \label{fig_framework3}
\end{figure}

The proposed multistage sequential decision-making model is shown in Fig~\ref{fig_framework3}. To achieve the goal of optimizing the diagnostic process, under our proposed model, the features of a patient's clinical data are classified into several different stages according to the relevant cost, invasiveness, time, patient's willingness to participate, etc. The earlier stages involve features that are patient-friendly, inexpensive, though potentially less informative. The latter stages could involve some more time-consuming, expensive, unpleasant, but potentially more informative examinations. As in group sequential clinical trials, if we can make a definitive yes or no decision at a given stage, then we will stop taking features/observations and make the corresponding decision; but if we are not sure, then we will move to the next stages by taking additional features/observations from the patient.

To be more specific, for the $k$-th stage, based on the cumulative features, an ordinal logistic regression model is trained as the classifier to estimate the probability of sickness for each patient ($\pi_{i}^{(k)}$). For every classifier, the corresponding optimal cutoff points, $L_k, U_k$, or $C_K$ for the last stage, will be augment clinical decision making: if $\pi_i^{(k)} < L_k$, patient $i$ will be labeled as healthy (``0") and will not undergo further intervention; if $\pi_i^{(1)} > U_k$, patient $i$ will be labeled as sick (``1"), they will proceed to the recommended intervention, such as ERCP or surgery; otherwise,  the patient $i$ will be labeled as indeterminate (``0.5") and undergo further diagnostic workup.
Following the last stage, as no further diagnostic examinations are available, the provider must make a treatment recommendations based on a single threshold, $C_{K}$, and the patient will be classified as either healthy or sick.
By quantifying the probability of being sick via predictive models based on currently available features, the diagnosis could be more objective and stable compared with an empirical estimation.

In the evaluation of a newly presenting patient,
this framework with trained parameters can be utilized to arrive at a diagnosis and recommend treatments accordingly. In such a way, the necessary features will be collected sequentially and actively.

In the proposed multistage model, ordinal logistic regression is adopted as the basic predictive model in each stage. A novel feature of this model is, the patients' groups in the different stages have a nested structure rather than independent. Further, the features in the current stage $k$ include all those of the former stages and newly added features in the given category $k$.
This allows us to consider the predictive models in different stages together rather than in an isolated fashion. The details will be discussed further in the parameter estimation discussion in Section 3.

\subsection{Data Organization for 2-stage Model as Example}
To illustrate the proposed active learning style multistage sequential decision-making model, the 2-stage decision-making model will be elaborated as an illustrative example.
As discussed earlier, all the capturable features are first divided into several categories, each corresponds to a certain stage in our multistage model.
The criteria for feature grouping are flexible and can be adjusted according to real scenarios case-by-case:
1) For the features inside a category, the cost of collection and the information contained therein should be similar. For example, all indices can be obtained from one examination such as a blood test.

   2) Each category should contain some useful information, and excess categories should be avoided. For example, there is no need to treat each variable as a separate category,
   especially less informative variables

To illustrate the data organization for the proposed model, without loss of generality, we can take the 2-stage sequential decision-making model
as an example as follows.
In the observed (training) data, denote by $N_1$ and $N_2$  the number of patients in  stage 1 and stage 2, respectively, and denote by $p_1$ and $p_2$ the number of features collected in the two stages. In the first stage, the observed features or explanatory variables  are denoted by $\boldsymbol{X}^{(1)}= [x_{1}^{(1)}, x_{2}^{(1)}, \ldots, x_{p_{1}}^{(1)}]$ and the corresponding decision variable.  $Y^{(1)} \in \{``0", ``0.5", ``1"\}$, where $Y_{i}^{(1)} = ``0.5"$ implies that the more features need to be collected on patient $i$ in stage 2 for diagnosis.
In stage 2,  the new observed features or explanatory variables  are denoted by $(x_{1}^{(2)}, x_{2}^{(2)}, \ldots, x_{p_{2}}^{(2)})$ and the corresponding decision variable $Y^{(2)} \in \{``0", ``1"\}.$ Note that the feature set from stage 1 will be used in stage 2 cumulatively, and thus the number of features in stage 2 is actually $p_{1}+p_{2}$.

Hence, for the 2-stage model, the (training) data can be summarized as:
\begin{itemize}
    \item Stage 1: $\left\{Y_{i}^{(1)}; x_{i,1}^{(1)}, x_{i,2}^{(1)}, \ldots, x_{i, p_{1}}^{(1)}\right\} =\left\{Y_{i}^{(1)}; \boldsymbol{X}_{i}^{(1)}\right\}, i=1,2, \ldots, N_{1}, \\
    Y_{i}^{(1)} \in \{{``0", ``0.5", ``1"}\};$

    \item Stage 2: $\left\{Y_{i}^{(2)}; x_{i,1}^{(1)}, x_{i,2}^{(1)}, \ldots, x_{i,p_{1}}^{(1)}; x_{i,1}^{(2)}, x_{i,2}^{(2)}, \ldots, x_{i,p_{2}}^{(2)}\right\}=\left\{\boldsymbol{Y}^{(2)}; \boldsymbol{X}^{(2)}\right\}, i=1,2, \ldots, N_{2},
    Y_{i}^{(2)} \in \{{``0", ``1"}\};$
\end{itemize}

For a new patient, after observing $\boldsymbol{X}_{new}^{(1)} = [x_{new,1}^{(1)}, x_{new,2}^{(1)}, \ldots, x_{new, p_{1}}^{(1)}],$ we need to predict $Y_{new}^{(1)} \in \{{``0", ``0.5", ``1"}\}.$ When we predict $Y_{new}^{(1)} \in \{``0", ``1"\},$ this is the final decision. When we predict  $Y_{new}^{(1)} = ``0.5",$ we will then take more observations from this patient in the stage 2, i..e, $[x_{new,1}^{(2)}, x_{new,2}^{(2)}, \ldots, x_{new, p_{2}}^{(2)}],$ and make another prediction  $Y_{new}^{(2)} \in \{``0", ``1"\},$ which will be our final decision.

\subsection{Our Proposed Model}

In our proposed model, the building block is an ordinal logistic regression model, which can serve as a filter in each stage to split patients into 3 groups based on the ordered probability of being sick. For better presentation, let us first provide a quick review of ordinal logistic regression, which can be seen a form of a generalized linear model (GLM) and an extension of binary logistic regression.
Under ordinal logistic regression model, we assume that there is an underlying continuous latent variable
\begin{equation}
Y^{*}_{i}=\boldsymbol{\beta}^\top\boldsymbol{x}_{i}+\boldsymbol{\varepsilon}_{i},
\label{eq:ystar}
\end{equation}
where the error term $\boldsymbol{\varepsilon}_{i}$ has a logistic distribution with the cumulative distribution function (CDF) $F(t) = P( \boldsymbol{\varepsilon}_{i} \le t) = \frac{1}{1+e^{-t}}.$ Of course, other distributions of  $\boldsymbol{\varepsilon}_{i}$'s are possible such as Gaussian/normal distribution. Then the observable response variable $Y_i$ will fall into several ordinal categories based on the values of latent variable $Y^{*}_{i}$ and the cutoff points.

Now we are ready to present our proposed 2-stage sequential decision-making model. We assume that there are two latent variables, $Y^{(1)*}$ and  $Y^{(2)*}$ for the two stages as defined in (\ref{eq:ystar}). Then the decision at each stage is given by
\begin{equation}
Y_{i}^{(1)}=\left\{\begin{array}{ll}
0, & \text { if } Y^{(1)*}_i<L_{1}^{*} \\
0.5, & \text { if } L_{1}^{*} \leq Y^{(1)*}_i<U_{1}^{*} \\
1, & \text { if } Y^{(1)*}_i \geq U_{1}^{*}
\end{array}\right.
,
\quad
Y_{i}^{(2)}=\left\{\begin{array}{ll}
0, & \text { if }  Y^{(2)*}_i <C_{2}^{*} \\
1, & \text { if }  Y^{(2)*}_i \geq C_{2}^{*}
\label{eq:our2model}
\end{array}\right.\end{equation}
Note that the cutoff points $\left\{ L_{1}^{*}, U_{1}^{*}, C_{2}^{*} \right\}$ applying on the latent variable $Y^{*}$. This is equivalent to applying the cutoff points $\left\{ L_{1}, U_{1}, C_{2} \right\}$ on the probability estimation $\pi_i^{(k)}$ in the Fig~\ref{fig_framework3} after going through a link function transformation.

Recall
that in our proposed multistage model, the output $Y_{i}^{(1)} \in \{{``0", ``0.5", ``1"}\}$ at the first stage represents healthy, indeterminate, and sick, respectively. At the last stage,  $Y_{i}^{(2)} \in \{{``0", ``1"}\}$ indicates healthy and sick. It is useful to mention that the outputs are ordered, i.e. $Y_i^{(1)} \leq Y_{i}^{(2)}.$

A subtlety in our proposed multistage model is the relationship of the parameter $\boldsymbol{\beta}$ in (\ref{eq:ystar}) across different stages.
A computationally simple baseline approach is to assume that they can be completely different, where we can simply apply ordinal logistic regression models to estimate parameters each stage individually and independently. Unfortunately, this does not consider the nested structure and potential correlations between the two adjacent stages, and thus might cause mis-understanding and confusing interpretation in practice. In our paper, we will make a further assumption that $\boldsymbol{\beta}^{(1)}$ and $\boldsymbol{\beta}^{(2)}$ enjoys the nested structure in the sense of having the same $\beta$ values with the same features, i.e.,
\[
\boldsymbol{\beta}^{(2)} =\left\{ \beta^{(1)}_1,\beta^{(1)}_2, \ldots, \beta^{(1)}_{p_1}; \beta^{(2)}_1, \beta^{(2)}_2, \ldots, \beta^{(2)}_{p_2} \right\}= \left\{ \boldsymbol{\beta}^{(1)}; \beta^{(2)}_1, \beta^{(2)}_2, \ldots, \beta^{(2)}_{p_2} \right\}.
\]
Under this common coefficients assumption, the newly added features from additional stages will not affect the weights of existing ones, which makes the model more interpretable for doctors and patients.

\section{Parameter Estimation}

At the high-level, we propose to use the maximum likelihood estimation (MLE) method to estimate not only the unknown $\beta$'s values in different stages, but also the cutoff values $\left\{ L_{1}^{*}, U_{1}^{*}, C_{2}^{*} \right\}$. Note that the coefficients $\beta$'s  can also be estimated by the penalized likelihood estimation framework such as LASSO for variable selection when there are a large of available features, which is beyond the scope of this article and will be presented elsewhere. Here we assume that we pre-screen all features, and only use those useful features in our model and medical decision making.

To better present our main ideas, let us begin with a single stage with a ordinal logistic regression model. To simplify the notation, we assume that the response $Y$ results from the underlying latent variable $Y^*$ by using cutoff points $\theta_{1} < \theta_{2}$:
\begin{equation}
Y=\left\{\begin{array}{lr}
0, \text { if }  \theta_{0} < Y^{*}  \leq \theta_{1} \\
0.5, \text { if } \theta_{1}<Y^{*} \leq \theta_{2} \\
1, \text { if } \theta_{2}<Y^{*} \leq \theta_{3} \\
\end{array}\right. ,
\label{eq:olr}
\end{equation}
where $\theta_0 = -\infty$ and $\theta_3 = \infty$ are introduced for the technical reason to simplify the notation.  Also note that at the last stage, we have $\theta_{1} = \theta_{2}$, which will force one to make a decision $\{Y=0\}$ or $\{Y=1\}.$ Under this notation,  the conditional distribution of $Y$ is given by
\begin{equation}
\begin{aligned}
P(Y=j \mid \boldsymbol{x}) &=P\left(\theta_{j-1}<Y^{*} \leq \theta_{j} \mid \boldsymbol{x}\right)
=P\left(\theta_{j-1}<\boldsymbol{\beta}^{\top}  \boldsymbol{x}+\varepsilon \leq \theta_{j}\right) \\
&=F\left(\theta_{j}-\boldsymbol{\beta}^{\top}  \boldsymbol{x}\right)-F\left(\theta_{j-1}-\boldsymbol{\beta}^{\top}  \boldsymbol{x}\right)
\end{aligned}
\label{eq:probn}
\end{equation}
where $\boldsymbol{\beta}, \theta_{j}$ are model coefficients and cutoff points to be estimated from the data, $P$ is the abbreviation of probability, and $F$ is the CDF of the error term $\boldsymbol{\varepsilon}$ which takes on the role of the inverse link function.

For a given stage or a given ordinal logistic regression model with $N$ subjects, the joint log-likelihood function becomes
\begin{equation}
\log \mathcal{L}\left(\boldsymbol{\beta}, \boldsymbol{\theta} \mid \boldsymbol{X}, Y\right)=
\sum_{i=1}^{N}
\sum_{j=1}^{J} I(y_{i}=j) \log \left[ F\left(\theta_{j}-\boldsymbol{\beta^{\top}} \boldsymbol{x}_{i}\right)-F\left(\theta_{j-1}-\boldsymbol{\beta^{\top}} \boldsymbol{x}_{i}\right) \right]
\notag
\end{equation}
It has been well investigated to maximize this log-likelihood objective function which possesses several nice properties.
McCullagh presented a Fisher scoring algorithm for ML estimation, expressing the likelihood function using cumulative probabilities. McCullagh showed that a sufficiently large $n$ guarantees a unique maximum of the likelihood \cite{mccullagh1980regression,agresti2003categorical}.
Burridge and Pratt showed that the log-likelihood is concave for many cumulative link models, including the logit, probit, and complementary log-log \cite{agresti2003categorical}. Iterative algorithms usually converge rapidly to the ML estimates.
The optimization computing algorithms are developed based on these useful properties. The existing approaches compute the MLE by iteratively reweighted least squares (IRLS). Actually, for packages in $R$ software, all GLM are fit using the IRLS  estimation method for the log-likelihood.

Now let us go back to the parameter estimation of our proposed multistage model. For simplicity,  denote $[\boldsymbol{\beta}, \boldsymbol{\theta}]$ as $\boldsymbol{\zeta}$, and $\log \mathcal{L}\left( \boldsymbol{\zeta} \mid \boldsymbol{X}, Y\right)$ as $\ell(\boldsymbol{\zeta})$.
Under our notation, the joint log-likelihood function for all samples in all stages will be as follows:
\begin{equation}
\begin{aligned}
\ell(\boldsymbol{\zeta})=
\sum_{k=1}^{K}
\sum_{i=1}^{N_k}
\sum_{j=1}^{J_{k}} I(y_{i}^{(k)}=j) \log \left[ F\left(\theta^{(k)}_{j}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right)-F\left(\theta^{(k)}_{j-1}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right) \right]
\notag
\end{aligned}
\end{equation}
where
$K$ is the number of stages in the proposed model, $N_{k}$ is the number of patients in the $k^{th}$ stage, and $J_k$ is the number of categories in the $k^{th}$ stage for which we have
\begin{equation}
J_{k}=\left\{\begin{array}{ll}
3, \text { if } 1 \leq k < K \\
2, \text { if } k=K
\end{array}\right.
\notag
\end{equation}

Given the joint likelihood function $\ell(\boldsymbol{\zeta})$, the next step is to simultaneously estimate $\boldsymbol{\zeta}$, all the parameters in the $K$ stages, that maximizes the likelihood function.
Here we propose to adopt the gradient descent algorithm and a key step is to compute the gradient for each parameter, which can be written as:
\begin{equation}
\frac{\partial \ell(\boldsymbol{\zeta})} { \partial \boldsymbol{\beta}^{(k)}}=
\sum_{i=1}^{N_{k}}
\sum_{j=1}^{J_{k}} I(y_{i}^{(k)}=j) \frac{ (-\boldsymbol{x}_{i}^{(k)}) \left[ f\left(\theta_{j}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right) - f\left(\theta_{j-1}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right) \right] }
{  F\left(\theta_{j}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right)-F\left(\theta_{j-1}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right) }
\notag
\end{equation}
\begin{equation}
\begin{aligned}
\frac{\partial \ell(\boldsymbol{\zeta})} { \partial \theta_{j}^{(k)}}=
\sum_{i=1}^{N_{k}}
\sum_{j=1}^{J_{k}} I(y_{i}^{(k)}=j)&\{
\frac{ f\left(\theta_{j}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right)}
{  F\left(\theta_{j}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right)-F\left(\theta_{j-1}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right) }  \\
&-\frac{ f\left(\theta_{j}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right)}
{  F\left(\theta_{j+1}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right)-F\left(\theta_{j}^{(k)}-\boldsymbol{\beta}^{(k)\top} \boldsymbol{x}_{i}^{(k)}\right)  }
\}
\end{aligned}
\notag
\end{equation}
where
$f(t) = F'(t)$ is the probability density function (PDF) of  $\boldsymbol{\varepsilon}$:

 f(t) = \frac{1}{(e^{t/2} + e^{-t/2})^2}$.

We should also remark some technical details in computation.
Note that the coefficient in stage $k$, $\boldsymbol{\beta}^{(k)}$, contains $\boldsymbol{\beta}^{(k-1)}$ and newly added coefficients for features in the $k^{th}$ stage. This makes the coefficients in different stages of different lengths. To address this issue, we propose finding the gradient in each stage $\frac{\partial \ell(\boldsymbol{\zeta})} { \partial \boldsymbol{\beta}^{(k)}}$ individually first. Then, adding them together with the front ends aligned and pad the shorter ends with 0's:
$\frac{\partial \ell(\boldsymbol{\zeta})} { \partial \boldsymbol{\beta}} = \left[\frac{\partial \ell(\boldsymbol{\zeta})} { \partial \boldsymbol{\beta}^{(1)}}; 0,0,...,0 \right] + ... + \left[ \frac{\partial \ell(\boldsymbol{\zeta})} { \partial \boldsymbol{\beta}^{(K)}}\right]
$.
In addition, the indicator function $I(y_i^{(j)} = j)$ can be replaced with a mask matrix with boolean values. Alternatively, the indicator function can be avoided by data preprocessing.
For instance, we can split the patients into different groups to mark the identity, e.g., using $G_{k,j}$ to denote the group of patients who have label $j$ in stage $k$, see \textit{Appendix A} for more details.

To numerically solve our proposed joint-optimization problem, we follow the framework of the existing $polr()$ function in the package MASS in R software, and adopt the $optim()$ function as optimization solver.
The $optim()$ is a function designed for general-purpose optimization based on Nelder--Mead, quasi-Newton and conjugate-gradient algorithms.
In $optim()$, the Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) method, an iterative method for solving unconstrained nonlinear optimization problems \cite{avriel2003nonlinear}, is adopted, since  it works best with analytic gradients based on our extensive experiences.

\section{Simulation Study}
In this section, a simulation is conducted to verify the effectiveness and efficiency of the proposed methodology.
First, a 2-stage synthetic dataset, $\{X^{(1)}, Y^{(1)}; X^{(2)}, Y^{(2)} \},$ is generated with prescribed coefficients $\boldsymbol{\beta}^{(1)}, \boldsymbol{\beta}^{(2)}$ and classification thresholds $\{ L^{*}_1, U^{*}_1, C^{*}_2\}$.
Then,
we conduct both a baseline approach
and the proposed approach on the synthetic dataset
to demonstrate the performance of the proposed method in the aspect of prediction accuracy, stability, interpretability, and relative efficiency

Let $N_{1}=10,000$ denote the total number of synthetic patients at the outset.
We design the features of patients as follows:
$X_1 \sim Bernoulli (0.3)$, $X_2 \sim N (-1, 1)$, $X_3 \sim N (1, 1)$, $X_4 \sim N (0, 2)$, $X_5 \sim Bernoulli (0.4)$, $X_6 \sim N (-1, 1)$, $X_7 \sim N (0, 1)$.
For stage 1, we design feature set as: $\boldsymbol{X}^{(1)}=\left\{ X_1, X_2, X_3, X_4 \right\}$.
For stage 2, we add 3 new features: $\boldsymbol{X}^{(2)}= \left\{ X_1, X_2, X_3, X_4; X_5, X_6, X_7 \right\}= \left\{\boldsymbol{X}^{(1)}; X_6, X_7, X_8 \right\}$.
The error term $\boldsymbol{\varepsilon}$ follows standard logistic regression for both stages:
$\boldsymbol{\varepsilon^{(1)}} \sim Logistic(0, 1), \quad \boldsymbol{\varepsilon^{(2)}} \sim Logistic(0, 1)$.

Note that if the error term $\boldsymbol{\varepsilon}$ follows logistic regression with other parameters,
the estimation of $\boldsymbol{\beta}$ and $\boldsymbol{\theta}$ will be scaled,
see Eq.~(\ref{eq:our2model}).
The coefficients are designed as: $\boldsymbol{\beta}^{(1)} = [2, 2, 2, 2], \quad \boldsymbol{\beta}^{(2)} = [\boldsymbol{\beta}^{(1)}, 4, 4, 4] = [2, 2, 2, 2, 4, 4, 4] $
which satisfys the common coefficients assumption.
Next, we can calculate the latent variable $Y^*$ according to Eq.~(\ref{eq:ystar}).
To cut the $Y^*$ into ordered categories, we design the cutoff points as $\boldsymbol{\theta}= \{ L_1^{*}, U_1^{*}, C_2^{*} \} = \{-2.2, 2.2,  0.5\}$.
The value $-2.2$ and $2.2$ are chosen since $U_1 = F(2.2) = 0.9$ and $L_1 = F(-2.2) = 0.1$, where $0.9$ and $0.1$ are the designed cutoff points on probability $Y$. Via such design,
we have $N_{2}  \approx \frac{1}{3}N_{1}$.
Given $\boldsymbol{\theta}$, the corresponding true label $Y$ can be obtained according to Eq.~(\ref{eq:our2model}).

In the baseline approach, the 2 stages proceed independently: in stage 1, a standard ordinal logistic regression is applied on $\{\boldsymbol{Y}^{(1)},\boldsymbol{X}^{(1)}\}$ via $polr()$ function in R; and in stage 2, a standard logistic regression is applied on $\{\boldsymbol{Y}^{(2)},\boldsymbol{X}^{(2)}\}$ via the $glm()$ function in R.
While in the proposed method, all the parameters in 2 stages are estimated together via MLE.
To provide reliable statistical inference, we adopt the Monte-Carlo method and run the data generation and parameter estimation 100 times

\begin{figure}[t]
    \centering
    \includegraphics[width =0.7\linewidth]{figures/meancoefsimu.png
    \caption{Mean coefficients estimation comparison among ground truth (black line), mean coefficients estimated from baseline approached individually (blue and green line), and mean coefficients estimated from proposed method jointly (red line). This plot demonstrates that both methods have little or no bias when estimation parameters, proving the effectiveness of both methods on estimating.}

    \label{fig:coefsimu}
\end{figure}

The comparison of the mean coefficients estimated from the two methods is visualized in Fig. ~\ref{fig:coefsimu}. Clearly, the mean coefficients estimated from the two methods are similar and both coincide well with ground truth, indicating that both methods have little or no bias on mean parameter estimation, proving the effectiveness of both methods on estimating.
To quantify the subtle difference,
we can utilize  mean square error (MSE) to measure the distance between the estimated coefficients and the ground truth.
For coefficients in stage 1, the MSE of the estimation from the baseline method and the proposed method  are $7.62 \times 10^{-6}$ and $2.25 \times 10^{-6}$, respectively. Meanwhile, for stage 2, the MSE of the baseline method and the proposed method are 0.0011 and 0.0003, respectively. This indicates a slight advantage of the proposed methods on mean coefficients estimation. See Appendix B for more comparison on mean prediction metrics.

Fig.~\ref{fig:stdcoefsimu} shows the standard deviation (std) of coefficients estimation from the two methods. It can be seen that the standard deviation of estimation from the proposed method is significantly smaller than those from baseline method in both stages.
To quantify this advantage, the relative efficiency of the proposed method with respect to the baseline method, defined as $e(T_1, T_2) = \frac{var(T_{2})}{var(T_{1})}$, is plotted in Fig. ~\ref{fig:relativeEfficiency}, where $T_1$ and $T_2$ represent the estimators for $\hat{\boldsymbol{\beta}}$ in the proposed method and baseline method, respectively.
It can be seen that the relative efficiency of the estimator in the proposed method for the coefficients in stage 1 is very big, as high as 19.38. And for the estimator for the coefficients in stage 2, the lowest relative efficiency is still $1.62$. This indicates that the efficiency of the parameter estimation is improved $62\% ~ 1838\%$ in the proposed methodology.
The reason for such improvement is straightforward: in the proposed methods, all available data is involved in the parameter estimation, while in the baseline method, only a small subset of the data is involved in later stages.

\begin{figure}[t]
\centering
\begin{minipage}[t]{0.48\linewidth}
\centering
\includegraphics[width=6cm]{figures/stdcoefsimu.png}
\caption{Standard deviation of coefficients estimation from two methods.
 \label{fig:stdcoefsimu}
\end{minipage}
\begin{minipage}[t]{0.48\linewidth}
\centering
\includegraphics[width=6cm]{figures/relaEffi.png}
\caption{Relative efficiency of the proposed method with respect to the baseline method.
\label{fig:relativeEfficiency}
\end{minipage}
\end{figure}

For the estimation of cutoff points, the baseline method and the proposed method give the mean estimations as $\{-2.201, 2.192,  0.536\}$ and $\{-2.199, 2.197,  0.506\}$, respectively.
And the standard of deviation of the estimations by the baseline and proposed methods are $\{0.061, 0.064,  0.178\}$ and $\{0.057, 0.062, 0.119\}$, respectively.
Compared with the ground truth, $ \{ \boldsymbol{\theta}\}_{ground ~ truth} =\{-2.2, 2.2, 0.5\}$,
both methods are functional, yet the proposed approach yields a lower mean deviation and standard deviation.

Overall, the proposed 2-stage model provides a more accurate and  efficient estimation since it involves all data in the parameter estimation.
Also, in implementation, the proposed algorithm enjoys an advantage with regard to computational convenience since we estimated all parameters in $\zeta$ at once. Additionally, the common coefficients assumption in the proposed method reduces the number of parameters to be estimated from $4+7+2+1 = 14$ to $7+2+1 = 10$ since the coefficients for features in stage 1 are not repeatedly estimated in later stages, but instead assumed to be consistent.

\section{Case Study}

In this section,  we analyze a real data set in \cite{zev2021creation} and demonstrate the usefulness of our proposed approach in the sense of being able to learn from data to effectively evaluate gallstone. In particular, the real clinical data is modeled as a 2-stage sequential decision-making model.
We have $N_{1} = 316$ patients in the first stage. Based on the missing profile, we filter the patients who have features in category 2 into the second stage and label them as $Y_{i}^{(1)}=0.5, i = 1,..., N_{2}$. As a result, $N_{2} = 189$.
Note that here we made the assumption that if a patient has available data in the next stage, the patient should be labeled as ``0.5" at the current stage. The reason for such assumption is that we believe that the provider only recommends further diagnostic interventions to a patient who has indeterminate presentation at the current stage (``0.5"). In other words, the real data comes from the provider's knowledge and experience, and the missing pattern is not missing at random but deliberately interfered with.
Note that features in $\boldsymbol{X}^{(1)}$ include hemolytic disease, gender, age, body mass index (BMI), total bilirubin, alanine transaminase (ALT), aspartate transaminase (AST), lipase, amalyse, gamma-glutamyl transferase (GGT), alkaline phosphatase. These features include demographic and laboratory tests. For features in $\boldsymbol{X}^{(2)}$, based on $\boldsymbol{X}^{(1)}$, it involves  imaging data, including CBD diameter and the presence of CBDS by ultrasound (US) and magnetic resonance imaging (MRI) examinations, respectively. See Table \ref{tbl:CBDS} in Appendix C for descriptive statistics
 \begin{figure}[t]
\centering

        \includegraphics[width=0.8\linewidth]{figures/coef32.png}

    \caption{Coefficients comparison between coefficients obtained from 2 methods for features in stages 1.

    }
    \label{fig:coef}
\end{figure}

Fig. ~\ref{fig:coef} visualizes the coefficients for features in category 1 estimated by baseline and proposed approach, respectively.
It can be seen that with the baseline approach, the coefficients estimated by the two individual models are similar for most features, $X_{2} \sim X_{11}$, but significantly different for $X_{1}$.
In the proposed 2-stage model, the coefficients estimated for feature $X_{1} \sim X_{11}$ seem to resemble a weighted average of the coefficients obtained from baseline individual models.
Also, compared with the baseline method, the coefficients estimated by the proposed method have no dramatic spikes. Also, with the common coefficients assumption, newly added features will not affect the weights of existing ones. More importantly, in the proposed method, no matter how many stages a provider utilized, the calculation can be done at once efficiently.

Our model and method can be useful for providing the treatment recommendation for a new patient, and a diagnostic mobile App was developed in 2021 in Apple Store based on this research. To be concrete, when a new patient presents, given the trained proposed method, the provider (e.g., clinician or doctor), can first recommend the patient undergo diagnostic tests included in category 1, and calculate the $\hat{Y}^{(1)}$ accordingly. Based on the $\hat{L}_1$ and $\hat{U}_1$, the patient will be labeled as healthy (``1"), sick (``0"), or indeterminate and proceed to the next stage (``0.5"). Then, only when the patient's presentation is indeterminate,  additional testing such as ultrasound will be undertaken.
In other words, only the necessary data are actively collected sequentially for sequential decision-making.
This approach simplifies and expedites the diagnostic workup.
Additionally, the provider can adjust the multistage sequential decision-making model flexibly by designing the number of stages and the features in each stage, also see Section 2 for the details on the feature grouping recommendations.

\section{Conclusion and Discussion}
In this study, to augment healthcare clinical decision-making based on limited examinations, a multistage sequential decision-making model is developed to actively collect only minimal necessary diagnostic data
To illustrate the proposed method, we optimized a 2-stage model on synthetic data in the simulation study and applied the proposed algorithm to pediatric patients with suspicion for obstructive common bile duct stones in an effort to more efficiently and accurately arrive at a diagnosis that allowed expedited intervention with ERCP.
In both the simulation study and real case study on collected pediatric data, the proposed method is more effective, stable, interpretable, and simpler in estimating the coefficients and cutoff points for all stages.
It estimates all parameters simultaneously
and fully considers the cumulative characteristics between these stages by making common coefficients assumption.

One limitation of the proposed methodology is that it involves all variables without selection.
A potential future improvement to the model will be the addition of  $L_{1}$ regularization term to the objective function to realize the variable selection.

Also, we can extend our proposed algorithm to satisfy more real-world cases with special requirements. For example, consideration of expense, delay to intervention, and need for invasive testing can be specified and minimized with modification to the objective function as
$\min _{N_{1}, N_{2}, \zeta}-\log \mathcal{L}(\zeta \mid X, Y)+\lambda\left(c_{1} N_{1}+c_{2} N_{2}\right)$,
where the $c_1, c_2$ are the relevant cost per patient in stage 1 and stage 2, respectively. $\lambda$ is the tuning parameter to adjust the weight between the negative likelihood and the relevant cost to the test. Note that here the $N_1, N_2$ can also be variables to be optimized over since the relevant cost depends on the number of patients in each stage.

In healthcare settings, the penalty for misdiagnosis can vary based on the patient and diagnosis. For example, false-positive detection may lead to unnecessary procedures or excess intervention, while false-negative detection may result in missed or delayed treatment. For such cases with weighted penalties on different types of misclassification, we can modify the objective function by specifying the penalties as
$\min _{\zeta}-\log \mathcal{L}(\zeta \mid X, Y)-\lambda_{1} \sum_{i} I\left(\left(y_{i}=0.5 \mid 1\right), \hat{y}_{i}=0\right)-\lambda_{2} \sum_{i} I\left(\left(y_{i}=0 \mid 0.5\right) \& \hat{y}_{i}=1\right)$,
where the $\lambda_1$ and  $\lambda_2$ are the tuning parameter for false-negative prediction and false-positive prediction, respectively.

The proposed algorithm is also applicable to a broad range of decision-making scenarios outside of healthcare as it utilizes active learning to drive the sequential collection of minimal requisite information, rather than all possible information simultaneously. This will be an area of future development.

redpajama_set_name: RedPajamaArXiv