\section{Introduction}

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pandemic has highlighted the importance of accurate classification of antibody test results. Most work has focused on labeling data as previously infected (seropositive) or na{\"i}ve.
Due to the deployment of SARS-CoV-2 vaccines in late 2020, there is a clear need for a classification scheme that correctly distinguishes between na{\"ive}, previously infected, and uninfected but vaccinated individuals. However, the traditional diagnostic classification methods of confidence intervals and receiver operating characteristics have no obvious extensions to a multiclass setting.

 Current multiclass applications in
  diagnostic classification are mostly limited to supervised learning and do not address the central role of mathematical modeling in
diagnostics.  Example studies include the application of support-vector machines to automatically sort endomysial autoantibody tests of celiac disease into one of four classes \citep{caetano2019automatic}, and another that trained deep neural networks to label resting-state functional magnetic resonance imaging results with one of six Alzheimer's disease state \citep{ramzan2020deep}.
However, these approaches may not accurately quantify precise training population characteristics or account for the role of prevalence,  both of which should inform the classification procedure. In contrast, modeling can overcome this limitation. Binary (two class) examples include two-dimensional (2D) modeling of antigen targets coupled with optimal decision theory \citep{patrone2021classification}, statistical modeling applied to either antibody or viral-load tests   \citep{bottcher2022statistical}, and an approach to the time-dependent problem for antibody measurements \citep{bedekar2022prevalence}. However, none of these works discussed multiclass extensions.

This paper uses mathematical modeling to fully address the task of multiclass classification in the context of diagnostic testing. We begin by showing that the notion of \textit{generalized prevalence}--the relative fraction of the population in each class--is fundamental for defining our objective function, the convex combination of false classification rates (Section \ref{sec:classif}). Minimization thereof yields optimal classification. Interestingly, we show that these prevalences can be computed without classification by solving a linear system.
 We validate our methods using a SARS-CoV-2 serological data set with na{\"i}ve, previously infected, and vaccinated classes \citep{ainsworth2020performance,wei2021antibody}\footnote{Certain commercial equipment, instruments, software, or materials are
identified in this paper in order to specify the experimental procedure adequately. Such identification is not intended to imply recommendation or endorsement by the National Institute of Standards and Technology, nor is it intended to imply that the materials or equipment identified are necessarily the best available for the purpose.} (Section \ref{sec:1D_ex}).
We then computationally validate the convergence of  generalized prevalence estimates to the true values in mean square and illustrate a generalization to 2D data (Section \ref{sec:comp}). Finally, the discussion includes further analysis of prevalence estimation,  extensions, and limitations (Section \ref{sec:disc}).

\section{Notation}

This work combines set and measure theory with applied diagnostics; readers will likely not be experts in both. In order for the ideas of this paper to be readily implemented by diagnostics experts and the applications understood by mathematicians, we provide baseline terminology from both fields.

\subsection{ Definitions from applied diagnostics}

\begin{itemize}

\item The na{\"i}ve class comprises individuals that have not been previously infected or vaccinated. In a binary classification, such samples are often referred to as `negative'.

\item The previously infected class comprises individuals with a prior infection but who are unvaccinated. In a binary classification, such samples are often referred to as `positive'.

\item The vaccinated class comprises individuals who have been inoculated against a disease without a prior infection.

\item  Training data correspond to samples for which the true classes are known.  Typically, such data are used to construct conditional probability models.

\item Test data correspond to samples for which the true classes are unknown or assumed to be unknown for validation purposes. Typically, a classification procedure is applied to such data.

\item  Generalized prevalence is the relative fraction of samples in a population that belong to each class.

\end{itemize}

\subsection{Definitions from measure theory}

\begin{itemize}

\item A set is a collection of objects, e.g.\ measurement values. A domain is a set in a continuous measurement space; see Figure  \ref{fig:phase} for an example.

\item The symbol $\mathbb{R}$ denotes the set of all real numbers. The symbol $\mathbb{R}^m$  denotes the real coordinate space of dimension $m$ consisting of all real-valued vectors of length $m$.

\item  The symbol $\in$ indicates set inclusion. The expression $\bm{r} \in A$ means $\bm{r}$ is in set $A$.

\item The symbol $\subset$ denotes the subset relationship of two sets. The expression $A \subset B$ means that all elements in $A$ are contained in $B$.

\item The use of a superscript $C$ denotes the complement of a set. The set $D^C$ contains all elements in the measurement space not in $D$.

\item  The symbol $\emptyset$ denotes the empty set, which contains no elements.

\item  The operator $\cup$ denotes the union of two sets. The set $C = A \cup B$ contains all elements in either $A$ or $B$ or both.

\item The operator $\cap$ denotes the intersection of two sets. The set $C = A \cap B$ contains all elements in both $A$ and $B$.

\item  The operator $\setminus$ denotes the set difference. The set $C = A \setminus B$ contains all objects in $A$ that are not also in $B$.
An equivalent interpretation is that $A \setminus B$ is the result of removing the common elements of $A$ and $B$ from $A$.

\item  The notation $A = \{\bm{r} : *\}$ defines the set $A$ as all $\bm{r}$ that satisfy condition $*$.

\end{itemize}

\subsection{Notation specific to this paper}

\begin{itemize}

\item The set $\Omega$ denotes the entire measurement space.

\item The label $C_j$ refers to the $j$th class.

\item The  generalized prevalence for class $C_j$ is denoted by $q_j$.

\item The set $D_j$ denotes a domain corresponding to $C_j$.

\item The use of a superscript $\star$ denotes an optimal quantity. For example, $D_j^{\star}$ could be an optimal classification domain corresponding to class $C_j$.

\end{itemize}

\section{Generalized prevalence estimation and multiclass classification}
\label{sec:classif}

Prevalence estimation and classification rely on the same framework of antibody measurements. For each individual or sample, we represent corresponding measurements as a vector $\bm{r} = (r_{1}, \ldots r_{m} ) \in \Omega \subset \mathbb{R}^m$.  Here, $\bm{r}$ could denote $m$ antibody types targeting different parts of a virus as measured in median fluorescence intensity (MFI). Let $P_j(\bm{r})$ describe the probability that a sample from class $C_j$ yields measurement value $\bm{r}$.
 These conditional probability density functions (PDFs) are assumed known in this section; their construction is considered for example serological data in Section \ref{sec:pdfs}.

 The generalized  prevalence $q_j$ is the relative fraction of the population corresponding to class $C_j$. In what follows we assume there are $n$ classes. The generalized prevalences must sum to 1, which implies
\begin{subequations}
\begin{eqnarray}
\sum_{j = 1}^n q_j = 1,
\label{eq:prev_cond_1} \\
q_k = 1 - \sum_{\substack{j = 1 \\ j \neq k}}^n q_j, \quad \text{ for } k \in \{1, \ldots, n\}.
\label{eq:prev_cond_2}
\end{eqnarray}
\label{eq:q_relation}
\end{subequations}
The probability density $Q(\bm{r})$ of a measurement $\bm{r}$ for a test sample is given by
\begin{equation}
Q(\bm{r}) = \sum_{j = 1}^n q_j P_j(\bm{r}).
\end{equation}
The product $q_j P_j(\bm{r})$ is the probability that a random sample both belongs to class $C_j$ and has measurement value $\bm{r}$;
thus, the expression for $Q$ is an instance of the law of total probability.
This quantity plays an important role in prevalence estimation and classification.

\subsection{Generalized prevalence estimation}
\label{sec:pool}
\label{sec:prev}

To demonstrate the importance of prevalence in diagnostic classification, consider the United States population's SARS-CoV-2 antibody response. In early 2020, most samples should have been classified as na{\"i}ve because the disease prevalence was small. By February 2022,
the disease prevalence was estimated at 57.7 \% \citep{clarke2022seroprevalence}; a significant fraction of samples should have been classified as previously infected. Crucially, \textit{the same measurement value may be classified differently depending on the disease prevalence}. This example shows that prevalence plays an integral role in classifying diagnostic tests and should be estimated before classification of test data.
We address this need by designing unbiased estimators for the prevalences $\{q_j\}$.

For $n$ classes, consider a partition $\{D_j\}$
that separates the measurement space $\Omega$ into $n$ nonempty domains. It is important to note that these $D_j$ are not classification domains.
Define
\begin{equation}
Q_{j} = \int_{D_j} Q(\bm{r}) d \bm{r} = \int_{D_j}  \sum_{k = 1}^n q_{k} P_k(\bm{r}) d \bm{r} = \sum_{k = 1}^n q_{k} \int_{D_j} P_k(\bm{r}) d \bm{r} = \sum_{k = 1}^n  P_{j,k}q_{k},
\end{equation}
 where
\begin{equation}
P_{j,k} = \int_{D_j} P_k(\bm{r}) d \bm{r}.
\end{equation}
 Writing this as a linear system yields
\begin{equation}
\left[ \begin{array}{c}
Q_1 \\ \vdots \\ Q_n
\end{array} \right] = \left[ \begin{array}{ccc}
P_{1,1} &  \ldots & P_{1,n} \\
\vdots & \ddots & \vdots \\
P_{n,1} & \ldots & P_{n,n}
\end{array} \right] \left[ \begin{array}{c}
q_1 \\ \vdots \\ q_n
\end{array} \right].
\label{eq:q_full}
\end{equation}
Taking $k = n$ in (\ref{eq:prev_cond_2}) implies
\begin{equation}
\left[ \begin{array}{c}
Q_1 \\ \vdots \\ Q_{n-1}
\end{array} \right]  - \left[ \begin{array}{c}
P_{1,n} \\ \vdots \\  P_{n-1, n}
\end{array} \right] = \left( \left[ \begin{array}{ccc}
P_{1,1} &  \ldots & P_{1, n-1}  \\
\vdots & \ddots & \vdots \\
P_{n-1,1} & \ldots & P_{n-1, n-1}
\end{array} \right] - \left[ \begin{array}{c}
P_{1,n}  \\ \vdots \\   P_{n-1,n}
\end{array} \right] \underbrace{[1, \ldots, 1]}_{n-1} \right) \left[ \begin{array}{c}
q_1 \\ \vdots \\ q_{n-1}
\end{array} \right].
\end{equation}
This yields the prevalences $\bm{q}$ as the solution to the system
\begin{subequations}
\begin{eqnarray}
\bm{q} = (\bm{P} - \bm{P_n})^{-1} \left( \bm{\overline{Q}} - \overline{\bm{P_n}} \right),
\label{eq:prev_est_m} \\
q_k \geq 0 \quad \text{ for } k = 1, 2, \ldots, n-1,
\label{eq:prev_est_const}
\end{eqnarray}
\label{eq:prev_est_eqn}
\end{subequations}
\noindent where $\bm{q}$ is the vector of length $n-1$ whose $j$th entry is $q_j$, $\bm{P}$ is the $(n-1) \times (n-1)$ matrix whose $(i,j)$th entry is $P_{i, j}$, $\bm{P_n}$ is the $(n-1) \times (n-1)$ matrix whose $(i,j)$th entry is $P_{i, n}$, $\bm{\overline{Q}}$ is the vector of length $n-1$ whose $j$th entry is $Q_j$, and $\overline{\bm{P_n}}$ is the vector of length $n-1$ whose $j$th entry is $P_{j, n}$. The last prevalence $q_n$ is found via  (\ref{eq:prev_cond_2}) with $k = n$.  We assume that the inverse of the matrix $\bm{P} - \bm{P_n}$ exists; Section \ref{sec:disc_lim_prev} further discusses the matrices $\bm{P}$ and $\bm{P} - \bm{P_n}$.

To estimate the generalized prevalences, estimate the $Q_j$ by $\hat{Q}_j$, where
\begin{equation}
 Q_j \approx \hat{Q}_j = \frac{1}{S} \sum_{i = 1}^S \mathbb{I}( \bm{r}_i \in D_j).
 \label{eq:Q_approx}
 \end{equation}
 Here, $S$ is the total number of samples and $\mathbb{I}$ denotes the indicator function. Substituting $\hat{Q}_j$ for $Q_j$ in (\ref{eq:prev_est_eqn}) yields an estimate $\hat{q}_k$ for $q_k$. When the PDFs $P_j(\bm{r})$ are known and $(\bm{P} - \bm{P_n})^{-1}$ exists, these estimates $\hat{q}_k$ are unbiased, i.e., $E[\hat{q}_k] = q_k$. This follows directly from the fact that $\hat{Q}_j$ is a Monte Carlo estimator of $Q_j$. Further, the  generalized prevalence estimates converge to the true values in mean square as the number of samples is increased \citep{caflisch1998monte}. This is illustrated in Section \ref{sec:conv_prev_est}.

We note that the  generalized prevalence estimates are not unique due to the arbitrariness of the $\{D_j\}$.
However, the non-uniqueness allows us to select any reasonable partition over which to find the estimators $\{ \hat{Q}_j\}$. This is discussed further in Section \ref{sec:disc_lim_prev}.

\subsection{Optimal classification}
\label{sec:opt_class}

Our task is to define a partition $\{D_j\}$ (not necessarily the same as for prevalence estimation) of the measurement space $\Omega$ such that each domain corresponds to one and only one class $C_j$. A measurement $\bm{r}$ is assigned to class $j$ if $\bm{r} \in D_j$. We require
\begin{subequations}
\begin{eqnarray}
 \mu_{j}\left(\bigcup_{k = 1}^n D_k \right) = 1 \quad \forall j \in \{1, \ldots, n\},
 \label{eq:req2} \\
 \mu_{\ell}(D_j \cap D_k) = 0 \text{ for } j \neq k, \quad \forall \ell = 1, 2, \ldots, n
 \label{eq:req1}
\end{eqnarray}
\label{eq:req}
\end{subequations}
where $ \mu_{\ell}(X) = \int_X P_{\ell}(\bm{r}) d\bm{r}$. Here, (\ref{eq:req2}) ensures that any sample can be classified and (\ref{eq:req1}) enforces single-label classification (up to sets of measure zero). To identify an optimal partition $\{D_j^{\star}\}$, we construct the loss function
\begin{equation}
\mathscr{L}(D_1, \ldots, D_n) = \sum_{j = 1}^n q_j \int_{\Omega \setminus D_j} P_j(\bm{r}) d\bm{r}.
\label{eq:loss}
\end{equation}
Here, (\ref{eq:loss}) is the  generalized prevalence-weighted convex combination of false classification rates as a function of the domains $D_j$. Intuitively, we expect that a sample with measurement $\bm{r}$ should be assigned to the class domain $D_j$ to which it has the highest probability of belonging; that is, the highest value of $q_j P_j(\bm{r})$ for all $j$. Accordingly, the loss function (\ref{eq:loss}) penalizes misclassified measurements $\bm{r}$ with high probability values.

To address situations in which a measurement has an equal highest probability of belonging to two or more classes, we introduce the following set for each class $C_j$:
\begin{equation}
\mathscr{E}_j =  \bigcup_{\substack{k = 1 \\ k \neq j}}^n \{ \bm{r}:  q_k P_k(\bm{r}) =  q_j P_j(\bm{r}) = \max_i q_i P_i(\bm{r})  \} .
\end{equation}
In most practical implementations all $ \mathscr{E}_j$ have measure zero, and
the domains
\begin{equation}
D_j^{\star} = \{\bm{r} : q_j P_j(\bm{r}) > q_k P_k (\bm{r}) \text{ for } k \neq j\}
\label{eq:opt_d}
\end{equation}
minimize the loss function $\mathscr{L}$ up to sets of measure zero.  The proof is shown in \ref{sec:app_a} and involves a straightforward application of set theory; see also \cite{williams2006gaussian} for similar ideas.

If $\mathscr{E}_j$ has nonzero measure, randomly assigning a measurement in $\mathscr{E}_j$ to one of the classes to which it has equal maximal probability of belonging does not effect the loss $\mathscr{L}$.
In this case, the optimal domains are generalized to
\begin{equation}
D_j^{\star} = \{\bm{r} : q_j P_j(\bm{r}) > q_k P_k (\bm{r}) \text{ for } k \neq j\} \cup Z_{\mathscr{E}_j},
\end{equation}
where $Z_{\mathscr{E}_j}$ is an element of a partition of $\mathscr{E}_j$ that we define iteratively as follows.
\begin{subequations}
\begin{eqnarray}
 Z_{\mathscr{E}_1} = \mathscr{E}_1, \\
 Z_{\mathscr{E}_k} = \mathscr{E}_k \mathbin{\big \backslash} \bigcup_{j = 1}^{k-1} \mathscr{E}_j, \quad k \in \{ 2, \ldots, n\} .
 \end{eqnarray}
 \end{subequations}
 This ensures that no measurement in a set $\mathscr{E}_j$ is assigned to more than one optimal domain. Note that by construction, $Z_{\mathscr{E}_n}$ is empty.

 Figure \ref{fig:phase} shows a 2D conceptual illustration of $\{ \mathscr{E}_j\}$, which are the lines delineating the optimal regions. In 2D, line segments have Lebesgue measure zero. Thus, classification follows  (\ref{eq:opt_d}).  Note that the ``multipoint" at which the lines meet has equal probability of belonging to all three classes. We discuss this further in Section \ref{sec:loc_acc}.

\begin{figure}[t]
\centering
\includegraphics[scale=.5]{7_19_22_phase.eps}
\caption{Illustration of classification domains $D_1, D_2$, and $D_3$ in which the sets of equal probability of a measurement belonging to two or more classes are shown as lines separating the optimal regions.}
\label{fig:phase}
\end{figure}

\section{Example applied to SARS-CoV-2 antibody data with three classes}
\label{sec:1D_ex}

To demonstrate the concepts developed in Section \ref{sec:classif}, we apply our methods to serological data with three classes. Publicly available data sets associated with \cite{ainsworth2020performance} and \cite{wei2021antibody} provide previously infected, na{\"i}ve, and vaccinated antibody measurements. The vaccine data \citep{wei2021antibody} are recorded for individuals that were innoculated with one of two vaccines. We refer to these as Vaccine A and Vaccine B and analyze the populations separately and together.  The studies provide SARS-CoV-2 anti-spike immunoglobulin (IgG) antibody measurements; see \ref{sec:app_MFI}  for measurement details. We use one-dimensional (1D) data to illustrate a straightforward multiclass example; Section \ref{sec:2D} demonstrates that our analysis holds for higher measurement dimensions.

All data are transformed to a logarithmic scale as follows:
\begin{equation}
r = \log_2(\tilde{r} + 2) - 1.
\label{eq:log_transform}
\end{equation}
Here, $\tilde{r}$ and $r$ represent the original and log-transformed values; $\tilde{r}$ has units of ng/mL and $r$ is nondimensional.
This transformation puts the data on the scale of bits and allows for better viewing of measurements that range over several decades of MFI values in the original units. \cite{wei2021antibody} truncated vaccinated samples, with lower and upper transformed limits of 1 and roughly 8.

Figure \ref{fig:1D_hist_data} shows a histogram of the data with the vaccinated category split by vaccine manufacturer (Figure \ref{fig:1D_hist_split}) and combined (Figure \ref{fig:1D_hist_all}).
Previously infected samples have the largest spike IgG antibody levels and na{\"i}ve samples the smallest; the vaccinated class falls in the middle. The vaccinated class overlaps with some na{\"i}ve and previously infected samples. Due to the truncation of vaccinated measurements, the corresponding right-most histogram bin contains many samples.
We separate the data into randomly generated training (80 \% of samples) and test (20 \%) populations.

\begin{figure}[t]
\centering
\subfloat[][Vaccine A and Vaccine B split]{\includegraphics[scale=.5]{3_18_22_3class_hist_2.eps}
\label{fig:1D_hist_split}}
\subfloat[][Vaccine A and B combined]{\includegraphics[scale=.5]{3_18_22_3class_hist_all.eps}
\label{fig:1D_hist_all}}
\caption{Histograms of previously infected, na{\"i}ve, and vaccinated data from Ainsworth et. al (2020) and Wei et al. (2021).
}
\label{fig:1D_hist_data}
\end{figure}

\subsection{Conditional probability distributions}
\label{sec:pdfs}

We fit probability distributions to the training data to model the na{\"i}ve, previously infected, and vaccinated antibody responses. For our purposes, we assume these are distinct classes; a sample belongs to one and only one of the three categories.
To construct the conditional PDF for each population, we select a parameterized model that empirically characterizes the shape and spread of the samples.
We determine parameters separately for the three training populations by maximum likelihood estimation (MLE).
The na{\"i}ve training population is fit to a Burr distribution
\begin{equation}
N(r) = \frac{ck}{\lambda} \left( \frac{r}{\lambda} \right)^{c-1} \left[ 1 +\left( \frac{r}{\lambda}\right)^c \right]^{-k-1},
\label{eq:N_eqn}
\end{equation}
which describes a right-skewed sample population.

The previously infected  training population is fit to a stable distribution described by characteristic function
\begin{equation}
\phi(r) = \exp\left\{i r \delta - |\gamma r|^{\alpha} \left[1 + \frac{2i}{\pi} \beta \text{sgn}(r) \log(\gamma r)\right] \right\}
\label{eq:P_eqn}
\end{equation}
for $\gamma \neq 1$.
Here, $i$ is the imaginary unit and sgn is the sign function, which returns +1, -1, or 0. This distribution describes a left-skewed sample population.

We fit the vaccinated training populations to an extreme-value distribution after observing the mostly symmetric shape of the data with a spike at the right truncation limit:
\begin{equation}
V(r) = \frac{1}{\sigma} \exp \l \frac{r - \mu}{\sigma} \r \exp \left[ - \exp \l \frac{r - \mu}{\sigma} \r \right].
\label{eq:V_eqn}
\end{equation}
We apply data censoring to better fit the truncated data; this is described in \ref{sec:app_trunc}.

\begin{figure}[h]
\centering
\includegraphics[scale=.5]{8_24_22_vaxA_pdf.eps}
\caption{Conditional PDFs for the three na{\"i}ve, previously infected, and vaccinated classes trained on the training data for a Vaccine A visualization of the vaccinated class. See Supplemental Figure S1 for the PDFs of the Vaccine B and combined visualizations. }
\label{fig:pdf}
\end{figure}

The analysis in this section is identical for all three visualizations of the vaccinated class. In what follows, we report all results but show only the Vaccine A figures as examples. Corresponding figures for the Vaccine B and combined visualizations of the vaccine class are left to the Supplemental Data.

Figure \ref{fig:pdf} shows the conditional PDFs, represented as continuous curves, trained on the three-class training data with a Vaccine A vaccinated class. The blue, red, and black curves correspond to the na{\"i}ve, previously infected, and vaccinated models.
The effect of truncating the data at the upper limit is visible in the right-most bin of the vaccinated class histogram; this is accounted for by the data-censoring. As a result, the vaccinated class PDF exhibits spikes at the upper and lower truncation values. This spike is an artifact of the original data collection process and not a typical problem.

\subsection{Generalized prevalence estimation}

Recall that prevalence estimation of test data requires a partition that separates the measurement space $\Omega$ into  $n$ nonempty domains. Here, the number of classes is $n = 3$. We create a partition using $k$-means clustering with $k = 3$, which assigns each measurement to the cluster with the closest mean. Figure \ref{fig:part} shows the partition for our test data set with a Vaccine A vaccinated class. The clustering separates the three populations reasonably well; see Section \ref{sec:disc_lim_prev} for the importance of this statement. The partition need not perfectly separate the data by class to estimate prevalences with high accuracy.
We estimate  generalized prevalences for the test data via (\ref{eq:prev_est_eqn}) and record true and estimated values in Table \ref{table:prev_est}.

\begin{figure}[h]
\centering
\includegraphics[scale=.5]{5_24_22_clust_AZ.eps}
\caption{Test data $k$-means partitioning with a Vaccine A vaccinated class for  generalized prevalence estimates. We use $k = 3$ classes; the clustered domains are labeled as $D_1, D_2$, and $D_3$. See Supplemental Figure S2 for the partitions of the Vaccine B and combined visualizations of the vaccinated class.}
\label{fig:part}
\end{figure}

\begin{table}[h]
\centering
\begin{tabular}{|l|r|rrr|rrr|r|}
\hline
 & \multicolumn{1}{l|}{} & \multicolumn{3}{c|}{\textbf{\begin{tabular}[c]{@{}c@{}}Vaccine data set\\ Estimated (true)  generalized prevalence\end{tabular}}} & \multicolumn{3}{c|}{\textbf{\begin{tabular}[c]{@{}c@{}}Errors\\ (\%)\end{tabular}}} & \multicolumn{1}{c|}{\textbf{\begin{tabular}[c]{@{}c@{}}Avg.\\ errors\\ (\%)\end{tabular}}} \\ \hline
 & \multicolumn{1}{l|}{\textbf{}} & \multicolumn{1}{c|}{\textit{A}} & \multicolumn{1}{c|}{\textit{B}} & \multicolumn{1}{c|}{\textit{All}} & \multicolumn{1}{c|}{\textit{A}} & \multicolumn{1}{c|}{\textit{B}} & \multicolumn{1}{c|}{\textit{All}} &  \\ \hline
\multirow{3}{*}{\textbf{Class}} & $N$ & \multicolumn{1}{r|}{0.523 (0.521)} & \multicolumn{1}{r|}{0.538 (0.531)} & \multicolumn{1}{r|}{0.445 (0.444)}  & \multicolumn{1}{r|}{0.377}& \multicolumn{1}{r|}{1.02} &  {0.223}  & 0.540 \\ \cline{2-9}
 & $P$ & \multicolumn{1}{r|}{0.300 (0.286)}  & \multicolumn{1}{r|}{0.313 (0.292)} &  {0.285 (0.244)} & \multicolumn{1}{r|}{4.69} & \multicolumn{1}{r|}{7.10} & {17.0} & 9.93 \\ \cline{2-9}
 & $V$ & \multicolumn{1}{r|}{0.177 (0.193)} & \multicolumn{1}{r|}{0.149 (0.177)} & {0.270 (0.312)} & \multicolumn{1}{r|}{7.99} & \multicolumn{1}{r|}{15.0} & {13.6} & 12.2 \\ \hline
\textbf{Avg.} &  & \multicolumn{1}{r|}{} & \multicolumn{1}{r|}{} &  & \multicolumn{1}{r|}{4.74}  & \multicolumn{1}{r|}{7.67} & {10.3} & 7.55 \\ \hline
\end{tabular}
\caption{ Estimated and true  generalized prevalences for the test data na{\"i}ve (N), previously infected (P), and vaccinated (V) classes.  Vaccine A (A) and Vaccine B (B) are considered separately and together (All).}
\label{table:prev_est}
\end{table}

\subsection{Optimal classification}

We classify the training data using known generalized prevalences via  (\ref{eq:opt_d}).
Figure \ref{fig:train} shows the optimal domains, labeled $D_N^{\star}$, $D_V^{\star}$, and $D_P^{\star}$, for a Vaccine A vaccinated class.
 For this 1D example with three classes, the optimal classification domain boundaries can be represented by upper and lower threshold levels. Samples with measurements below the smaller level are classified as na{\"i}ve, samples with measurements between the thresholds as vaccinated, and samples with measurements above the larger level as previously infected. All three populations have overlapping PDFs, which reduces classification accuracy.

 Accurate classification of test data is possible with reasonably close prevalence estimates. We classify the test data using estimated generalized prevalences and display the optimal classification domains for a Vaccine A vaccinated class in Figure \ref{fig:test_opt}.

\begin{figure}[h]
\centering
\subfloat[][Training data]{
\includegraphics[scale=.5]{8_24_22_train_opt_AZ.eps}\label{fig:train}}
\subfloat[][Test data]{\includegraphics[scale=.5]{8_24_22_test_opt_AZ.eps}\label{fig:test_opt}}
\caption{Training (a) and test (b) data with a Vaccine A vaccinated class with optimal decision thresholds using a known prevalence.
Vertical dashed lines indicate the optimal decision boundaries. The optimal na{\"i}ve, vaccinated, and previously infected domains are labeled $D_N^{\star}$, $D_V^{\star}$, and $D_P^{\star}$. See Supplemental Figures S3 and S4 for the optimal domains for Vaccine B and combined visualizations of the vaccinated class.
}
\label{fig:train_test}
\end{figure}

\begin{table}[h]
\centering
\begin{tabular}{|c|r|r|}
\hline
\multicolumn{1}{|l|}{} & \multicolumn{1}{c|}{\textbf{\begin{tabular}[c]{@{}c@{}}Training\\ classification \\ error (\%)\end{tabular}}} & \multicolumn{1}{c|}{\textbf{\begin{tabular}[c]{@{}c@{}}Test\\ classification\\ error (\%)\end{tabular}}} \\ \hline
\textbf{Vaccine A} & 7.87  & 5.08 \\ \hline
\textbf{Vaccine B} & 7.07 &  5.46 \\ \hline
\textbf{Combined} & 7.90 & 4.78 \\ \hline
\end{tabular}
\caption{Classification errors for training data using a known prevalence and test data with an estimated prevalence. Vaccine A and Vaccine B are considered separately and together (Combined).}
\label{table:train_test_classif}
\end{table}

Training and test data classification errors are recorded in Table \ref{table:train_test_classif}. Taken over the three considerations of the vaccinated class, the average error for the training data is 7.61 \% and the same for the test data is 5.11 \%.

\section{Computational validation}
\label{sec:comp}

We numerically demonstrate two important features of our generalized prevalence estimation and multiclass optimal classification procedures. First, we show the convergence of our prevalence estimates to the true values as the number of samples is increased. Second, we present a 2D tri-class problem to show how the method generalizes to higher dimensional measurement spaces.

\subsection{Convergence of prevalence estimates}
\label{sec:conv_prev_est}

We use our probability models (\ref{eq:N_eqn})-(\ref{eq:V_eqn}) to generate synthetic data sets whose relative frequencies match the  generalized prevalences of the \cite{ainsworth2020performance} and \cite{wei2021antibody} data.
We systematically increase the number of synthetic data points used while holding generalized prevalences fixed to study the effect of sample size on prevalence convergence. For each number of points used, we generate 1000 synthetic data sets and compute statistics on our results.

\begin{figure}[t]
\centering
\subfloat[][1000 synthetic samples]{\includegraphics[scale=.5]{8_26_22_syn_vax_A_1000.eps}
\label{fig:syn_AZ_data}}
\subfloat[][ Boxchart of estimate statistics]{\includegraphics[scale=0.5]{5_24_22_prev_est_converge_AZ.eps}
\label{fig:syn_AZ_boxchart}} \\
\subfloat[][ Convergence of the error in mean square]{\includegraphics[scale=0.5]{5_24_22_prev_est_conv_rate_AZ.eps}
\label{fig:syn_AZ_conv}}
\begin{minipage}[b]{18em}
\caption{Prevalence estimation convergence for synthetic data using a Vaccine A vaccinated class (1000 simulations).
 In (b), the boxes display the median and upper and lower quartiles as the line inside the box and its top and bottom edges. The whiskers show non-outlier maximum and minimum values; outliers vary from the median by more that 1.5 times the difference between the upper and lower quartiles, and are shown as circles. In (b) and (c), the subscripts $N$, $P$, and $A$ denote na{\"i}ve, previously infected, and Vaccine A vaccinated. The number of samples is $S$.}
\label{fig:syn_AZ}
\end{minipage}
\end{figure}

\begin{table}[t]
\centering
\begin{tabular}{|l|r|r|r|r|}
\hline
 & \multicolumn{1}{l|}{} & \multicolumn{1}{c|}{\textbf{Na{\"i}ve}} & \multicolumn{1}{c|}{\textbf{Previously infected}} & \multicolumn{1}{c|}{\textbf{Vaccine A}} \\ \hline
\textbf{\begin{tabular}[c]{@{}l@{}}True generalized\\ prevalences\end{tabular}} & \multicolumn{1}{l|}{} & 0.521 & 0.286 & 0.193 \\ \hline
\multirow{4}{*}{\textbf{Number of samples}} & $10^2$ & 0.518 $\pm$ 0.0212 & 0.289 $\pm$ 0.0211 & 0.193 $\pm$ 0.0299 \\ \cline{2-5}
 & $10^3$ & 0.522 $\pm$ 0.0067 & 0.286 $\pm$ 0.0061 & 0.193 $\pm$ 0.0093 \\ \cline{2-5}
 & $10^4$ & 0.522 $\pm$ 0.0021 & 0.285 $\pm$ 0.0020 & 0.193 $\pm$ 0.0030 \\ \cline{2-5}
 & $10^5$ & 0.522 $\pm$ 0.0007 & 0.285 $\pm$ 0.0007 & 0.193 $\pm$ 0.0011 \\ \hline
\end{tabular}
\caption{ Generalized prevalence estimate means and standard deviations taken over 1000 simulations of synthetic data generated from the probability models for increasing number of samples for a Vaccine A vaccinated class.}
\label{table:prev_conv}
\end{table}

Figure \ref{fig:syn_AZ} shows our analysis for the Vaccine A vaccinated class.
Figure \ref{fig:syn_AZ_boxchart} shows a boxchart of the statistics for using $10^2, 10^3, 10^4$, and $10^5$ samples. The estimates have more outliers and variation when few samples are used, which decreases as the number of points is increased. Even with few samples, the median  generalized prevalence estimates are close to the true  generalized prevalences. Table \ref{table:prev_conv} records the mean and standard deviations of our results. Even for only 1000 samples, our estimates agree with the true  generalized prevalences with roughly 2 \% relative error.

Figure \ref{fig:syn_AZ_conv} plots the standard deviation of the  prevalence estimate error on a log-log scale against the number of samples.
The standard deviation should decrease with the inverse square root of the number of samples \citep{caflisch1998monte}, which is plotted for comparison. Our empirical convergence rates all agree with the theory through 10,000 samples; the rate is maintained for the Vaccine A vaccinated class.

\subsection{Generalization to higher dimensions}
\label{sec:2D}

\begin{figure}[t]
\centering
\subfloat[][Probability model contours and synthetic data]{\includegraphics[scale=.5]{9_28_21_RBD_N_3_cat_contours.eps}}
\subfloat[][$k$-means clustering]{\includegraphics[scale=.5]{9_28_21_RBD_N_3_class_clust.eps}} \\
\subfloat[][Optimal classification domains]{
\includegraphics[scale=.5]{9_28_2_RBD_N_3_cat_opt_data.eps}}
\begin{minipage}[b]{18em}
\caption{(a) Level sets of conditional PDFs with example synthetic data, (b) $k$-means clustering, (c) optimal classification domains with estimated generalized prevalences. In (c), the subscripts $N$, $P$, and $V$ denote na{\"i}ve, previously infected, and vaccinated.
}
\label{fig:3_contours}
\end{minipage}
\end{figure}

We now explore a 2D synthetic numerical validation of  generalized prevalence estimation and multiclass optimal classification. See \cite{luke2022improving} for a discussion of the implications of higher-dimensional modeling on diagnostic testing accuracy. The synthetic values we use are modeled off the receptor-binding domain (RBD) and nucleocapsid (N) SARS-CoV-2 antibody targets; together these form a measurement double $\bm{r}$. Details about the models and information about the data are given in \ref{sec:app_b}. Figure \ref{fig:3_contours}a shows an example of 2D synthetic antibody measurements with na{\"i}ve, previously infected, and vaccinated classes with true  prevalences of 0.3, 0.2, and 0.5. The conditional PDFs are shown as contour lines of constant probability. We use 1000 total synthetic samples.

To quantify uncertainty in the prevalence estimates, we randomly generate 1000 synthetic sets of samples using fixed prevalences. We then partition the measurement space via $k$-means clustering using one synthetic sample set (see Figure \ref{fig:3_contours}b), fix the partition, and use (\ref{eq:prev_est_eqn}) to generate prevalence estimates for all sets.
The results are shown in Table \ref{table:q_stats}.  Figure \ref{fig:q_hist} shows histograms of the  generalized prevalence estimates and true values, which fall within the middle of each distribution.  We classify using these estimated prevalences via (\ref{eq:opt_d}) and find an average error of 1.58 \%. Figure \ref{fig:3_contours}c shows example optimal classification domains. The gold region is the previously infected domain, the purple is the vaccinated, and the remainder of the measurement space, colored in light blue, defines the na{\"i}ve domain. For this example,
the false classification rate is  1.8 \%.

\begin{table}[h]
\centering
\begin{tabular}{|c|r|r|r|r|}
\hline
\multicolumn{1}{|l|}{} & \multicolumn{1}{|l|}{\textbf{True val}} & \multicolumn{1}{c|}{$\bm{\mu}$} & \multicolumn{1}{c|}{$\bm{\sigma}$} & \multicolumn{1}{c|}{\textbf{CV} $\bm{\left( \frac{\sigma}{\mu}\right)}$} \\ \hline
$\bm{q_1}$     & 0.3        & 0.300                          & $9.6 \times 10^{-3}$   & 0.0319                                                    \\ \hline
$\bm{q_2}$     & 0.2        & 0.200                          & $7.4 \times 10^{-3}$                                          & 0.0371                                                     \\ \hline
$\bm{q_3}$    & 0.5         & 0.500                           & $6.3 \times 10^{-3}$                               & 0.0126                                                     \\ \hline
\end{tabular}
\caption{ Statistics for the data shown in Figure \ref{fig:q_hist}. The true values of the prevalences are given along with the mean $\mu$, standard deviation $\sigma$, and coefficient of variation (CV) of the estimates. }
\label{table:q_stats}
\end{table}

\begin{figure}[t]
\centering
\subfloat[][Na{\"i}ve]{\includegraphics[scale=.5]{7_22_22_q1_hist.eps}}
\subfloat[][Previously infected]{\includegraphics[scale=.5]{7_22_22_q2_hist.eps}} \\
\subfloat[][Vaccinated]{\includegraphics[scale=.5]{7_22_22_q3_hist.eps}}
\begin{minipage}[b]{18em}
\caption{(a-c): Histograms of  generalized prevalence estimates from 1000 synthetic data sets. Ten bins are used for each histogram and the true prevalence is shown as a vertical red line.}
\label{fig:q_hist}
\end{minipage}
\end{figure}

\section{Discussion}
\label{sec:disc}

\subsection{Limiting cases of prevalence estimation and implications for assay design}
\label{sec:disc_lim_prev}

An interesting observation of our prevalence estimation scheme is that the structure of the matrix underpinning the linear system encodes information about overlap between populations. As such, the matrix potentially informs best practices for prevalence estimation.
Further, our method extends the binary procedure of \cite{patrone2021classification} and may provide insight into the simpler setting. Here we examine limiting cases of prevalence estimation and connect characteristics of the matrix $\bm{P}$ to assay accuracy.

 We explore interpretations of equivalent definitions of singularity of the matrix $\bm{P}$. Recall that the quantity $P_{j,k}$ gives the probability density of class $k$ falling in domain $D_j$. If all elements of a row (column) of the matrix $\bm{P}$ are zero, the probability of any measurement value falling in (belonging to) the corresponding domain (class) is zero. If the columns of $\bm{P}$ are linearly dependent, the probability of a sample belonging to class $C_k$ having a measurement in domain $D_j$ is a linear combination of the probabilities of samples belonging to all other classes having measurements in domain $D_j$.
This occurs for a choice of partition where all points fall in a single domain $D_j$. In this extreme case, there is an apparent dependence (in the linear algebra sense) of the measurement values of different classes.
As a related example, for the 1D SARS-CoV-2 antibody data from \cite{ainsworth2020performance} and \cite{wei2021antibody}, we can construct a partition where one trial domain is empty, both $\bm{P}$ and $\bm{P} - \bm{P_n}$ are singular, and therefore prevalence estimation is not possible.

To avoid this situation, one should select nonempty trial domains, i.e., training data should lie in each element of the partition.
In the limiting case that $P_{ij} = 0$, the measurement of a sample in class $C_j$ has zero probability of falling in domain $D_i$. The most extreme separation of training data occurs when the PDFs have nonzero support only on mutually exclusive elements of the partition. In this setting, the matrix $\bm{P}$ is a permutation matrix, and the prevalence estimates are merely the relative fractions $\bm{\hat{Q}}$ of measurements in each domain. If the partition elements are correctly matched to the classes, this extreme separation corresponds to a perfect assay because there are no misclassifications.
We note that the matrices that result from a $k$-means partition of the 1D SARS-CoV-2 tri-class data are close to permutation matrices; one example is
\begin{equation}
\bm{P} = \left[\begin{array}{ccc}
0.9749 & 0.0058 & 0.1055 \\
0.0237 & 0.0495 & 0.8101 \\
0.0015 & 0.9447 & 0.0844
\end{array} \right].
\label{eq:ex_P}
\end{equation}
Selection of trial domains with a high degree of class separation may be a key to our low-error prevalence estimates.

We speculate that under certain conditions a matrix $\bm{P}$ that is a permutation matrix may be optimal in the sense that it minimizes the prevalence estimate error.
In 1D, it may be possible to construct an optimization in terms of the samples assigned to each element of the partition, such as
\begin{equation}
\argmin_{f_1(x), f_2(x), \bm{P_{\pi}}} ||\bm{P_{\pi} P} - \bm{I} ||_2^2,
\end{equation}
where $f_1(x)$ and $f_2(x)$ are indicator functions determining which samples are assigned to elements 1 and 2 of the partition (without loss of generality).
Here, $\bm{P_{\pi} P}$ is the row permutation of $\bm{P}$ closest to the identity matrix, $\bm{I}$. For the matrix $\bm{P}$ given by (\ref{eq:ex_P}), for example, $\bm{P_{\pi}} = [\bm{e}_1; \bm{e}_3; \bm{e}_2]$, where $\bm{e}_j$ is the $j$th standard basis vector.

We leave a search for the minimum prevalence error estimate to future work; see \cite{patrone2022minimizing} for an approach to the binary case. The extension of their work to the multiclass setting is not obvious because the objective function to minimize can be generalized in many different ways.

As a final note on extreme cases of prevalence estimation, in expectation, the problem is unconstrained. The constrained problem may be a viable alternative when it is known that one prevalence is close to zero.

\subsection{Local accuracy}
\label{sec:loc_acc}

Recall that in Section \ref{sec:opt_class} we needed to consider sets of measurements with equal probability of belonging to more than one class. This concept is related to \textit{local accuracy}, $Z$, which compares the probability that a test sample belongs to a particular class and has measurement $\bm{r}$ to the measurement density of a test sample with measurement $\bm{r}$.
We generalize the binary version from \cite{patrone2022holdout} to the multiclass setting:
\begin{equation}
Z(\bm{r}, D_1, \ldots, D_n) = \frac{q_k P_k(\bm{r})}{Q(\bm{r})} = \frac{q_k P_k(\bm{r})}{\sum_{j = 1}^n q_j P_j(\bm{r})}, \quad \bm{r} \in D_k,
\end{equation}
where $\{D_k\}$ partitions $\Omega$. Let $Z^{\star}(\bm{r}) = Z(\bm{r}, D_1^{\star}, \ldots, D_n^{\star})$ be the local accuracy of the optimal solution to the multiclass problem.  It is straightforward to show that $1/n \leq Z^{\star} \leq 1$. Due to optimality, if $\bm{r} \in D_k^{\star}$, we have $ q_k P_k(\bm{r}) \geq q_j P_j(\bm{r})$ for $j \neq k$.
Then
\begin{equation}
n q_k P_k (\bm{r})  \geq \sum_{j = 1}^n q_j P_j(\bm{r}) = Q(\bm{r}),
\end{equation}
and so $q_k P_k(\bm{r})/Q(\bm{r}) =  Z^{\star}(\bm{r}) \geq 1/n$ for $r \in D_k^{\star}$. $Z^{\star}$ is maximized at 1 when $Q(\bm{r}) = q_k P_k(\bm{r})$ for $\bm{r} \in D_k^{\star}$.
 In the multiclass setting, we have $Z^{\star} = 1/n$ when
\begin{equation}
q_1 P_1(\bm{r}) = \ldots = q_n P_n(\bm{r}).
\end{equation}
We will refer to such an $\bm{r}$, if it exists, as a multipoint of the optimal domains. The lower bound on $Z^{\star}$ is only attained at a multipoint. To see this, consider some measurement $\bm{v}$ that is not a multipoint.  Then there exist $j, k \in \{1, \ldots, n\}$, $j \neq k$, such that $q_j P_j(\bm{v}) < q_k P_k(\bm{v})$. Then, since the classification is optimal, $\bm{v} \not \in D_j^{\star}$ and $\bm{v} \in D_m^{\star}$ for some $m$ (it may be that $m = k$). Clearly, $q_j P_j(\bm{v}) < q_m P_m(\bm{v})$. Further, $q_{\ell} P_{\ell} (\bm{v}) \leq q_m P_m(\bm{v})$ for $\ell \neq j$. It follows that
\begin{equation}
Q(\bm{v}) = \left[ \sum_{\substack{i = 1 \\ i \neq j}}^n q_i P_i(\bm{v}) \right] +  q_j P_j(\bm{v}) < (n-1) q_m P_m(\bm{v}) + q_m P_m (\bm{v})
\end{equation}
and so $Q(\bm{v}) < n q_m P_m(\bm{v})$, which gives $1/n < Z^{\star}(\bm{v})$ for a non-multipoint.

 The concept of local accuracy could be used to  decide which values to hold out in an indeterminate class in order to meet a global accuracy target \citep[see][]{patrone2022holdout}.
For any measurement, we can compute what the local accuracy would be if we chose to assign it to each class in turn. Using the SARS-CoV-2 tri-class example as an illustration, conducting this procedure on a measurement $\bm{r}$ may result in, say, similarly high local accuracies for the previously infected and vaccinated classes, but a low local accuracy for the na{\"i}ve class. In this situation and without an optimal classification scheme, we may not feel confident labeling the sample as previously infected or as vaccinated, since their probabilities are similar, but we can say the sample is almost certainly not na{\"i}ve. This leads naturally to the observation that any subset of classes can be combined to make a new class. In particular, we can reduce the problem to a binary classifier. For our serological example, perhaps it is desirable to consider previously infected and vaccinated samples together, or equally possible that it is difficult to tell them apart for a particular assay, and so our goal becomes to classify them separately from na{\"i}ves. This reduction of the problem size by combining classes is in a sense a projection onto a lower class space. Specifically, consider
\begin{equation}
Q(\bm{r}) = \sum_{j = 1}^n q_j P_j(\bm{r}) = \underbrace{\sum_{j = 1}^k q_j P_j (\bm{r})}_{\tilde{q}_1 \tilde{P}_1(\bm{r})} + \underbrace{\sum_{j = k+1}^n q_j P_j (\bm{r})}_{\tilde{q}_2 \tilde{P}_2(\bm{r})},
\end{equation}
where $\tilde{P}_1(\bm{r})$ and $ \tilde{P}_2(\bm{r})$ are newly-created PDFs with associated prevalences $\tilde{q}_1$ and $\tilde{q}_2 = 1 - \tilde{q}_1$. The task becomes to find $\tilde{q}_1$, for which there may be an optimal strategy.

The diagnostic community's current analog to local accuracy is the concept of a likelihood ratio (LR), calculated as $s_e/(1 - s_p)$ for a previously infected test result, where $s_e$ and $s_p$ represent sensitivity and specificity. The previously infected LR can be interpreted as the ratio of probabilities of correctly to incorrectly predicting a previously infected result \citep{riffenburgh2011statistics}. These values use average population information through $s_e$ and $s_p$ values, which may not always be available or representative. In contrast, local accuracy uses local information, since the latter is conditioned on knowing individual measurement values.

\subsection{Extensions}

Multiclass methods are readily equipped to handle further stratification of antibody data, such as by age group, biological sex, or coronavirus disease of 2019 (COVID-19) booster status. An additional class could be added for individuals who are both vaccinated and previously infected. Studies have demonstrated a greater antibody response post-vaccination for previously-infected versus COVID-19 na\"{i}ve recipients \citep{dalle2021serology,narowski2022sars}; this could allow for these populations to be distinguished by our classification scheme. Further, we minimize the prevalence-weighted combination of misclassifications, but the optimization problem can be rewritten for any desired objective function. Reformulations include ``rule-in'' or ``rule-out'' tests that meet desired sensitivity or specificity targets \citep{florkowski2008sensitivity}. Our methods may even be generalizable to multi-label classification, in which a sample can be assigned to more than one class; we anticipate challenges designing the corresponding optimization problem. Finally, the methods presented here can be applied to any setting where class size estimation and population labeling are required; an example is cell sorting in flow cytometry.

\subsection{Limitations}

Model selection is inherently subjective;  \cite{schwartz1967estimation} showed that the error goes to zero as more data points are added. As the number of antibody measurements increases, corresponding to viewing the data in higher dimensions, additional modeling choices become available.
 \cite{patrone2021classification} suggest the possibility of minimizing misclassifications over a family of models; see also \cite{smith2013uncertainty} for a discussion of model form errors. Classification accuracy and prevalence estimation of the 1D data sets from \cite{ainsworth2020performance} and \cite{wei2021antibody} suffer from overlap between their spike IgG values. If more measurements were available per sample, modeling the data in a higher dimension could improve class separation and thereby lower error rates \citep[see][]{luke2022improving}. Further, our models do not account for time-dependence. This concept is important when classifying antibody tests, which are known to have a half life on the order of several months post infection or vaccination \citep{xia2021longitudinal,kwok2022waning}. See \cite{bedekar2022prevalence} for a time-dependent approach to the binary setting.

\subsection{Implications for assay developers}

We have solved the multiclass diagnostic classification problem, which was previously unresolved. Antibody measurements from vaccinated individuals can now be distinguished from previously infected and na{\"i}ve samples.

Our work is the first to obtain unbiased predictions of the relative fractions of vaccinated, previously infected, and na{\"i}ve individuals in a population. These estimates are improved as more samples are added.  Best practices for conducting these predictions include dividing the range of all possible measurement values into nonempty regions that create separation between samples of neighboring regions. This can be easily achieved using pre-defined clustering algorithms. Our procedure hinges on selecting probability distributions to model training populations, which can be conducted automatically for measurements of a single antibody target in several open-source programming languages. Our classification scheme is also easily implementable, and can be modified to prioritize specificity if desired. Regardless of the reformulation, the error is minimized by construction.

\section{Acknowledgements}

This work is a contribution of the National Institute of Standards and Technology and is not subject to copyright in the United States. R.L. was funded through the NIST PREP grant 70NANB18H162. The aforementioned funder had no role in study design, data analysis, decision to publish, or preparation of the manuscript. Use of data provided in this paper has been
approved by the NIST Research Protections Office (IRB no. ITL 2020 2057). The authors wish to thank Drs. Daniel Anderson and Eric Shirley for useful discussions during preparation of this manuscript.

\section{Data Availability}

Analysis scripts and data developed as a part of this work are available upon reasonable request. Original data are provided in \cite{ainsworth2020performance} and \cite{wei2021antibody}.

\section{Declarations of Competing Interests}

The authors have no competing interests to declare.

\setcounter{equation}{0}


\section{Introduction}
\label{sec:intro}

Recent advances in colloidal synthesis and microscopy techniques have dramatically improved our ability to characterize how particle ordering and relaxation in thermal liquids varies with particle shape.
For example, high-quality monodisperse colloidal rods and cylinders can now be produced \cite{dugyala13}, and their liquid-state dynamics (in suspensions) can be observed using confocal microscopy and allied techniques \cite{weeks00,gasser01}.
Of particular interest is the recent experimental observation of {``}liquid glass{''}.
Roller \textit{et al.}\ found \cite{roller21} that suspensions of ellipsoidal colloids with aspect ratio $\alpha_{\rm R} = 3.5$ exhibited two distinct glass transitions at packing fractions $\phi_{\rm g}^{\rm rot}$ and $\phi_{\rm g}^{\rm trans}$.
In the liquid glass state ($\phi_{\rm g}^{\rm rot}  <  \phi < \phi_{\rm g}^{\rm trans}$), particles rotations' are arrested, but they remain free to translate within  locally-nematic precursor domains.
This phenomenon had been predicted by simulations \cite{letz00}, but had previously only been experimentally observed in quasi-2D systems \cite{zheng11,zheng14,mishra13}.
Experiments like Ref.\ \cite{roller21} offer both a new avenue for understanding the physics of anisotropic molecular glassformers and an obvious motivation for theoretical studies of related models.

For hard ellipses and ellipsoids of revolution with aspect ratio $\alpha$, complex liquid-state dynamics are expected for packing fractions in the range $\phi_{\rm {o}}(\alpha) < \phi < \phi_{\rm g}(\alpha)$, where $\phi_{\rm {o}}(\alpha)$ is the ``onset'' density \cite{sastry98,chaudhuri10} and $\phi_g(\alpha)$ is either the rotational or translational glass transition density.
Evaluating $\phi_{\rm {o}}(\alpha)$ and $\phi_{\rm g}(\alpha)$ using simulations is very { computationally expensive} \cite{pfleiderer08,shen12,davatolhagh12,xu15}.
A more easily obtained yet physically relevant packing fraction that lies between $\phi_{\rm {o}}(\alpha)$ and $\phi_g(\alpha)$ is the random sequential adsorption (RSA) density $\phi_{\rm s}(\alpha)$, the maximum density at which impenetrable ellipses of aspect ratio $\alpha$ can be packed under a protocol that sequentially inserts them with random positions and orientations.
The differences $\phi_{\rm g}(\alpha) - \phi_{\rm s}(\alpha)$ and $\phi_{\rm J}(\alpha) - \phi_{\rm s}(\alpha)$ are of particular interest because they indicate how much packing efficiency can be gained (before glass formation and jamming, respectively) by allowing particles to move freely while remaining positionally and orientationally disordered.
For example, in the $\alpha \to \infty$ limit one expects $\phi_{\rm J}(\alpha) = \phi_{\rm s}(\alpha)$ because particle rotations are completely blocked \cite{philipse96,desmond06}.

While jamming of ellipses and ellipsoids is now fairly well understood \cite{donev04,donev07,delaney05,schreck09}, RSA of these systems remains relatively poorly characterized.
$\phi_{\rm s}(\alpha)$ values for ellipses have been reported only for $1 \leq \alpha \leq 5$ \cite{sherwood90,viot92,ciesla16,haiduk18}, and detailed characterization of RSA ellipse packings' structure has only been performed for the aspect ratio $\alpha \simeq 1.85$ that maximizes $\phi_{\rm s}$ \cite{ciesla16}.
Thus there is a need to substantially expand our knowledge of these packings.

In this paper, we characterize saturated RSA ellipse packings over a wider range of aspect ratios and in much greater detail than has been previously attempted.
First we determine their packing fractions $\phi_{\rm s}(\alpha)$ to within $\sim 0.1\%$ for $1 \leq \alpha \leq 10$.
Then we characterize their positional-orientational order using several metrics.
We find a previously-unreported structural transition at $\alpha = \alpha_{\rm TS} \simeq 2.4$.
For $1 < \alpha < \alpha_{\rm TS}$ ($\alpha > \alpha_{\rm TS}$), packings have an excess of tip-to-side (side-to-side) contacts.
The peak prevalence of the favored type of contact, an effective order parameter for these systems, increases exponentially with $|\alpha - \alpha_{\rm TS}|$.
We also show that (i) local nematic order and structures comparable to the precursor domains observed in experiments \cite{roller21,zheng11} gradually emerge as $\alpha$ increases beyond $3$, and (ii) the increasing size of single-layer lamellae that are randomly oriented within these packings makes their long-wavelength density fluctuations increase rapidly with $\alpha$ for $\alpha \gtrsim 5$.

\section{Generating Saturated Packings}

Saturated RSA packings of anistropic 2D particles are generated by placing them with random positions and orientations, typically within square domains of size $L \times L$, until no more particles can be inserted.
In practice,  RSA packing generation's inherently slow kinetics \cite{swendsen81} make achieving complete saturation for $L$ that are large enough to minimize finite-size effects extremely difficult.
Therefore we employ an efficient protocol that generates packings that are demonstrably within $< 1\%$ (and for most $\alpha$, much less than $1\%$) of their saturation densities, for system sizes that are in the $L \to \infty$ limit.
Specifically, we use $L = 1000{\sigma}${, where $\sigma$ is the length of the ellipses' minor axes.
This} is the same {$L$} employed in Refs.\ \cite{ciesla16,haiduk18} and is sufficiently large that finite-size errors on $\phi_{\rm s}$ should be less than $0.1\%$ \cite{ciesla16}.
{ Below, we set $\sigma = 1$.}

We divide our periodic domains into $\Lambda = \rm{floor}(L/\alpha)^2$ linked cells of size $(L/\Lambda) \times (L/\Lambda)$.
Here $\rm{floor}[x]$ rounds $x$ downward to the nearest integer, e.g.\ $\rm{floor}[4.6] = 4$.
Thus, when insertion of an ellipse $i$ at position $\vec{r}_i$ is attempted, only ellipses in the cell containing $\vec{r}_i$
and its 8 neighboring cells need be checked for overlap.
Checks for overlap with the set $\{ j \}$ of all ellipses that have already been inserted in these cells (at positions $\{ \vec{r}_j \}$) include three steps:
(1) Since we assume ellipses' minor axes have unit length while their major axes have length $\alpha$, any center-to-center distances $r_{ij} = |\vec{r}_j - \vec{r}_i| < 1$ imply overlap, and the insertion attempt is rejected.
(2) If $r_{ij} > \alpha$, the ellipses do not overlap and the code continues to the next $j$.
(3) If $1 \leq r_{ij} \leq \alpha$, overlap is possible.
 We determine whether the ellipses overlap using  Zheng and Palffy-Muhoray's exact expression \cite{zheng07} for their distance of closest approach $d_{\rm cap}$.
{ For monodisperse ellipses with unit-length minor axes, their expression reduces to
\begin{equation}
d_{cap} = \frac{ d' }{\sqrt{ 1 -  (1 - \alpha^{-2})(\hat{k_i}\cdot\hat{r}_{ij})^2 } },
\label{eq:dcap}
\end{equation}
where $d'$ is obtained using a complicated formula involving $r_{ij}$ and the ellipses' orientation vectors $\hat{k}_i$ and $\hat{k}_j$.
Note that $\hat{k}_j$ does not appear in Eq.\ \ref{eq:dcap} because it has been absorbed into $d'$; cf.\ Eqs.\ 33, 36 of Ref.\ \cite{zheng07}. }
If $d_{cap} > r_{ij}$, the ellipses overlap and the insertion attempt is rejected.
Otherwise the ellipse is inserted and the next insertion attempt begins.

We examined 81 different particle aspect ratios over the range $1 \leq \alpha \leq 10$.
To allow extrapolation of the runs' progress to the infinite-time limit, we tracked the packing fractions $\phi(\alpha,t)$ after $t$ insertion attempts per unit area had been completed.
Each run continued until $2.5\times10^5$ trials per unit area had been attempted; this was sufficient to reach well into the asymptotic $t^{-1/3}$ regime \cite{swendsen81} for all $\alpha$.
Then we determined the RSA densities $\phi_{\rm s}(\alpha)$ by fitting the { results to}
\begin{equation}
\phi(\alpha,t) = \phi_{\rm s}(\alpha) \left(1 - \left[ \frac{t}{\tau(\alpha)}\right]^{-1/3} \right) {,}
\label{eq:t13}
\end{equation}
{ where $\tau(\alpha)$ is a ``time'' characterizing the $\alpha$-dependent RSA kinetics.}

\begin{figure}[h]
\includegraphics[width=3in]{fig1.eps}
\caption{RSA kinetics for ellipses.  Results are shown for  $\alpha = 1.3$, $1.85$, $2.4$, $3.5$, $5$, $6.5$, $8$, and $10$.}
\label{fig:phivst}
\end{figure}

Figure \ref{fig:phivst} shows {$\phi(\alpha,t)$}
for selected $\alpha$.
There is some crossing of the curves at intermediate $t$ for $\alpha < 2$ because $\phi_{\rm s}(\alpha)$ is non-monotonic, but as expected for anisotropic particles \cite{vigil89}, $\tau$ increases with $\alpha$.
All $\alpha \leq 10$ had $\tau < 0.13$ and hence had $\phi_{\rm s}(\alpha)/\phi(\alpha, 2.5\times10^5) - 1 < .0082$.
In other words, at {the end of our packing-generation runs,}
all $\alpha \leq 10$ have $\phi$ that are within less than $1\%$ of $\phi_{\rm s}(\alpha)$.

\section{Structure of Saturated Packings}

Figure \ref{fig:phis} shows the extrapolated $\phi_{\rm s}(\alpha)$.
Results for all $\alpha \leq 5$ agree with previous studies.
In particular, they are consistent with the well-known disk RSA packing fraction $\phi_{\rm s}(1) = \phi_{\rm disks} = .54707$ \cite{feder80,zhang13} and the known density maximum at $\alpha \simeq 1.85$, i.e.\ $\phi_{\rm s}(1.85) \simeq .584$ \cite{sherwood90,viot92,ciesla16,haiduk18}.
For large aspect ratios, our results show a surprisingly slow decrease of $\phi_s$ with increasing $\alpha$.
Specifically, since  the area swept out by ellipses as they rotate about their centers is $A_{\rm sw} = \pi\alpha^2/4$ whereas the area of the ellipses themselves is $\pi \alpha/4$, Onsager's classical arguments \cite{onsager49} predict $\phi_{\rm s} \sim A/A_{\rm sw} = 1/\alpha$.
Our data show that this asymptotic regime is not reached until $\alpha \gg 10$.

\begin{figure}[h!]
\includegraphics[width=3.25in]{fig2.eps}
\caption{Packing fractions of saturated {RSA} ellipse packings..
Red circles show the $\phi_{\rm s}(\alpha)$ estimated from Eq.\ \ref{eq:t13}, while the blue curve shows Eq.\ \ref{eq:analyticphis}.  The dashed gray line shows $\phi = 0.775\alpha^{-1/5}$ and is included only to illustrate the lack of convergence to $\alpha^{-1}$ scaling.}
\label{fig:phis}
\end{figure}

We find that  the extrapolated $\phi_{\rm s}(\alpha)$ are well fit by
\begin{equation}
\phi_{\rm s}(\alpha) = \phi_{\rm disks} \times \displaystyle\frac{1 + (3/8)\ln(\alpha) +  (17/25)(\alpha-1)}{1 + (80/99)(\alpha - 1) + (\alpha-1)^2/96}.
\label{eq:analyticphis}
\end{equation}
$\phi_{\rm s}(\alpha) > \phi_{\rm disks}$ for all $1 < \alpha \lesssim 4.40$, indicating that particles' anisotropy enhances their packability over this range of aspect ratios.
Equation \ref{eq:analyticphis} is an empirical formula designed to capture the known trends: convergence to $\phi_{\rm disks}$ as $\alpha \to 1$, and $1/\alpha$ scaling at large $\alpha$ {\cite{philipse96,desmond06}}.
Its mean and rms fractional deviations from the extrapolated $\phi_{\rm s}(\alpha)$ are respectively $0.02\%$ and $0.06\%$.
We do \textit{not} claim that Eq.\ \ref{eq:analyticphis} is an exact expression valid for all $\alpha$,
{ or even that its functional form is the same as that of the ``true'' $\phi_s(\alpha)$ which could be obtained given infinite computer power.}

Next we turn to detailed analyses of these packings' positional and orientational order.
Figure \ref{fig:gands} shows the pair correlation functions $g(r)$ and structure factors $S(q)$ for the same eight aspect ratios highlighted in Fig.\ \ref{fig:phivst}.
Three major trends are apparent.

First, all results for $\alpha \leq 2.2$ are consistent with previous studies \cite{ciesla16,haiduk18}.
The $g(r)$ exhibit a primary nearest-neighbor peak at $r_{\rm nn} \simeq (\alpha + 1)/2$ and a much smaller second-nearest-neighbor peak at $r_{\rm sn} = 2 +  \mathcal{O}(\alpha)$.
Both peaks rapidly broaden as $\alpha$ increases.
The primary peak's position indicates that it corresponds to tip-side contacts.
A similar peak is present in RSA packings of low-$\alpha$ spherocylinders; this feature is likely universal for packings of convex anisotropic objects \cite{ricci94,ciesla16}.
The peak at $r_{\rm sn}$ vanishes for $\alpha \gtrsim 1.9$; ellipses with $\alpha \gtrsim 1.9$'s failure to form  an approximately-isotropic second coordination shell during RSA helps explain why their $\phi_s$ is maximized at this $\alpha$.

Second, { for all $4 \lesssim \alpha \lesssim 8$,} $g(r) < 1$ for all $r < (\alpha + 1)/2$.
This ``correlation hole'' \cite{degennes79} can be understood in terms of two-body impenetrability constraints.
Suppose $\Delta \theta$ is the difference in the orientation angles of two ellipses.
If the ellipses are in contact, the maximal distance $r_{\rm max}$ between their centers ranges from $1$ for $\Delta\theta = 0$, which corresponds to side-to-side contacts of perfectly aligned ellipses, to $(\alpha + 1)/2$ for $\Delta\theta = 90^\circ$, which corresponds to tip-to-side contacts of perpendicularly oriented ellipses.
Therefore, for any pair of ellipses (contacting or not) with $r < (\alpha + 1)/2$, the allowed $\Delta \theta$ are restricted, increasingly so as $r$ decreases towards $1$.
This tends to reduce the number of ellipse pairs with center-to-center distances $r < (\alpha + 1)/2$ relative to the number that would be present in an ideal gas.
The correlation hole emerges when this effect overwhelms the natural tendency of particles in RSA packings to form coordination shells \cite{swendsen81}.
{ It disappears when a side-to-side contact peak at $r \simeq r_{\rm ss}$ [defined by $g(r_{\rm ss}) > 1$] emerges as $\alpha$ increases beyond $8$.}

\begin{figure}[h!]
\includegraphics[width=3in]{fig3.eps}
\caption{Pair correlation functions $g(r)$ [panel (a)] and structure factors $S(q)$ [panel (b)] for the same systems highlighted in Figure \ref{fig:phivst} { \cite{footSofQvec}}.  Dotted curves show results for disks.  The inset to panel (b) shows estimated $S_0 = \lim_{q \to 0} S(q)$, here defined as $S(0.25)$.}
\label{fig:gands}
\end{figure}

Third, examining the $S(q)$ curves reveals how density fluctuations evolve with increasing $\alpha$.
Low-$\alpha$ systems' $S(q)$ are qualitatively similar to disks' $S(q)$.
They have prominent peaks at $q \simeq 2\pi/r_{\rm nn}$ and $q \simeq 4\pi/r_{\rm nn}$, which respectively correspond to the first- and second-nearest neighbor peaks shown in Fig.\ \ref{fig:gands}(a).
Two features which are not readily apparent from Fig.\ \ref{fig:gands}(a) are:\ (i) intermediate-wavelength ($1 \lesssim q \lesssim 4$) density fluctuations increase monotonically with $\alpha$, and (ii) longer-wavelength density fluctuations increase rapidly -- as shown in the inset, roughly exponentially -- for $\alpha \gtrsim 5$.
We will show below that the second trend corresponds to the emergence of lamellae that grow larger with increasing $\alpha$.

\begin{figure*}[htbp]
\includegraphics[width=7in]{fig4-v2.pdf}
\caption{Representative snapshots and $g(r,\Delta\theta)$ for $\alpha = 1.3$, $1.85$, $2.4$, and $3.5$ (top rows) and $\alpha = 5$, $6.5$, $8$, and $10$ (bottom rows) from left to right.  The snapshots show only ellipses whose \textit{centers} lie within the central $50\times 50$ regions of the overall $L\times L \equiv 1000\times 1000$ domains.}
\label{fig:snaps}
\end{figure*}

The abovementioned trends can be better understood by visualizing saturated RSA packings and their positional-orientational order.
Figure \ref{fig:snaps} shows snapshots of typical packings for the eight $\alpha$ highlighted above, along with results for the pair correlation function $g(r,\Delta\theta)$, which is the ratio of the number of ellipse pairs with center-to-center distance $r$ and orientation-angle difference $\Delta\theta$ to the number that would be present in an ideal gas of these particles.
In other words $g(r,\Delta\theta)$ is just the generalization of the standard pair correlation function $g(r)$ to include orientation-angle differences { \cite{footSofQvec}}.

Results for $\alpha \lesssim 2$ are consistent with previous studies \cite{sherwood90,viot92,ciesla16,haiduk18}, but shed new light on low-aspect-ratio RSA ellipse packings' structure because these systems'  $g(r,\Delta\theta)$ has not been previously reported.
The snapshots show the expected excess of tip-side contacts, and the $g(r,\Delta\theta)$ plots illustrate the degree to which this favored ordering diminishes { as} $r$ and $\Delta\theta$ vary away from the ideal tip-side configuration [$r_{\rm its} = (\alpha+1)/2,\ \Delta\theta = 90^\circ$].
The primary peaks in $g(r,\Delta\theta)$ occur at distances $r_{\rm ts} = r_{\rm its} + \mathcal{O}(\alpha)$ that increase with $\alpha$, but remain near $\Delta\theta = 90^\circ$.
Secondary peaks at distances $r_{\rm 2}(\alpha, \Delta\theta)$ that increase with increasing $\alpha$ and decreasing $\Delta\theta$, roughly mirroring the small-$r$ excluded-area cutoffs $r_{\rm min}(\alpha,\Delta\theta)$ below which $g(r,\Delta\theta) = 0$, are also evident.
These reflect the fact that there are at least two distinct favored near-neighbor configurations for any given $\Delta\theta$.

The snapshot and $g(r,\Delta\theta)$-plot for $\alpha = 3.5$ illustrate how the structures of saturated ellipse packings with $\alpha = 3-4$ differ from their lower-aspect-ratio counterparts.
Local nematic precursors comparable to those observed in both previous Monte Carlo simulations of hard ellipses and spherocylinders \cite{ricci94,xu13b} and experiments on quasi-2D colloidal-ellipsoid suspensions  \cite{zheng11,zheng14,zheng21,tan21} are readily visible.
Such precursor domains are believed to be largely responsible for these systems unusual liquid-state dynamics, and in particular for their ability to exhibit separate rotational and translational glass transitions \cite{wang18b,zheng11,zheng14,zheng21,roller21}.

The bottom rows of Fig.\ \ref{fig:snaps} show that for $\alpha \gtrsim 5$, saturated RSA ellipse packings possess a common structure that evolves quantitatively but not qualitatively with increasing $\alpha$.
These packings' $g(r,\Delta\theta)$ are simpler, with little structure aside from their side-to-side-contact peaks near $(r_{\rm ss},0)$.
The snapshots indicate that packings are largely composed of randomly-oriented, single-layer lamellae.
In contrast with the longer-range nematic order that develops in { sufficiently} dense systems of ellipses and other anisotropic particles in thermodynamic equilibrium \cite{ricci94,xu13b,bautista14b}, these lamellae appear have a kinetic origin; see below for a potential explanation.
The combination of random lamellar orientation and increasing lamellar width (i.e.\ increasing particle length) with increasing $\alpha$ directly explains these packings' rapidly increasing long-wavelength density fluctuations [Fig.\ \ref{fig:gands}(b)].
Local particle density is high inside lamellae and low where lamellae of different orientations intersect (or, more precisely, where they would intersect if they were able to grow further).
Since (i) the size { of the lamellae} and (ii) the density disparity between the lamellae and their intersection regions both increase rapidly with $\alpha$, so does $S_0$.
The same effect is seen in RSA of zero-thickness line segments \cite{ziff90}.

One natural order parameter for these packings is $g_{\rm max} = \max[g(r,\Delta\theta)]$, the maximal excess positional-orientational correlation, here defined to occur for pairs with  center-to-center distance $r_{\rm max}$ and orientation-angle difference $\Delta\theta_{\rm max}$.
Comparing the snapshots and $g(r,\Delta\theta)$ shown in the top rows of Fig.\ \ref{fig:snaps} highlights our central result:\ the existence of a previously-unreported transition in RSA ellipse packings' structure.
 $g_{\rm max}$ decreases rapidly with increasing $\alpha$ over the range $1 < \alpha \leq \alpha_{\rm TS} \simeq 2.4$ as the peaks in $g(r)$ and $g(r,\Delta\theta)$ [respectively at $r_{\rm nn}$ and $r_{\rm ts}$] broaden and shift outward.
As $\alpha$ increases beyond $\alpha_{\rm TS}$, the positions of the peaks in $g(r,\Delta\theta)$ change discontinuously.
Specifically, $r_{\rm max}$ decreases from $r_{\rm ts} \simeq (\alpha+1)/2$ to $r_{\rm ss} \simeq 1.1$, and $\Delta\theta_{\rm max}$ decreases from just below $90^\circ$ to just above $0^\circ$.
Although the heights of the peaks in  $g(r,\Delta\theta)$ at $r_{\rm ts}$ and $r_{\rm ss}$ both evolve continuously with increasing $\alpha$, this discontinuous change in ($r_{\rm max}$, $\Delta\theta_{\rm max}$) can be regarded as a structural transition from tip-side- to side-side-dominated contact { \cite{footTStrans}}.

\begin{figure}[htbp]
\includegraphics[width=3in]{fig5.eps}
\caption{ Exponential increase in maximal positional-orientational correlations with increasing $|\alpha-\alpha_{\rm TS}|$.  Dotted lines show Eq.\ \ref{eq:fitg}.  The inset shows the orientation-angle differences $\Delta\theta_{\rm max}$ at which the maxima occur, highlighting the transition from tip-side- to side-side-dominated contact.}
\label{fig:gmax}
\end{figure}

For $\alpha > \alpha_{\rm TS}$, $g_{\rm max}$ increases rapidly with $\alpha$ as side-to-side nearest-neighbor contacts with $r \simeq r_{\rm ss}$ become increasingly prevalent.
We find that in fact $g_{\rm max}$ increases exponentially with $|\alpha - \alpha_{\rm TS}|$.
As shown in Figure \ref{fig:gmax},  $g_{\rm max}(\alpha)$ is well fit by
\begin{equation}
g_{\rm max} \simeq \Bigg{ \{ } \begin{array}{lcr}
g_0 \exp\left( \displaystyle\frac{\alpha_{\rm TS}^* - \alpha}{\alpha^*_{\rm low} } \right) & , & \alpha < \alpha_{\rm TS} \\
\\
g_0 \exp\left( \displaystyle\frac{\alpha_{\rm TS}^* - \alpha}{ \alpha^*_{\rm high} }  \right) & , & \alpha > \alpha_{\rm TS}
\end{array},
\label{eq:fitg}
\end{equation}
for all $\alpha \lesssim 5$, with $g_0 \simeq 2.70$, $\alpha_{\rm TS}^* \simeq 2.36$, $\alpha^*_{\rm low} \simeq 1.79$ and $\alpha^*_{\rm high} \simeq 2.90$.
Systems with $\alpha = \alpha_{\rm TS}^*$ can be considered maximally { locally} disordered.
They are frustrated in the sense that tip-side- and side-side-contact-dominated ordering are mutually incompatible { yet equally likely.}

While the values of the parameters $g_0$, $\alpha_{\rm TS}^*$, $\alpha^*_{\rm low}$ and $\alpha^*_{\rm high}$ depend on how the binning used to calculate $g(r,\Delta\theta)$ is performed, the functional form of Eq.\ \ref{eq:fitg} is robust.
Rigorously explaining this result is challenging.
In particular, two-body geometry alone does not explain it, because although the range of $\Delta\theta$ [$0 \leq \Delta\theta \leq \Delta\theta^*(r)$] over which a new ellipse can be inserted with a center-to-center distance $r \lesssim r_{\rm ss}$ from its nearest neighbor decreases with increasing $\alpha$, this decrease would only produce $g_{\rm max} \sim 90^\circ/\Delta\theta^*(r) \sim \alpha$.
We suspect that, at least for $\alpha > \alpha_{\rm TS}$, the exponential increase of $g_{\rm max}$ arises from a combination of geometry and kinetics.
The locally densest configurations of ellipses with $r \simeq r_{\rm ss}$ are single-layer lamellae.
As RSA proceeds, the gaps into which new particles may be inserted increasingly correspond to configurations that increase the length of such lamellae, i.e.\ grow the lamellae at their ends.
Thus $g_{\rm max}$ increases much faster than linearly with $\alpha$, falling below its initial exponential trend only when lamellar growth becomes limited by their intersection with other, differently oriented lamellae, i.e.\ at $\alpha \simeq 5$.
{ For $\alpha > 5$ the growth of $g_{\rm max}$ with $\alpha$ appears to be linear.}

\section{Discussion and Conclusions}

Hard ellipse liquids' dynamics exhibit multiple nontrivial features arising from the asymmetry of their constituent particles \cite{zheng11,zheng14,mishra13,zheng21,shen12,xu15,wang18b,shen12,davatolhagh12}.
Thus they provide a minimal, single-parameter model for understanding aspects of the  dynamics of small-molecule glassforming liquids that are tightly coupled to particle asymmetry, such as translational-rotational decoupling \cite{chong09}.
They also provide an easily-visualized system that can be used as a basis for understanding their more complicated three-dimensional-ellipsoidal counterparts \cite{mishra14,mishra15,roller21,roller20,pal20,pal21}.

In this paper, we characterized the structure of asymptotically saturated RSA ellipse packings, over a far wider range of aspect ratios and in substantially more detail than had been previously reported.
We obtained a simple analytic expression for the saturation density $\phi_s(\alpha)$ [Eq.\ \ref{eq:analyticphis}] that captures results for all $\alpha \leq 10$ to within less than $0.1\%$ and has the correct $1/\alpha$ scaling \cite{philipse96} in the large-$\alpha$ limit.
The crossover to this asymptotic scaling was slower than expected; while it must eventually occur since $\phi_s$ is necessarily below the jamming density $\phi_J(\alpha)$ \cite{donev04,donev07,delaney05,schreck09}, we showed that it does so only at some $\alpha \gg 10$.

We also characterized how positional-orientational correlations in these packings evolve with increasing $\alpha$.
Our main result for relatively low aspect rations ($\alpha \leq 5$) was finding a previously-unreported structural transition from tip/side- to side/side-contact-dominated structure at $\alpha = \alpha_{\rm TS} \simeq 2.4$.
At this aspect ratio, the peak value $g_{\rm max}$ of packings' positional-orientational pair correlation function $g(r,\Delta\theta)$ is minimal, and packings can be considered maximally frustrated owing to competition of these two, mutually incompatible types of structural order.
$g_{\rm max}$ increases exponentially with increasing $|\alpha - \alpha_{\rm TS}|$, presumably owing primarily to kinetic effects.

{ $\alpha \simeq 2.4$ is also the aspect ratio above which ellipses possess an equilibrium nematic liquid-crystalline phase \cite{bautista14b}, suggesting that the abovementioned structural transition is intimately connected to the equilibrium isotropic-nematic transition.
Explaining such a connection -- if indeed one exists -- will be nontrivial because saturated RSA packings of anistropic 2D particles often have local structure that is substantially different than that of their equilibrium counterparts at the same $\phi$ \cite{ricci94}.
One potential strategy for obtaining such an explanation is extending replica-trick-based techniques, which were very recently employed to successfully predict RSA sphere packings'  structure \cite{jadrich22}, to anisotropic particles.}

RSA ellipse packings with $\alpha > 5$ had not been previously investigated.
We found that these packings possess a common structure that varies quantitatively but not qualitatively with increasing $\alpha$:\ they are primarily composed of randomly oriented single-layer lamellae.
$g_{\rm max}$ falls below the abovementioned exponential trend, presumably because  lamellar growth is limited by ``collisions'' with other, differently oriented lamellae.
The same collisions drive a rapid increase in long-wavelength density fluctuations as the size of and density disparity between lamellae and their intersection regions increase with increasing $\alpha$.

In conclusion, the detailed characterization of saturated RSA ellipse packings performed in this work, in addition to being of interest in its own right, lays the groundwork for followup studies of the structure and dynamics of thermally-equilibrated, glassy and jammed ellipses.
For example,
it would be interesting to check whether a comparable structural transition and comparable behavior of $g_{\rm max}$ are present in jammed ellipse packings, and if so, at which aspect ratio it occurs.
Our results suggest that followup studies examining  $g(r,{\Delta}\theta)$ and/or systems with $\alpha > 5$ may prove especially fruitful.

This material is based upon work supported by the National Science Foundation under Grant No.\ DMR-2026271.