\section{Supplemental material}
\subsection{A. Proof of Theorem~\ref{main_thm}}

\begin{lemma}\label{monomial}
    For $g \in \mathcal{M}_m$,
    \begin{equation}
        [B_n(g)]_{\vb*{n'}\vb*{n}}=\delta_{\sigma_g(\vb*{n'}),\vb*{n}}  \prod_{i=0}^{m-1} [D(g)]_{ii}^{n_i}.
        \label{monomial_formula}
    \end{equation}
\end{lemma}
\begin{proof}
Because $B_n$ is a group homomorphism, $B_n(g)= B_n(P_{\sigma_g})B_n(D(g))$.
From the definition,
\begin{equation}\label{permutation_formula}
    B(P_{\sigma_g})\ket*{\vb*{n}}=\prod_{i}\frac{1}{\sqrt{n_i!}}\qty\Big(\sum_{j}(a^\dag_j [P_{\sigma_g}]_{ji}))^{n_i}\ket*{\vb*{0}}
    =\prod_{i}\frac{1}{\sqrt{n_i!}}(a^\dag_{\sigma_g(i)})^{n_i}\ket*{\vb*{0}} = \ket*{{\sigma_g}^{-1}(\vb*{n})}
\end{equation}
and
\begin{equation}\label{diagonal_formula}
    B(D(g))\ket*{\vb*{n}}=\prod_{i}\frac{1}{\sqrt{n_i!}}\qty\Big(\sum_{j}(a^\dag_j [D(g)]_{ji}))^{n_i}\ket*{\vb*{0}}
    =\prod_{i}[D(g)]_{ii}^{n_i}\ket*{\vb*{n}}.
\end{equation}
Therefore, we have Eq.~\eqref{monomial_formula}.
\end{proof}

\begin{lemma}\label{projector}
Let
\begin{equation}
    \rho^G_\lambda = \frac{1}{\abs{G}}\sum_{g\in G} \lambda(g)^{-1} B(g)
\end{equation}
for an Abelian group $G$. Then, $\rho^G_\lambda$ is the projector onto $V^G_\lambda$.
\end{lemma}
\begin{proof}
For any state $\ket*{\Psi} \in V^G_\lambda$, $\rho^G_\lambda \ket*{\Psi}=\ket*{\Psi}$ holds.
It also holds that
\begin{align}
    B(g) \rho^G_\lambda &= \frac{1}{\abs{G}}\sum_{h\in G} \lambda(h)^{-1} B(gh)
    =\frac{1}{\abs{G}}\sum_{h\in G} \lambda(g^{-1}h)^{-1} B(h) =
    \lambda(g)\rho^G_\lambda
\end{align}
for any $g \in G$.
Therefore, $\rho^G_\lambda \ket*{\Psi} \in  V^G_\lambda$ for any state $\ket*{\Psi}$, and $(\rho^G_\lambda)^2 =\rho^G_\lambda$ .
\end{proof}
Lemma~\ref{projector} can be regarded as a special case of the projection formula in representation theory~\cite{Fulton1991-gz}, where the considered group is Abelian.
\begin{proof}[Proof of Theorem~\ref{main_thm}]
(i) We have
\begin{align}
    \norm{\rho_\lambda^G \ket*{\Psi}}^2 &= \sum_{\vb*{n},\vb*{n'}}  c_{\vb*{n'}}^*  c_{\vb*{n}} \bra{\vb*{n'}} \rho_\lambda^G \ket*{\vb*{n}}
    = \sum_{\vb*{n},\vb*{n'}}  c_{\vb*{n'}}^*  c_{\vb*{n}} \frac{1}{\abs{G}}\sum_{g\in G}\lambda(g)^{-1} \bra{\vb*{n'}} B(g) \ket*{\vb*{n}} \notag \\
    &= \frac{1}{\abs{G}}\sum_{g\in G} \lambda(g)^{-1}\sum_{\vb*{n}}c_{\sigma^{-1}_g(\vb*{n})}^* c_{\vb*{n}}\prod_{i=0}^{m-1} D(g)_{i,i}^{n_i},
\end{align}
where the second and third equalities hold from Lemmas~\ref{projector} and~\ref{monomial}, respectively.

(ii) If $\ket*{\Psi_\text{in}}\in V_\lambda^G$ holds, $\ket*{\Psi_\text{out}}=B(U)\ket*{\Psi_\text{in}}$ satisfies
\begin{equation}
    B(UgU^\dag)B(U)\ket*{\Psi_\text{in}}=B(U)B(g)\ket*{\Psi_\text{in}}=\lambda(g)B(U)\ket*{\Psi_\text{in}}
\end{equation}
for any $g\in G$.
Thus, letting $G'=UGU^\dag$ and $\lambda'(UgU^\dag)=\lambda(g)$ for $g \in G$, it holds that  $\ket*{\Psi_\text{out}} \in V^{G'}_{\lambda'}$.
The converse also holds.

(iii) $\bra{\vb*{n}}\ket*{\Psi}=\bra{\vb*{n}}\rho_\lambda^G \ket*{\Psi}=0$ if $\rho_\lambda^G\ket*{\vb*{n}}=0$.
For any $g\in G$,
\begin{equation}
    \lambda(g)\bra{\vb*{n}}\ket*{\Psi}=\bra{\vb*{n}}B(g)\ket*{\Psi}=B(D(g))_{\vb*{n},\vb*{n}}\bra{\sigma_g(\vb*{n})}\ket*{\Psi},
\end{equation}
where Lemma~\ref{monomial} is applied.
Thus, $\abs{\bra{\vb*{n'}}\ket*{\Psi}}=\text{const.}$ for any $\ket*{\vb*{n'}}\in \mathcal{O}_{\vb*{n}}$.
\end{proof}
\begin{proof}[Proof of Corollary~\ref{interference}]
Because $\ket*{\Psi_\text{out}} \in V^{G'}_{\lambda'}$ from Thm.~\ref{main_thm}(ii), Thm.~\ref{main_thm}(iii) yields $\bra{\vb*{n}}\ket*{\Psi_\text{out}}=0$ when $\rho_\lambda^G\ket*{\vb*{n}}=0$.
This condition is equal to Eq.~\eqref{condition} from Thm.~\ref{main_thm}(i), when $\ket*{\Psi}=\ket*{\vb*{n}}$ and $G$ consists of diagonal matrices.
\end{proof}
\begin{proof}[Proof of Corollary~\ref{stab_meas_formula}]
Note that any Fock state is a joint eigenstate of $B(G')$ when $G'$ consists of diagonal matrices, as seen from Eq.~\eqref{permutation_formula}.
Then, the projector can be expressed as $\rho^{G'}_{\lambda'} = \sum_{\ket*{\vb*{n}}\in V^{G'}_{\lambda'}} \dyad{\vb*{n}}$.
From the definition, the probability of obtaining the measurement results $g=\lambda(g)$ for each $g\in G$ is
\begin{equation}
    \sum_{\ket*{\vb*{n}}\in V_{\lambda'}^{G'}} \abs{\bra{\vb*{n}} B(U)\ket*{\Psi}}^2,
\end{equation}
which is equal to
\begin{equation}
    \bra{\Psi} B(U^\dag) \rho^{G'}_{\lambda'} B(U)\ket*{\Psi}=\norm{\rho^{G}_{\lambda} \ket*{\Psi}}^2.
\end{equation}
\end{proof}

\subsection{B. Proof of Theorem~\ref{bell_thm}}
\begin{definition}\label{Dicke}
    The Dicke states are
    \begin{equation}
        \ket*{D^m_{k}} = \binom{m}{k}^{-\frac{1}{2}} \sum_{\substack{x\in \{0,1\}^{m}\\ wt(x)=k}} \ket*{x},
    \end{equation}
    where $wt(x)$ is a Hamming weight of $x$.
\end{definition}
\begin{lemma}\label{probability}
    Assume an Abelian transitive subgroup of symmetry group $\Sigma \subset S_m$, where there exists unique $\sigma \in \Sigma$ such that $\sigma(0)=i$ for each $i\in \{0,\cdots, m-1\}$.
    Letting ${P}_j(\theta)$ be the operator acting on the $j$th modes as ${P}_j(\theta)\ket*{{k}}_j=e^{ik\theta}\ket*{k}_j$,
    \begin{equation}\label{prob_Dicke_formula}
        \norm*{\rho^{\Sigma}_1 {P}_0(\theta) \ket*{D^m_{k}}}^2 = \abs{\frac{(m-k)+e^{i\theta}k}{m}}^2,
    \end{equation}
    where $\rho^{\Sigma}_1$ is the projector onto the joint eigenstate of $\{B(P_\sigma)|\sigma \in \Sigma\}$ with eigenvalues of $1$.
\end{lemma}
\begin{proof}
From the definition, the Dicke states are invariant under any permutation of the modes; that is,
\begin{equation}
    B(P_\sigma)\ket*{D^m_{k}}=\ket*{D^m_{k}}
\end{equation}
for any $\sigma \in S_m$.
Further, because ${P}_i(\theta)$ acts on the $i$th modes,
\begin{equation}
    B(P_\sigma){P}_i(\theta)B(P_\sigma)^\dag={P}_{\sigma(i)}(\theta),
\end{equation}
from Eq.~\eqref{permutation_formula}.
From these equations and Lemma~\ref{projector},
\begin{equation}\label{cal1}
    \rho^{\Sigma}_1 {P}_0(\theta) \ket*{D^m_{k}} = \frac{1}{\abs{\Sigma}}\sum_{\sigma\in \Sigma} {P}_{\sigma(0)}(\theta) \ket*{D^m_{k}}.
\end{equation}
The group $\Sigma$ being an Abelian transitive group implies that there is only one element $\sigma \in \Sigma$ such that $\sigma(0)=i$ for each $i\in \{0,\cdots, m-1\}$, and that $\abs{\Sigma}=m$.
Therefore, Eq.~\eqref{cal1} is equivalent to
\begin{equation}
    \frac{1}{m}\sum_{i=0}^m {P}_i(\theta) \ket*{D^m_{k}} = \frac{(m-k)+e^{i\theta}k}{m}\ket*{D^m_{k}}.
\end{equation}
Therefore, Eq.~\eqref{prob_Dicke_formula} holds.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{bell_thm}]
From the main text, the remainder of the proof of Eq.~\eqref{measurment_basis} involves showing that the LOC consisting of $I_2 \otimes F_m$, $2m$ detectors, and ancilla photons $\ket*{\beta^-}^{\otimes m-1}$ corresponds to the measurement on $\ket*{\beta^\pm}$ expressed as
\begin{equation}\label{measurment_basis2}
    \{\sqrt{1-m^{-1}}\bra{\beta^+}, K'_0, \cdots, K'_m \},
\end{equation}
where
\begin{equation}\label{Kraus'}
    K'_k = \sqrt{\binom{m}{k}}\frac{1}{2^{m/2}}\qty[\bra{\beta^-} + \frac{(m-2k)}{m}\bra{\beta^+}].
\end{equation}
Here, the LOC $I_2 \otimes F_m = F_m \oplus F_m$ preserves the number of photons incident on each block of the LOC.
Thus, we can distinguish the input state according to these numbers.
For simplicity, we introduce the notation $\ket*{\overline{0}}=\ket*{20}$ and $\ket*{\overline{1}}=-\ket*{02}$.
Let $\overline{P}_i(\theta)$ be the operator acting as $\overline{P}_i(\theta)\ket*{x}_i=e^{ix\theta} \ket*{x}_i$ for $x=\{\overline{0},\overline{1}\}$.
Then, we can expand the state as
\begin{equation}
    \ket*{\beta^+}\ket*{\beta^-}^{\otimes m-1} = \frac{1}{2^{m/2}}\overline{P}_0(\pi)(\ket*{\overline{0}}+\ket*{\overline{1}})^{\otimes m}
    =\frac{1}{2^{m/2}} \sum_{k=0}^{m}\binom{m}{k}^{\frac{1}{2}} \overline{P}_0(\pi) \ket*{\overline{D}_k^m},
    \label{expansion}
\end{equation}
where $\ket*{\overline{D}_k^m}$ is the Dicke state in Def.~\ref{Dicke} defined for $\{\ket*{\overline{0}},\ket*{\overline{1}}\}^{\otimes m}$.
The number $k$ in Eq.~\eqref{expansion} is obtained from the measurement result.
The LOC $I \otimes X_m$ corresponds to the permutation of $(0\cdots m-1)$ of the modes of $\{\ket*{\overline{0}},\ket*{\overline{1}}\}^{\otimes m}$.
Thus, Lemma~\ref{probability} yields
\begin{equation}\label{probability_cal}
    \norm*{\rho^{G}_1 {\overline{P}}_0(\pi) \ket*{\overline{D}^m_{k}}}^2 = \norm*{\rho^{C_n}_1 {P}_0(\pi) \ket*{D^m_{k}}}^2 = \qty\Big(\frac{m-2k}{m})^2,
\end{equation}
where $G=\langle I \otimes X_n\rangle$ and $C_m = \langle (0 \cdots m-1)\rangle$.
Because $\ket*{\beta^-}^{\otimes m} \in V_1^G$, we can distinguish $\ket*{\beta^+}$ from $\ket*{\beta^-}$ when we measure $I \otimes X_n \neq 1$, the total probability of which is
\begin{equation}
    \frac{1}{2^m}\sum_{k=0}^{m}\binom{m}{k}\qty\Big(1-\qty\Big(\frac{m-2k}{m})^2)=1-\frac{1}{m}.
\end{equation}
In contrast, when we obtain $I \otimes X_n = 1$, we can expect the corresponding Kraus operators to be Eq.~\eqref{Kraus'}, from Eq.~\eqref{probability_cal}.
To completely identify the measurement basis, we need to perform a kind of measurement tomography, in which we consider multiple bases as the input states and reconstruct the measurement corresponding to this setup.
For input states
\begin{equation}
    \frac{1}{\sqrt{2}}(\ket*{\beta^-}_0\pm i \ket*{\beta^+}_0) = \frac{e^{\pm i \pi/4}}{\sqrt{2}} P_0(\mp \frac{\pi}{2})_0 (\ket*{\overline{0}}+\ket*{\overline{1}}),
\end{equation}
We can calculate the measurement probabilities from Lemma~\ref{probability} in the same manner.
Therefore, we can ensure that the considered setup corresponds to the measurement represented by Eq.~\eqref{measurment_basis2}.
Thus, the entire setup corresponds to the quantum measurement on two qubits represented by Eq.~\eqref{measurment_basis}.
\end{proof}
Following Ref.~\cite{Kim2022-oj}, we introduce the relative entropy of entanglement of quantum measurements as an entanglement measure.
Computing this measure is generally difficult; however, we can explicitly calculate it in the special case as follows.
\begin{definition}\label{measure_def}
    For the $d$-dimensional systems $A$ and $B$, consider a quantum measurement on the composite system $AB$, which is described by positive operator-valued measure (POVM)  elements as $M_{AB}=\{M_x\}$.
    The relative entropy of entanglement of a quantum measurement $M_{AB}$ is
    \begin{equation}
        E(M_{AB})=\frac{1}{d^2} \min_{F_{AB} \in \text{sepM(A:B)}}\sum_x S(M_x||F_x),
    \end{equation}
    where $S( \cdot || \cdot )$ is the quantum relative entropy and $\text{sepM(A:B)}$ is the set of separable measurements on $A$ and $B$~\cite{Watrous2018-la}.
\end{definition}
\begin{lemma}\label{RE_lemma}
    Assume that the rank of all POVM elements is $1$ in Def.~\ref{measure_def}. Then,
    \begin{equation}\label{RE_formula}
        E(M_{AB}) = \frac{1}{d^2} \sum tr(M_x) E_s \qty\Big(\frac{M_x}{\tr(M_x)}),
    \end{equation}
    where $E_s(\cdot)$ is the entanglement entropy of quantum states.
\end{lemma}
\begin{proof}
    Let $M_x =p_x \rho_x$ with $p_x=\tr(M_x)$, and $F_x=q_x \sigma_x$ with $q_x=\tr(F_x)$.
    Then,
    \begin{align}\label{RE_calculation}
        E(M_{AB})
        &= \frac{1}{d^2} \min_{F_{AB} \in \text{sepM(A:B)}}\sum_x \qty\Big{p_x S(\rho_x||\sigma_x) + p_x \log{\frac{p_x}{q_x}}} \notag \\
        &=\frac{1}{d^2} \min_{F_{AB} \in \text{sepM(A:B)}} \qty\Big{\sum_x p_x S(\rho_x||\sigma_x) + S(p||q)},
    \end{align}
    where $S(p||q)$ is the classical relative entropy of $p=\{p_x\}$ with respect to $q=\{q_x\}$.
    $S(p||q)$ takes the minimum value of $0$ when $q=p$.
Here, $\rho_x$ is a pure state, and $S(\rho_x||\sigma_x)$ takes the minimum value of $E_s(\rho_x)$~\cite{Vedral1998-fg}.
    Letting $\mathcal{I}$ and $\Delta$ be the identity and completely dephasing channel for any fixed basis, respectively, the state at the minimum is expressed as $\sigma_x=(\mathcal{I} \otimes \Delta)(\rho_x)$.
    Letting $F_{AB}=\{F_x\}=\{p_i (\mathcal{I} \otimes \Delta)(\rho_x)\}_x$, $F_{AB}$ satisfies the conditions for being a quantum measurement, especially as
    \begin{equation}
        \sum_x F_x = \sum_x p_x (\mathcal{I} \otimes \Delta)(\rho_x) = (\mathcal{I} \otimes \Delta)(I) = I,
    \end{equation}
    from $\sum_x p_x \rho_x = \sum_x M_x = I$.
    Thus, the quantum measurement $F_{AB}$ simultaneously minimizes each term in Eq.~\eqref{RE_calculation} and, hence, Eq.~\eqref{RE_formula} holds.
\end{proof}

We use the following Lemma~\ref{known_thm} without proof.
\begin{lemma}\label{known_thm}
    For the binary entropy function $H_b$,
    \begin{equation}\label{inquality1}
        4p(1-p) \leq H_b(p) \leq (4p(1-p))^{1/ln4}
    \end{equation}
    for $0 \leq p \leq 1$~\cite{Topsoe2001-wd}, and
    \begin{equation}\label{inquality2}
        \sqrt{\frac{n}{8 k (n-k)}}  < \binom{n}{k}2^{-nH_b(k/n)} \leq  \sqrt{\frac{n}{\pi k (n-k)}}
    \end{equation}
    for integers $n$ and $1 \leq  k \leq n-1$~\cite{Cover2006-eh}.
\end{lemma}
\begin{lemma}\label{F_function_lemma}
    Let
    \begin{equation}
        F(x)=-x^2\log_2{\frac{x^2}{x^2+(1-x)^2}}-(1-x)^2\log_2{\frac{(1-x)^2}{x^2+(1-x)^2}}
    \end{equation}
    and
    \begin{equation}\label{f_function_def}
        f(x)=\frac{4x^2(1-x)^2}{x^2+(1-x)^2}.
    \end{equation}
    Then,
    \begin{equation}\label{inequality2}
        f(x) \leq F(x) \leq \frac{1}{2}
    \end{equation}
    for $0\leq x \leq 1$ and $f'(x) \geq 0$ for $0\leq x \leq \frac{1}{2}$~\footnote{We could not find the proof of $F'(x) \geq 0$. Therefore, we introduced the function $f$ and proved $f'(x) \geq 0$ instead.}.
\end{lemma}
\begin{proof}
With the binary entropy function $H_b(x)$, $F$ can be expressed as
\begin{equation}
    F(x)=(x^2+(1-x)^2)H_b\qty\Big(\frac{x^2}{x^2+(1-x)^2}).
\end{equation}
For $0\leq x \leq 1$, it holds that  $(4x(1-x))^{1/ln4}\leq 2\sqrt{x(1-x)}$.
Thus, Eq.~\eqref{inquality1} in Lemma~\ref{known_thm} yields the following inequality:
\begin{equation}\label{F_ub_formula}
    f(x) \leq F(x) \leq 2x(1-x) \leq \frac{1}{2}.
\end{equation}
The first derivative of $f$ is
\begin{equation}
    f'(x)=\frac{8x(1-x)((1-x)-x)(1-x(1-x))}{(x^2+(1-x)^2)^2} \geq 0
\end{equation}
for $0\leq x \leq \frac{1}{2}$.
\end{proof}
\begin{proof}[Proof of Corollary~\ref{bell_values}]
The average success probability $P_m$  of Bell-state discrimination is calculated as the probability that the maximally mixed input state is projected to one of the Bell states.
For even $m$, the case of $k=\frac{m}{2}$ is special because the corresponding Kraus operator becomes the projection onto a Bell state $\ket*{\phi^-}$, with probability $\binom{m}{m/2}2^{-m}$.
Therefore, $P_m$ is calculated as in Eq.~\eqref{success_probability} in total.
The relative entropy of entanglement $E_m$ can be calculated as in Eq.~\eqref{entanglement_measurement} using Lemma~\ref{RE_lemma}.

From Eq.~\eqref{inquality1} in Lemma~\ref{known_thm}, $P_m$ are upper bounded as
\begin{equation}
    P_m \leq P_u(m) = \frac{1}{4}\qty(3-\frac{1}{m}+\sqrt{\frac{4}{\pi m}}).
\end{equation}
From the calculation, $P_u'(x)<0$ for $x>\pi$.
Numerically, we confirmed that $P_u(44) \simeq 0.7868$ and that the maximum value of $P_m$ for $m<44$ is
\begin{equation}
    P_8 = \frac{403}{512} \simeq 0.7871.
\end{equation}
Thus, $P_8$ is the maximum value of $P_m$, from $P_8>P_u(44)$.
Moreover, $P_m \rightarrow \frac{3}{4}$ for $m \rightarrow \infty$ are proved by lower-bounding $P_m$ in the same manner.

Letting
\begin{equation}
    S(m) = \frac{1}{2^{m-1}}\sum_{k=0}^m \binom{m}{k}F\qty\Big(\frac{k}{m}),
\end{equation}
$S(m)\rightarrow 1$ yields $E_m \rightarrow 1$ for $m\rightarrow \infty$.
Because $F(x)\leq \frac{1}{2}$ for $0\leq x \leq 1$, from Lemma~\ref{F_function_lemma}, it holds that
\begin{equation}
    S(m) \leq \frac{1}{2^m}\sum_{k=0}^m \binom{m}{k} = 1.
\end{equation}
Introducing a real number $0< c < \frac{1}{2}$, we obtain the lower bound of $S_m$ as
\begin{equation}
    S _m > \frac{1}{2^{m-1}} \sum_{k=\lceil cm \rceil}^{m-\lceil cm \rceil} \binom{m}{k}F\qty\Big(\frac{k}{m}) \\
    > f(c) \frac{1}{2^{m-1}} \sum_{k=\lceil cm \rceil}^{m-\lceil cm \rceil} \binom{m}{k},
\end{equation}
where $f$ is defined in Eq.~\eqref{f_function_def}, and the second inequality holds from $F(x) \geq f(x)$ and $f'(x) \geq 0$ for $0\leq x< \frac{1}{2}$ in Lemma~\ref{F_function_lemma}.
Further, Eq.~\eqref{inquality1} in Lemma~\ref{known_thm} implies the inequality
\begin{equation}
    \binom{m}{k} \leq 2^{m H_b(\frac{k}{m})}
\end{equation}
for $0\leq k \leq n$.
Thus,
\begin{equation}
    \sum_{k=0}^{\lceil cm \rceil-1} \binom{m}{k}  \leq c m 2^{m H_b(c)}.
\end{equation}
Finally, we have the following lower bound:
\begin{equation}\label{lower_bound}
    S(m) > 2f(c)(1-2cm 2^{-m(1-H_b(c))}).
\end{equation}
It holds that  $H_b(c)<1$ for $c\neq \frac{1}{2}$ and $f(\frac{1}{2})=\frac{1}{2}$.
Thus, Eq.~\eqref{lower_bound} implies that $S(m)\rightarrow 1$ for $m\rightarrow \infty$ and $c\rightarrow \frac{1}{2}$.
\end{proof}


\section{Introduction}\label{sec1}

Many authors submitting \( \sin \cos \tan \inf_{x} \) to NJD journals use \LaTeXe\ to
prepare their papers. This paper describes the
\textsf{WileyNJD-v2.cls} class file which can be used to convert
articles produced with other \LaTeXe\ class files into the correct
form for publication in \emph{Wiley NJD Journals}.

The \textsf{WileyNJD-v2.cls} class file preserves much of the standard
\LaTeXe\ interface so that any document which was produced using
the standard \LaTeXe\ \textsf{article} style can easily be
converted to work with the \textsf{WileyNJD-v2} style. However, the
width of text and typesize will vary from that of
\textsf{article.cls}; therefore, \emph{line breaks will change}
and it is likely that displayed mathematics and tabular material
will need re-setting.

In the following sections we describe how to lay out your code to
use \textsf{WileyNJD-v2.cls} to reproduce the typographical look of
\emph{Wiley NJD Journals}.

\subsection{Procedure to install fonts}

\begin{enumerate}
\item All font files are available under the Lato-Fonts folder
\item Font installer is available under the same folder Windows-Lato-fontinstaller.exe
\item Execute (double click the EXE file) the EXE file that will install all fonts/map files to your local drive.
\end{enumerate}

\subsection{The Three Golden Rules}

Before we proceed, we would like to
stress \emph{three golden rules} that need to be followed to
enable the most efficient use of your code at the typesetting
stage:
\begin{enumerate}
\item[(i)] keep your own macros to an absolute minimum;

\item[(ii)] as \TeX\ is designed to make sensible spacing
decisions by itself, do \emph{not} use explicit horizontal or
vertical spacing commands, except in a few accepted (mostly
mathematical) situations, such as \verb"\," before a
differential~d, or \verb"\quad" to separate an equation from its
qualifier;

\item[(iii)] follow the \emph{NJD} reference style.
\begin{enumerate}[a.]
\item Chemistry --- Use the ``APA'' option as \verb"\documentclass[APA]{WileyNJD-v2.cls}"

\end{enumerate}
\end{enumerate}

\section{Getting Started} The \textsf{WileyNJD-v2.cls} class file should run
on any standard \LaTeXe\ installation. If any of the fonts, class
files or packages it requires are missing from your installation,
they can be found on the \emph{\TeX\ Live} CD-ROMs or from CTAN.

LaTeX document class options

\begin{flushleft}
\begin{enumerate}[a.]
\item LATO1COL--- For LATO font large one column layout use the ``LATO1COL'' option as \verb"\documentclass[APA,LATO1COL]{WileyNJD-v2}"
\item LATO2COL--- For LATO font large two column layout use the ``LATO2COL'' option as \verb"\documentclass[APA,LATO2COL]{WileyNJD-v2}"
\item LATOSMALL--- For LATO font small layout use the ``LATOSMALL'' option as \verb"\documentclass[APA,LATOSMALL]{WileyNJD-v2}"
\end{enumerate}
\end{flushleft}

\section{The Article Header Information}
The heading for any file using \textsf{WileyNJD-v2.cls} is shown in
Figure~\ref{F1}.

\subsection{Remarks}

\begin{enumerate}[(I).]

\item Use \verb"\title{<title> \protect\thanks{<title footnotes}}" for article title and title footnote.
\item Use \verb"\authormark{}" for running heads.

\item Note the use of \verb"\author[<link>]{<name>}" and \verb"\address[<link>]{<name>}" to
link names and addresses. The author for correspondence is marked
by ``*'' and \verb"\corres{}" is used to give that
author's address, which will be printed besides abstract, prefaced by
`Correspondence to:'.

\item For submitting a double-spaced manuscript, add
\verb"doublespace" as an option to the documentclass line. \verb"\documentclass[doublespace]{WileyNJD-v2}"

\item Use \verb"\presentaddress{}" for present address.

\item In abstract \verb"\abstract[<title>]{abstract paragraph}" use optional parameter for title followed by abstract paragraph.

\item For Key words use \verb"\keywords{}".

\item For how to site use \verb"\jnlcitation{\cname{\author{<author name>},"

\verb"\ctitle{<title>}, \cjournal{<Journal name>}, \cvol{<vol>}.}}".

\item For title page abbreviations use \verb"\footnotetext{<\textbf{Abbreviation title:} Abreviations>}"

\item Use \verb"\articletype{<article category}" for article header information

\item Use \verb"\received{<received date>} \revised{<revised date>} \accepted{<accepted date>}" for history dates.

\end{enumerate}

\begin{figure}[p]
\setlength{\fboxsep}{0pt
\setlength{\fboxrule}{0pt
\begin{center}
\begin{boxedverbatim}

\documentclass[APA,LATO1COL]{WileyNJD-v2}

\articletype{Article Type

\received{26 April 2016}
\revised{6 June 2016}
\accepted{6 June 2016}

\begin{document}

\title{<Initial cap, lower case>\protect\thanks{<title footnote.>}}

\author[<address link>]{<Author name><corresponding author*>}

\author[<address link>,<address link>]{Author Name}

\authormark{AUTHOR ONE \textsc{et al}}

\address[<address link>]{\orgdiv{<Org Division>}, \orgname{<Org name>},
\orgaddress{\state{<State name>}, \country{<Country name>}}}
\address[<address link>]{\orgdiv{<Org Division>}, \orgname{<Org name>},
 \orgaddress{\state{<State name>}, \country{<Country name>}}}

\corres{<corresponding author link*> <author name, address.
 \email{<authorone@email.com>}}

\presentaddress{<Present address>}

\abstract[<Abstract heading>]{<Abstract paragraph>}

\keywords{<keyword1>, <keyword2>,...}

\jnlcitation{\cname
\author{<aurhor name>},
\author{<aurhor name>},
\author{<aurhor name>},
\author{<aurhor name>},  and
\author{<aurhor name>}} (\cyear{<year>}),
\ctitle{<journal title>}, \cjournal{<journal name>} <year> <vol> Page <xxx>-<xxx>}

\footnotetext{\textbf{<abbreviation head:>} <abbreviations> ..}

\maketitle

\section{Introduction}
.
.
.
\end{boxedverbatim}
\end{center}
\caption{Example for title page.\label{F1}}
\end{figure}

\section{The Body of the Article}

\subsection{Section headings}

\begin{enumerate}[(H1)]
\item Section --- use \verb"\section{}"
\item SubSection--- use \verb"\subsection{}"
\item SubSubSectioin--- use \verb"\subsubsection{}"
\item Paragraph--- use \verb"\paragraph{}"
\item Subparagraph--- use \verb"\subparagraph{}"
\end{enumerate}

\subsection{Mathematics} \textsf{WileyNJD-v2.cls} makes the full
functionality of \AmS\/\TeX\ available. We encourage the use of
the \verb"align", \verb"gather" and \verb"multline" environments
for displayed mathematics.

\subsection{Figures and Tables}

\textsf{WileyNJD-v2.cls} uses the
\textsf{graphicx} package for handling figures.

Figures are called in as follows:
\begin{verbatim}
\begin{figure}
\centering
\includegraphics{<figure name>}
\caption{<Figure caption>}
\end{figure}
\end{verbatim}

The standard coding for a table is shown in Figure~\ref{F2}.

\begin{figure}[h]
\setlength{\fboxsep}{0pt
\setlength{\fboxrule}{0pt
\begin{center}
\begin{boxedverbatim}
\begin{table}
\caption{<Table caption>}
\centering
\begin{tabular}{<table alignment>}
\toprule
<column headings>\\
\midrule
<table entries
(separated by & as usual)>\\
<table entries>\\
.
.
.\\
\bottomrule
\end{tabular}
\begin{tablenotes}
\item Source: xxx.
\item[1] xxx.
\item[2] xxx.
\end{tablenotes}
\end{table}
\end{boxedverbatim}
\end{center}
\caption{Example for table layout.\label{F2}}
\end{figure}

\subsection{Cross-referencing}
The use of the \LaTeX\ cross-reference system
for figures, tables, equations, etc., is encouraged
(using \verb"\ref{<name>}" and \verb"\label{<name>}").

\subsection{Box text}

\begin{verbatim}
\begin{boxtext}
\section*{<title>
Paragraph
\end{boxtext}
\end{verbatim}

\subsection{List items}

\subsubsection{Enumerate list styles}
\begin{verbatim}
\begin{enumerate}[1]
\item
\end{enumerate}

\begin{enumerate}[1.]
\item
\end{enumerate}

\begin{enumerate}[(1)]
\item
\end{enumerate}

\begin{enumerate}[I]
\item
\end{enumerate}

\begin{enumerate}[i]
\item
\end{enumerate}

\begin{enumerate}[a]
\item
\end{enumerate}

\end{verbatim}

\subsubsection{Bullet list styles}

\begin{verbatim}
\begin{itemize}
\item
\end{itemize}
\end{verbatim}

\subsubsection{Description list}

\begin{verbatim}
\begin{description}
\item[<entry>] description text.
\end{description}
\end{verbatim}

\subsection{Enunciations}

\begin{verbatim}
\begin{theorem}[<Theorem subhead>]\label{thm1}
<theorem text>.
\end{theorem}

\begin{proposition}[<proposition subhead>]\label{pro1}
<proposition text>.
\end{proposition}

\begin{definition}[<definition subhead>]\label{dfn1}
<definition text>.
\end{definition}

\begin{proof}
<proof text>.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm1}]
<proof text>.
\end{proof}

\end{verbatim}

\subsection{Program codes}

Use \verb+\begin{verbatim}...\end{verbatim}+ for program codes without math. Use \verb+\begin{alltt}...\end{alltt}+ for program codes with math. Based on the text provided inside the optional argument of \verb+\begin{code}[Psecode|Listing|Box|Code|+\hfill\break \verb+Specification|Procedure|Sourcecode|Program]...+ \verb+\end{code}+ tag corresponding boxed like floats are generated. Also note that \verb+\begin{code}[Code|Listing]...+ \verb+\end{code}+ tag with either Code or Listing text as optional argument text are set with computer modern typewriter font.  All other code environments are set with normal text font. Refer below example:

\begin{verbatim}
\begin{lstlisting}[caption={Descriptive Caption Text},label=DescriptiveLabel]
for i:=maxint to 0 do
begin
{ do nothing }
end;
Write('Case insensitive ');
WritE('Pascal keywords.');
\end{lstlisting}
\end{verbatim}

\subsection{Acknowledgements} An Acknowledgements section is started with \verb"\ack" or
\verb"\acks" for \textit{Acknowledgement} or
\textit{Acknowledgements}, respectively. It must be placed just
before the References.

\subsection{Bibliography}

\begin{enumerate}[1]
\item Use \verb"
\section{Introduction}\label{introduction}}

Smooth-in-time accelerated failure time (AFT) models can estimate
absolute risks by modeling the hazard directly through a user-specified
baseline hazard distribution \citep{kleinbaum2012survival}. Cox
proportional hazards models are used more often than AFT models, causing
analyses to be based on hazard ratios and relative risks rather than on
survival curves and absolute risks \citep{hanley2009}. The
identification of an appropriate distribution for the baseline hazard in
an AFT model may be difficult for common diseases that have many
interacting risk factors, or a Cox model where the disease pathogenesis
may change with age, making it difficult to maintain the proportional
hazards assumption. For example, previous studies of breast cancer
incidence have discovered time-varying interactions with covariates of
interest, such as tumor size \citep{coradini2000time}. One approach to
provide flexibility in the baseline hazard involves using the basis of
splines on time in our model \citep{royston2002flexible}. However,
regression-based models are limited in that they require prior knowledge
about potential time-varying interactions and their quantitative
effects.

Neural networks provide a data-driven approach to approximating
interaction terms. For example, DeepSurv is a neural network-based
proportional hazards model that implements the Cox partial
log-likelihood as a custom loss function \citep{katzman2018DeepSurv},
resulting in a stepwise absolute risk curve that cannot accommodate
time-varying interactions. Compared to Cox regression, DeepSurv shows
better performance on the Study to Understand Prognoses Preferences and
Risk Treatments (SUPPORT) dataset \citep{knaus1995SUPPORT}. To handle
non-proportional hazards, a modification to the loss function was
proposed \citep{faraggi1995neural}. To remove the need for tuning the
number of layers and nodes, the concept of extreme learning machines has
been applied to a Cox neural network model \citep{wang2019extreme}.

As an alternative method that assumes a baseline hazard distribution,
DeepHit specifies each survival time of interest in the model and
directly estimates survival curves, rather than deriving a hazard
function \citep{lee2018DeepHit}. It assumes an inverse Gaussian
distribution as the baseline hazard and it outperforms DeepSurv on the
Molecular Taxonomy of Breast Cancer International Consortium (METABRIC)
dataset \citep{curtis2012genomic}.

Providing further flexibility, Deep Survival Machines (DSM) is a
parametric survival model using neural networks with a mixture of
distributions as the baseline hazard \citep{dsmPaper}. On both the
SUPPORT and METABRIC datasets, DSM outperforms DeepSurv and DeepHit
\citep{dsmPaper}. However, like DeepHit and DeepSurv, DSM cannot model
time-varying interactions.

We note that these alternative neural network approaches all require
custom loss functions \citep{katzman2018DeepSurv} \citep{lee2018DeepHit}
\citep{dsmPaper}. DeepHit introduces a hyperparameter weighing its two
loss functions (negative log-likelihood and ranking losses) while DSM
requires a two-phase learning process and user implementations for
distributions beyond Log-Normal or Weibull \citep{lee2018DeepHit}
\citep{dsmPaper}. Regression-based approaches require prior
specification of all interaction terms, which makes it challenging to
model covariate effects that change over time. The current neural
network models provide flexibility at the cost of opacity, while
regression models provide clarity at the cost of flexibility.

In this article, we propose Case-Base Neural Networks (CBNN) as a method
that models time-varying interactions and a flexible baseline hazard
using commonly available neural network components. Our approach to
modeling the full hazard uses case-base sampling \citep{hanley2009}.
This sampling technique allows probabilistic models to predict survival
outcomes. As part of the case-base framework, we use transformations of
time as a feature (covariate) to specify different baseline hazards. For
example, by including splines of time as covariates, we can approximate
the Royston-Parmar flexible baseline hazard model
\citep{royston2002flexible}, however, this still requires explicit use
of time-varying interactions. CBNN can model both without extra tuning
parameters.

In Section \ref{methods}, we describe how case-base sampling and neural
networks are combined both conceptually and algebraically, along with
our hyperparameter choices and software implementation. Section
\ref{sims}, describes our metrics and compares the performance of CBNN,
DeepSurv, DeepHit, DSM, Cox regression and case-base using logistic
regression (CBLR) on simulated data. Section \ref{casestudies} describes
the real-data analysis, while Section \ref{discussion} explores the
implications of our results and contextualizes them within neural
network survival analysis in a single event setting.

\hypertarget{methods}{%
\section{Case-base neural network methodology, metrics and
software}\label{methods}}

In this section, we define case-base sampling, which converts the total
survival time into discrete person-specific moments (person-moments).
Then, we detail how neural networks can be used within this framework,
explicitly incorporating time as a feature while adjusting for the
sampling bias. Finally, we report on the software versions used. An R
package is available for use at
\url{https://github.com/Jesse-Islam/cbnn}. The entire code base to
reproduce the figures and empirical results in this paper is available
at \url{https://github.com/Jesse-Islam/cbnnManuscript}.

\hypertarget{case-base-sampling}{%
\subsection{Case-base sampling}\label{case-base-sampling}}

Case-base sampling is an alternative framework for survival analysis
\citep{hanley2009}. In case-base sampling, we sample from the continuous
survival time of each person in our dataset to create a \emph{base
series} of \emph{person-moments}. This \emph{base series} complements
the \emph{case series}, which contains all person-moments at which the
event of interest occurs.

For each person-moment sampled, let \(X_i\) be the corresponding
covariate profile \(\left(x_{i1},x_{i2},...,x_{ip} \right)\), \(T_i\) be
the time of the person-moment and \(Y_i\) be the indicator variable for
whether the event of interest occurred at time \(T_i\). We estimate the
hazard function \(h(t \mid X_i)\) using the sampled person-moments.
Recall that \(h(t \mid X_i)\) is the instantaneous potential of
experiencing the event at time \(t\) for a given set of covariates
\(X_i\), assuming \(T_i \geq t\).

Now, let \(b\) be the (user-defined) size of the \emph{base series} and
let \(B\) be the sum of all follow-up times for the individuals in the
study. If we sample the \emph{base series} uniformly across the study
base, then the hazard function of the sampling process is equal to
\(b/B\). Therefore, we have the following equality
\footnote{We are abusing notation here, conflating hazards with probabilities. For a rigorous treatment, see Saarela \& Hanley (2015) section 3 \cite{saarela2015}.}:
\begin{align}\label{eqn:main}
\frac{P\left(Y_i=1 \mid X_i, T_i\right)}{P\left(Y_i = 0 \mid X_i, T_i\right)} = \frac{h\left(T_i \mid X_i\right)}{b/B}.
\end{align}

To provide some intuition for this equation, the hazard is the rate at
which the event (\(Y_{i}=1\)) occurs at a given time, conditional on
your survival up to the follow-up-time of interest. As we are looking at
person-moments, we may retrieve the rate at which an event does not
occur (\(Y_{i}=0)\) at a given time by taking the number of samples in
\(b\) and dividing it by the total number of moments in the study base
\(B\) at the current follow-up-time. by taking the ratio. of these two metics,
we end up with our odds on the left hand side. The odds of a person-moment
being a part of the \emph{case series} is the ratio of the hazard
\(h(T_i \mid X_i)\) and the uniform rate \(b/B\). Though an intuitive approach
is used to describe this equation, a complete derivation of this equation
can be found at \citep{saarela2015}.

Sampling techniques are known to result in a loss of information. To
address this concern, the original case-base methodology paper uses a
measure for relative information when comparing two averages, based on
the size of our base-series and case-series. If we let the size of
\(b=100c\), then the variances and covariances are expected to be
proportional to \(\frac{1}{c}+\frac{1}{100c}\) rather than
\(\frac{1}{c}+\frac{1}{\inf}\), about 1 percent larger. If there is a
concern with a 1 percent inflation of these estimates, a larger ratio
than 100 may be used. As there are concerns about dependency due
multiple person-moments coming from each individual, the information
loss was studied in Saarela 2015 (\citep{saarela2015}), where they found
a ratio of 100:1 base-series to case series gives an efficiency
comparable to Mantel-Haenszel and conditional logistic regression
methods \citep{saarela2015}.

\begin{align}\label{eqn:offset}
\log \left( h\left(t \mid X_i\right)\right) = \log \left(\frac{P\left(Y_i = 1 \mid X_i, t\right)}{P\left(Y_i = 0 \mid X_i, t\right)}\right) + \log\left(\frac{b}{B}\right).
\end{align}

To estimate the correct hazard function, we adjusting for the bias
introduced when sampling a fraction of the study base \(B\) by adding an
offset term \(\log\left(\frac{B}{b} \right)\) \eqref{eqn:offset}. Next,
we propose using neural networks to model the odds.

\hypertarget{neural-networks-to-model-the-hazard-function}{%
\subsection{Neural networks to model the hazard
function}\label{neural-networks-to-model-the-hazard-function}}

After case-base sampling, we pass all features, including time, into any
user-defined feed-forward neural network, to which an offset term is added,
then passed through a sigmoid activation function (Figure
\ref{fig:NNarch}). We use the sigmoid activation function as it is the
inverse of the log-odds and commonly available in popular neural network
packages (keras, pytorch, etc) \citep{keras} \citep{pytorch}. which we can use
to calculate the hazard. The general form for the neural network using
CBNN is:

\begin{align}\label{eqn:nnProb}
P\left(Y=1|X,T\right)&=\mathrm{sigmoid}\left(f_{\theta}(X, T) + \log\left(\frac{B}{b}\right) \right),
\end{align}

where \(T\) is a random variable representing the event time, \(X\) is
the random variable for a covariate profile, \(f_{\theta}(X, T)\)
represents any feed-forward neural network architecture,
\(\log\left(\frac{B}{b}\right)\) is bias term set by case-base sampling,
\(\theta\) is the set of parameters learned by the neural network and
\(\mathrm{sigmoid}(x)=\frac{1}{1+e^{-x}}\). By approximating a
higher-order polynomial of time using a neural network, the baseline
hazard specification is now data-driven, where user-defined
hyperparameters such as regularization, number of layers and nodes
control the flexibility of the hazard function. We provide a detailed
description of the choices we made in the next sub-section.

\begin{figure}

{\centering \includegraphics[width=1\linewidth]{figures/nnarch2.PNG}

}

\caption{Steps involved in CBNN from case-base sampling to the model framework we use for training. The first step is case-base sampling, completed before training begins. Next, we pass this sampled data through a feed-forward neural network. We add the offset and pass that through a sigmoid activation function, whose output is a probability. Once the neural network model completes its training, we can convert the probability output to a hazard, using it for our survival outcomes of interest.}\label{fig:NNarch}
\end{figure}

The following derivation shows how our probability estimate is converted
to odds: \begin{align*}
 \log\left( h(t \mid X) \right) &= \log\left(\frac{\mathrm{sigmoid}\left(f_{\theta}(X, T) + \log\left(\frac{B}{b}\right)\right)}{1-\mathrm{sigmoid}\left(f_{\theta}(X, T) + \log\left(\frac{B}{b}\right)\right)}\right) + \log\left(\frac{b}{B}\right) \\
 &= \log\left( \frac{\frac{\exp\left(f_{\theta}(X, T) + \log\left(\frac{B}{b}\right)\right)}{\exp\left(f_{\theta}(X, T) + \log\left(\frac{B}{b}\right)\right)+1}}{1-\frac{\exp\left(f_{\theta}(X, T) + \log\left(\frac{B}{b}\right)\right)}{\exp\left(f_{\theta}(X, T) + \log\left(\frac{B}{b}\right)\right)+1}}\right) + \log\left(\frac{b}{B}\right) \\
 &= \log\left(\exp\left( f_{\theta}(X, T) + \log\left(\frac{B}{b}\right) \right) \right) + \log\left(\frac{b}{B}\right) \\
 &= f_{\theta}(X, T) + \log\left(\frac{B}{b}\right) + \log\left(\frac{b}{B}\right) \\
&= f_{\theta}(X, T).
\end{align*}

We use binary cross-entropy as our loss function \citep{gulli2017}:
\begin{align*}
L(\theta)=-\frac{1}{N} \sum^{N}_{i=1} y_{i} \cdot \log(\hat{f}_{\theta}(x_{i}, t_{i}) ) + (1-y_{i} )\cdot \log(1-\hat{f}_{\theta}(x_{i}, t_{i}) ),
\end{align*} where \(\hat{f}_{\theta}(x_{i}, t_{i})\) is our estimate
for a given covariate profile and time, \(y_{i}\) is our target value
specifying whether an event occurred and \(N\) represents the number of
individuals in our training set.

Backpropagation with an appropriate minimization algorithm (e.g.~Adam,
RMSPropagation, stochastic gradient descent) is used to optimize the
parameters in the model \citep{gulli2017}. For our analysis, we use Adam
as implemented in Keras \citep{gulli2017}. Note that the size of the
\emph{case series} is fixed as the number of events, but we can make the
\emph{base series} as large as we want. A ratio of 100:1 \emph{base
series} to \emph{case series} is sufficient \citep{hanley2009}. We pass
our feed-forward neural network through a sigmoid activation function
(Figure \ref{fig:NNarch}). Finally, we can convert this model output to
a hazard. When using our model for predictions, we manually set the
offset term to 0 in the new data, as we account for the bias during the
fitting process.

Since we are directly modeling the hazard, we can readily estimate the
risk function (\(F\)) at time \(t\) for a covariate profile \(X\),
viz.\begin{align}\label{eqn:ci2}
F\left(t\mid X\right)& = 1 - \exp\left(-\int_{0}^{t}h(u|X) \,\textrm du\right).
\end{align} We use a finite Riemann sum \citep{hughes2020calculus} to
approximate the integral in \eqref{eqn:ci2}.

\hypertarget{hyperparameter-selection}{%
\subsection{Hyperparameter selection}\label{hyperparameter-selection}}

Neural networks are flexible when defining the architecture and
optimization parameters. These hyperparameter decisions can affect the
estimated parameters and were chosen during a set of initial simulations
to determine if CBNN can learn complex interactions in practice. We set
a max number of epochs to be \(2000\), batch size as \(512\), learning
rate as \(10^{-3}\), decay as \(10^{-7}\), patience as \(10\) epochs,
\(\{50,50,25,25\}\) nodes in each hidden layer with \(50\%\) dropout at
each layer, a minimum delta loss on the validation set of \(10^{-7}\)
over 10 epochs and Adam \citep{gulli2017} as our optimizer. These
choices may permit the model to approximate higher-order interactions
while preventing over-fitting \citep{srivastava2014dropout}. We fix a
train-validation-test split that allows us to update the weights with a
subset of the data (training), assess performance at each epoch
(validation) and gauge performance of the final model (test) for each
method. We select the best weights after training based on validation
loss \citep{gulli2017}.

\hypertarget{software-implementation}{%
\subsection{Software implementation}\label{software-implementation}}

R \citep{Rsoft} and python \citep{py} are used to evaluate methods from
both languages. We fit the Cox model using the \textbf{survival} package
\citep{survpkg}, the CBLR model using the \textbf{casebase} package
\citep{cbpkg}, the DeepSurv model using the \textbf{survivalmodels}
package \citep{survmods}, the DeepHit model using \textbf{pycox}
\citep{lee2018DeepHit} and the DSM model using
\textbf{DeepSurvivalMachines} \citep{dsmPaper}. We made the components
of CBNN using the \textbf{casebase} package for the sampling step and
the \textbf{keras} \citep{keras} package for our neural network
architecture. The \textbf{simsurv} package \citep{simsurv} is used for
our simulation studies, while \textbf{flexsurv} \citep{flexsurv} is used
to fit a flexible baseline hazard using splines for our complex
simulation. We use the implementation of \(C_{IPCW}\) from the python
package \textbf{sksurv} \citep{sksurv}. The \textbf{riskRegression}
package \citep{riskRegression} is used to get the Index of Prediction
Accuracy (IPA metric). Both metrics are described in detail in the
following section. We modify the \textbf{riskRegression} package to be
used with any user supplied risk function \(F\). To ensure that both R
and Python-based models are running in unison on the same data through
our simulations and bootstrap, we use the \textbf{reticulate} package
\citep{reticulate}.

\hypertarget{sims}{%
\section{Simulation studies}\label{sims}}

In this section, we use simulated data to evaluate the performance of
CBNN and compare our approach with existing regression-based (Cox, CBLR)
and neural network-based (DeepHit, DeepSurv, DSM) methods. We specify a
linear combination of each covariate as the linear predictor in
regression-based approaches (Cox, CBLR), which contrasts with neural
network approaches that allow for non-linear interactions. We simulate
data under a simple exponential model and a complex baseline hazard with
time-varying interactions, each with 2000 individuals and 10\% random censoring. For both
settings, we simulate three covariates:

\[
z_{1} \sim \textrm{Bernoulli}(0.5) \qquad \qquad
z_{2} \sim \begin{cases}
 N(0,0.5) & \textrm{if } z_{1}=0\\
 N(10,0.5) & \textrm{if } z_{1}=1
\end{cases} \qquad \qquad
z_{3} \sim \begin{cases}
 N(8,0.5) & \textrm{if } z_{1}=0\\
 N(-3,0.5) & \textrm{if } z_{1}=1.
\end{cases}
\]

The DeepHit-specific hyperparameter alpha is set to \(0.5\) (equal
weight between its negative log-likelihood and ranking losses
\citep{lee2018DeepHit}). We modify the \textbf{DeepSurvivalMachines}
\citep{dsmPaper} package to include dropout and a minimum delta loss
during the fitting process. For DSM, we define a mixture model of six
Weibull distributions for the baseline hazard. All other hyperparameters
are held constant across all neural network methods in both the
simulation studies and real data applications.

Besides the methods mentioned above, we include the Optimal model in our
comparisons using CBLR. That is, we include the exact functional form of
the covariates in a CBLR model (referred to as Optimal for simplicity).
We calculate \(t\)-based 95\% confidence intervals using 100
replications of the simulated data. For all analyses, we use 80\% for
training and 20\% for the test set. 20\% of the training set is kept for
validation at each epoch. We predict risk functions \(F\) using
\eqref{eqn:ci2} for individuals in the test set, which are used to
calculate our \(C_{IPCW}\) and IPA scores.

\hypertarget{performance-metrics}{%
\subsection{Performance metrics}\label{performance-metrics}}

We use two metrics to assess the performance of the different methods of
interest: 1) (IPA) \citep{kattan2018index} and 2) inverse probability
censoring weights-adjusted concordance index (\(C_{IPCW}\))
\citep{uno2011}, which we define below. We may be interested in an
absolute risk a specific number of years in the future. For the sake of
comparison, we plot each metric over time \(t\) so that the performance
is transparent. As we go further along follow-up-time, we expect the
performance of all models to drop as there are fewer individuals left in
the study by then. For \(C_{IPCW}\), the score demonstrates performance
up to time \(t\) \citep{uno2011}, while IPA is performance at time
\(t\).

\hypertarget{index-of-prediction-accuracy-ipa}{%
\subsubsection{Index of prediction accuracy
(IPA)}\label{index-of-prediction-accuracy-ipa}}

The IPA is a function of the Brier score (\(BS(t)\)) \citep{graf1999},
which is defined as \begin{align}
BS(t)=\frac{1}{N}\sum^{N}_{i=1}\left(\frac{\left(1 - \widehat{F}(t \mid X_{i})\right)^{2}\cdot I(T_{i}\leq t,\delta_{i}=1)}{\widehat{G}(T_{i})} + \frac{\widehat{F}(t\mid X_{i})^{2}\cdot I(T_{i}>t)}{\widehat{G}(t)}\right),
\end{align} where \(\delta_{i}=1\) shows individuals who have
experienced the event, \(N\) represents the number of samples in our
dataset over which we calculate \(BS(t)\), \(\widehat{G}(t)=P[c>t]\) is
a non-parametric estimate of the censoring distribution, \(c\) is
censoring time and \(T_{i}\) is an individual's survival or censoring
time. The Brier score provides a score that accounts for the information
loss because of censoring. There are three categories of individuals
that may appear within the dataset once we fix our \(t\) of interest.
Individuals who experienced the event before \(t\) are present in the
first term of the equation. The second term of the equation includes
individuals who experience the event or are censored after \(t\). Those
censored before \(t\) are the third category of people. The inverse
probability censoring weights (IPCW) adjustment (\(G(\cdot)\)) is to
account for these category three individuals whose information is
missing. The IPA score as a function of time is given by \begin{align}
\textrm{IPA}(t) &= 1-\frac{BS_{model}(t)}{BS_{null}(t)}, \nonumber
\end{align} where \(BS_{model}(t)\) represents the Brier score over time
\(t\) for the model of interest and \(BS_{null}(t)\) represents the
Brier score if we use an unadjusted Kaplan-Meier (KM) curve as the
prediction for all observations \citep{kattan2018index}. Note that IPA
has an upper bound of one, where positive values show an increase in
performance over the null model and negative values show that the null
model performs better. These scores show how performance changes over
follow-up time.

A potential artifact of IPA is that the score is unstable at earlier and
later survival times. This is because of a near equivalent Brier score
among each model and the null model. At small values, a difference of
\(0.1\) creates a more significant fold change than at larger values. As
the Brier score is potentially small at the start and at the end of
their respective curves, The IPA score may be unstable at the same
locations.

\hypertarget{inverse-probability-censoring-weights-adjusted-concordance-index}{%
\subsubsection{Inverse probability censoring weights-adjusted
concordance
index}\label{inverse-probability-censoring-weights-adjusted-concordance-index}}

The \(C_{IPCW}\) is a non-proper, rank-based metric that does not depend
on the censoring times in the test data \citep{uno2011}. The
\(C_{IPCW}\) is given by

\begin{align} \label{eq:cidx}
C_{IPCW}(t) &= \frac{\sum^{N}_{i=1}\sum^{N}_{j=1}\delta_{i}\left\{\widehat{G}(T_{i})\right\}^{-2} I(T_{i}<T_{j},T_{i}<t) I\left(\widehat{F}(t|X_{i})>\widehat{F}(t|X_{j})\right)}{\sum^{N}_{i=1}\sum^{N}_{j=1}\delta_{i}\left\{\widehat{G}(T_{i})\right\}^{-2} I(T_{i}<T_{j},T_{i}<t)}.
\end{align} where the follow-up period of interest is (0,\(t\)),
\(I(\cdot)\) is an indicator function, \(\widehat{F}(X_{i},t)\) is the
risk function estimated for everyone in the study at time \(t\) and
\(C_{IPCW}\) can compare the performance of different models, where a
higher score is better. Note that the \(C_{IPCW}\) may produce
misleading performance, as it ranks based on survival times, not event
status \citep{cindexfails2019}. This metric is considered an unbiased
population concordance measure because of the IPCW adjustment
\citep{uno2011}.

\begin{figure}

{\centering \includegraphics[width=1\linewidth]{analyses/figures/megaPlot}

}

\caption{Summarizes the simple simulation (A, B), complex simulation (C, D), SUPPORT case study (E, F) and METABRIC case study (G, H) results. The first row shows the IPA score for each model in each study over follow-up time. Negative values mean our model performs worse than the null model and positive values mean the model performs better. The second row demonstrates the $C_{IPCW}$ score for each model in each study over follow-up time. A score of 1 is the maximum performance for either metric. Each model-specific metric in each study shows a 95-percent confidence interval over 100 iterations. The models of interest are case-base with logistic regression (CBLR), Case-Base Neural Networks (CBNN), Cox, DeepHit, DeepSurv, Deep Survival Machines (DSM), Optimal (a CBLR model with the exact interaction terms and baseline hazard specified) and Kaplan-Meier (to serve as a baseline, predicting the average for all individuals). Note than Cox and CBLR are superimposed, due to very similar prediction outcomes.}\label{fig:megaPlot}
\end{figure}

\hypertarget{simple-simulation-constant-baseline-hazard}{%
\subsection{Simple simulation: constant baseline
hazard}\label{simple-simulation-constant-baseline-hazard}}

We simulate data from a simple model that primarily depends on a
constant baseline hazard:
\[\log h(t \mid X_i) = \beta_{{1}}z_{1}+\beta_{{2}}z_{2}+\beta_{{3}}z_{3},\]
The covariate effects are given by
\(\beta_{{1}}=0.1, \beta_{{2}}=0.1, \beta_{{3}}=0.1\). Once we simulate
survival times, we introduce 10\% random censoring.

\hypertarget{performance-comparison-in-simple-simulation}{%
\subsubsection{Performance comparison in simple
simulation}\label{performance-comparison-in-simple-simulation}}

Figure \ref{fig:megaPlot} A, B and Table \ref{tab:megaTable} A show the
results for the simple simulation. The regression-based methods (CBLR,
Cox, Optimal) outperform the neural network ones in the simple
simulation setting. Among the neural network approaches, CBNN
outperforms all other methods in terms of both IPA and \(C_{IPCW}\)
(Figure \ref{fig:megaPlot} A, B). Specifically, we see CBNN is
consistent across time with smaller confidence intervals compared to
DeepHit, DeepSurv and DSM.

In this simple setting, the regression models are much closer to the
Optimal model, while the neural network models perform worse than the KM
null model. The wide confidence bands in \(C_{IPCW}\) suggest the neural
network models may be over-parameterized.

\hypertarget{complex-simulation-flexible-baseline-hazard-time-varying-interactions}{%
\subsection{Complex simulation: flexible baseline hazard, time-varying
interactions}\label{complex-simulation-flexible-baseline-hazard-time-varying-interactions}}

This simulation demonstrates performance with the presence of a complex
baseline hazard and a time-varying interaction. Originally used to show
the spline-based hazard model proposed by Royston and Parmar
\citep{royston2002flexible}, the breast cancer dataset provides a
complex hazard from which we simulate, available in the
\textbf{flexsurv} R package \citep{flexsurv}. To increase the complexity
of our data-generating mechanism for this simulation, we design the
model as follows: \begin{align}
\log h(t \mid X_i) =\sum_{i=1}^{5} (\gamma_{i} \cdot \psi_{i}) + \beta_{{1}} (z_{1}) + \beta_{{2}} (z_{2})+ \beta_{{3}} (z_{3})+ \tau_{1} ( z_{1} \cdot z_{2} \cdot t)+ \tau_{2} ( z_{1} \cdot z_{3})+ \tau_{3} (z_{2} \cdot z_{3}), \nonumber
\end{align}

\FloatBarrier
\begin{table}
\caption{Four tables representing performance at certain follow-up times for the simple simulation, complex simulation, SUPPORT and METABRIC. Each table shows performance for each method in each study at $25\%$, $50\%$. $75\%$ and $100\%$ of follow-up time. The bold elements show the best model for each study, at each follow-up time of interest. These tables are included to provide exact measures at certain intervals. The models of interest are: Cox, case-base with logistic regression (CBLR), DeepSurv, DeepHit, Case-Base Neural Network (CBNN), Optimal and Deep Survival Machines (DSM).}
\label{tab:megaTable}

\begin{center}\includegraphics[width=1\linewidth]{analyses/figures/megaTable} \end{center}

\end{table}

where
\(\gamma_{1}=3.9, \gamma_{2}=3, \gamma_{3}=-0.43, \gamma_{4}=1.33,\gamma_{5}=-0.86, \beta_{{1}}=1, \beta_{{2}}=1, \beta_{{3}}=1, \tau_{1}=10, \tau_{2}=2, \tau_{3}=2\)
and \(\psi\) are basis splines. The \(gamma\) coefficients are obtained
from an intercept-only cubic splines model with three knots using the
\emph{flexsurvspline} function from the \textbf{flexsurv} package
\citep{flexsurv}. Note that we fix these values for the analysis. The
\(\beta\) coefficients represent direct effects, \(\tau_{2}\) and
\({\tau_3}\) represent interactions and \(\tau_{1}\) is a time-varying
interaction.

\hypertarget{performance-comparison-in-complex-simulation}{%
\subsubsection{Performance comparison in complex
simulation}\label{performance-comparison-in-complex-simulation}}

Figure \ref{fig:megaPlot} C, D and Table \ref{tab:megaTable} B show the
performance over time on a test set in the complex simulation. Apart
from the Optimal regression model, CBNN outperforms the competing models
when examining IPA and up to the 75-percentile of follow-up time for
\(C_{IPCW}\). We expected the Optimal model to perform best in both
metrics. However, this was not the case for \(C_{IPCW}\) and may be due
to an artifact of concordance-based metrics, where a misspecified model
may perform better than a correctly specified one
\citep{cindexfails2019}. We attribute the performance of CBNN to its
flexibility in modeling time-varying interactions and baseline hazard,
flexibility the other neural network models do not have.

\hypertarget{casestudies}{%
\section{Application to SUPPORT and METABRIC data}\label{casestudies}}

Our complex simulation demonstrates the superior performance of CBNN in
ideal conditions with clean data. To obtain a more realistic performance
assessment, we compared models using two real datasets with a
time-to-event outcome. The first case study examines the SUPPORT dataset
\citep{knaus1995SUPPORT}. The second case study examines the METABRIC
dataset \citep{curtis2012genomic}. We use the same hyperparameters as in
the simulation studies. As we do not know the true model for the real
data, we exclude the Optimal model. We split the datasets keeping 20\%
of the observations as a test set. 20\% of the training set is kept
aside for validation at each epoch. We predict risk functions for
everyone in the test set, which is used to calculate our metrics. We
conduct 100 bootstrap re-samples for the real data applications to
obtain confidence intervals.

\hypertarget{performance-evaluation-using-the-support-dataset}{%
\subsection{Performance evaluation using the SUPPORT
dataset}\label{performance-evaluation-using-the-support-dataset}}

The SUPPORT dataset tracks the time until death for seriously ill
patients at five American hospitals \citep{knaus1995SUPPORT}. We use a
pre-processed version of the dataset made available in the DeepSurv
package \citep{katzman2018DeepSurv}. This dataset contains 9104 samples
and 14 covariates (age, sex, race, number of comorbidities, presence of
diabetes, presence of dementia, presence of cancer, mean arterial blood
pressure, heart rate, respiration rate, temperature, white blood cell
count, serum's sodium and serum's creatinine)
\citep{katzman2018DeepSurv}. Patients with missing features were
excluded and 68.10\% of the patients died during the 5.56-year study
period \citep{katzman2018DeepSurv}.

Figure \ref{fig:megaPlot} E, F and Table \ref{tab:megaTable} C
demonstrates the performance over time on a test set. The regression
models (CBLR, Cox) perform best considering IPA, followed by CBNN from
the \(25^{th}\) to \(100^{th}\) percentile of follow-up time. We note a
drop in performance for CBNN before the \(25^{th}\) percentile of
follow-up. For \(C_{IPCW}\), CBNN outperforms the competing models
consistently over follow-up time. Note that performance is similar for
all models, aside from DeepSurv whose \(C_{IPCW}\) is lower than the
rest (Figure \ref{fig:megaPlot} E, F and Table \ref{tab:megaTable} C).

\hypertarget{performance-evaluation-using-the-metabric-dataset}{%
\subsection{Performance evaluation using the METABRIC
dataset}\label{performance-evaluation-using-the-metabric-dataset}}

METABRIC is a 30-yearlong study aiming to discover the molecular drivers
of breast tumors, following 2000 individuals with breast cancer until
death \citep{curtis2012genomic}. They described these growths as primary
invasive breast carcinomas, with the goal of discovering both genetic
and clinical risk factors for breast cancer survival
\citep{curtis2012genomic}. We used the processed dataset made available
through DeepSurv \citep{katzman2018DeepSurv}, which includes 1980
patients of which 57.72\% die due to breast cancer within a median 10
years of follow-up \citep{katzman2018DeepSurv}. There are 10 covariates
in total: time, 4 RNA-Seq gene expressions (MKI67, EGFR, PGR and ERBB2)
and 5 clinical features (hormone treatment indicator, radiotherapy
indicator, chemotherapy indicator, ER-positive indicator and age at
diagnosis) \citep{katzman2018DeepSurv}.

Figure \ref{fig:megaPlot} G, H and Table \ref{tab:megaTable} D shows the
performance on a test set over time. The IPA scores suggest that
regression models outperform competing models on this dataset, as all
the neural network models are equal to or perform worse than KM over
follow-up time. Our CBNN model is comparable to KM until around the
50-percentile of follow-up time, after which CBNN, DeepSurv and DSM drop
in performance. The \(C_{IPCW}\) produces a different ranking. Our CBNN
model outperforms the other models up to the \(25^{th}\) percentile of
follow-up time, whereas DeepHit performs best from the \(50^{th}\) to
\(75^{th}\) percentile of follow-up. With this metric, CBNN and DeepHit
outperform the regression models. This disagreement between IPA and
\(C_{IPCW}\) may be due to the misspecification issue of
concordance-based metrics \citep{cindexfails2019}. The neural network
models may be over-parameterized, as shown by the wide confidence bands
in \(C_{IPCW}\).

\hypertarget{discussion}{%
\section{Discussion}\label{discussion}}

CBNN models survival outcomes by using neural networks on case-base
sampled data. We incorporate follow-up time as a feature, providing a
data-driven estimate of a flexible baseline hazard and time-varying
interactions in our hazard function. The three competing neural network
models we evaluated cannot model time-varying interactions by design
\citep{dsmPaper} \citep{katzman2018DeepSurv} \citep{lee2018DeepHit}. DSM
requires a mixture component distributions for a flexible baseline
hazard to be fit \citep{dsmPaper}. As our goal is to fix the design and
compare performance, we did not change any shared hyperparameters. With
our choice of shared hyperparameters and model design, DSM did not
converge in the complex simulation and SUPPORT case study. Despite this
limitation, we include this method when it did function as it can fit
flexible baseline hazards. Compared to CBNN, the remaining two models
also have limitations. DeepSurv is a proportional hazards model and does
not estimate the baseline hazard \citep{katzman2018DeepSurv}. DeepHit
requires an alpha hyperparameter, is restricted to a single distribution
for the baseline hazard and models the survival function directly
\citep{lee2018DeepHit}. The alternative neural network methods match on
time, while CBNN models time directly.

To assess performance among these models, we use both IPA and
\(C_{IPCW}\) metrics. Concordance-based measures are commonly used in
survival analysis to compare models and we opt to keep them in our
analyses. However, \(C_{IPCW}\) is a non-proper metric and may cause
misspecified models to appear better than they should
\citep{cindexfails2019}. Therefore, we contextualize our \(C_{IPCW}\)
results in relation to IPA, a proper scoring rule. The model rankings
between \(C_{IPCW}\) and IPA differed for both the complex simulation
and METABRIC application. Wider \(C_{IPCW}\) confidence intervals for
DeepHit and DeepSurv show potential misspecification of the models
(Figure \ref{fig:megaPlot} D, H).

With this interpretation of \(C_{IPCW}\) in mind, we assess two
simulations with minimal noise and two case studies in real scenarios.
The simple simulation demonstrates potential pitfalls associated with
neural network models, particularly overparameterization. All neural
network-based approaches performed worse than null KM model (IPA score),
while the regression-based performed better. We attribute this to
potential overparameterization in the neural network models, as the wide
confidence intervals suggest over-fitting, even with dropout. Both
DeepSurv and DSM are affected, while DeepHit and CBNN are less so. This
is a limitation to our strategy for evaluating the methods, which uses a
fixed study design across all assessments.

From this baseline benchmark, we move on to a complex simulation that
requires a method that can learn both time-varying interactions and have
a flexible baseline hazard. Here, CBNN demonstrates a distinct advantage
over all other methods. The regression models show improved performance
over the null KM model, while the competing neural network models
performed worse. Based on our complex simulation results (Figure
\ref{fig:megaPlot} C, D and Table \ref{tab:megaTable} B), CBNN
outperforms the competitors when time-varying interactions and a complex
baseline hazard are present. This simulation shows how CBNN can perform
under ideal conditions, while the following two analyses on real data
serve to assess its performance in realistic conditions.

In the SUPPORT and METABRIC case studies, flexibility in both
interaction modeling and baseline hazard did not improve CBNNs relative
performance, suggesting that this flexibility does not aid prediction in
either case study. From a biological perspective, the 30-year follow-up
time in the METABRIC study may contain competing causes of death. The
causes at the start of the study may not match the causes towards the
end, potentially explaining the drop in performance as we reach later
survival times with fewer individuals. The baseline hazard is unlikely
to cause this drop as DSM is the most flexible competing model in our
comparison. Over-fitting is also unlikely given the tight confidence
intervals. Further research is required to determine the cause of the
drop in performance seen in the METABRIC case study. In both case
studies, CBNN outperforms the competing neural network methods.

The way neural networks are used in this paper are not directly
interpretable. This is by design, as we only wish to compare to
competitors based on predictive ability. All feed-forward neural network
models can make use of different architectures to improve performance or
interpretability. While a deep neural network as used in this paper is
not directly interpretable, techniques such as Local Interpretable
Model-Agnostic Explanations (LIME) \citep{ribeiro2016model} may be used
to interpret the relevance of each feature to a given prediction. An
example of architectural interpretability is a convolutional neural
network for images, which provide positional information about the
relevant pixels in an image \citep{keras}. Note that neither of these
are specific to CBNN. As CBNN is concerned with learning time-varying
interactions, a concern may be interpreting which time-varying
interactions have been learned. Other secondary methods after the
fitting process may be used to retrieve said information
\url{https://ui.adsabs.harvard.edu/abs/2017arXiv170504977T/abstract}.

\section{Conclusions}\label{sec5}

CBNN outperforms all competitors in the complex
simulation, demonstrating its value in survival settings that may
involve time-varying interactions and a complex baseline hazard. Once we
perform case-base sampling and adjust for the sampling bias, we can use
a sigmoid activation function to predict our hazard function. Our
approach simplifies the incorporation of censored individuals, allowing
survival outcomes to be treated as binary ones. Forgoing the requirement
of custom loss functions, CBNN only requires the use of standard
components in machine learning libraries (specifically, the add layer to
adjust for sampling bias and the sigmoid activation function). Due to
the simplicity in its implementation and by extension user experience,
CBNN is both a user-friendly approach to data-driven survival analysis
and is easily extendable to any feed-forward neural network framework.

\hypertarget{data-and-code-availability-statement}{%
\subsection*{Data and code availability
statement}\label{data-and-code-availability-statement}}
\addcontentsline{toc}{section}{Data and code availability statement}

The pre-processed data for the SUPPORT case study can be found at
\url{https://github.com/jaredleekatzman/DeepSurv/tree/master/experiments/data/support}.
The pre-processed data for the METABRIC case study can be found at
\url{https://github.com/jaredleekatzman/DeepSurv/tree/master/experiments/data/metabric}.
The data is accessed using \url{https://github.com/havakv/pycox}. The
code for this manuscript and its analyses can be found at
\url{https://github.com/Jesse-Islam/cbnnManuscript}. The software
package making CBNN easier to use can be found at
\url{https://github.com/Jesse-Islam/cbnn}.

\hypertarget{acknowledgements}{%
\subsection*{Acknowledgements}\label{acknowledgements}}
\addcontentsline{toc}{section}{Acknowledgements}

We would like to thank Dr.~James Meigs, the project
leader of these awards, for his support and helpful discussions. The
work was also supported as part of the Congressionally Directed Medical
Research Programs (CDMRP) award W81XWH-17-1-0347. We would also like to
thank Dr.~James Hanley for his support and discussions while extending
the case-base methodology.

\subsection*{Financial disclosure}

This work was supported by subcontracts from UM1DK078616 and R01HL151855
to Dr.~Rob Sladek.

\subsection*{Conflict of interest}

The authors declare no potential conflict of interests.