\section{Introduction and Related Work}
\vspace{-0.3cm}

Current medical diagnoses are most often based on the combination of several data inputs by medical experts, typically including (i) clinical history, interviews, and physical exams, (ii) laboratory tests, (iii) electrophysiological signals, and medical images. Advances in the machine learning (ML) community have highlighted the potential of ML to contribute to the field of Computer-Aided Diagnosis (CAD), for which we distinguish two main classes of methods.


\xhdr{Single modality analysis} Most recent efforts in the ML community have focused on analyzing a single medical data source --- called a ``modality". For instance, the parsing of electronic health records (EHR) for clinical history, patients' interviews, and physical exams has shown huge potential for the diagnosis of a broad range of diseases, from coronary artery disease to rheumatoid arthritis \cite{Harrell1984RegressionPrediction, Kurt2008ComparingDisease, Carroll2011NaiveArthritis., Hippisley-Cox2013PredictingScore, Rahimian2018PredictingRecords}. ML algorithms have been developed to process various types of laboratory results, from urine steroid profiles to gene expression data \cite{Huang2018NovelClassification, Wilkes2018UsingProfiles, Demirci2016ArtificialAlgorithms, Agrahari2018ApplicationsMalignancies}. Others have shown success in automatically processing and classifying medical signals such as electrocardiograms (ECG) \cite{Acharya2017AutomatedNetwork, Kiranyaz2016Real-TimeNetworks, Rad2017ECG-BasedAnalysis, Pawiak2018NovelSystem}, or electroencephalograms (EEG) \cite{Richhariya2018EEGMachine}. Finally, a large body of work has focused on medical imaging analysis encompassing a vast number of tasks, such as automatic extraction of diagnostic features \cite{Thomaz2017FeatureImages, Bar2018ChestTraining}, segmentation of anatomical structures \cite{Moeskops2016AutomaticNetwork, Li2015AutomaticNetworks, Shelhamer2017FullySegmentation, Milletari2016V-Net:Segmentation, Dolz20183DStudy, Christ1919AutomaticFields}, or direct diagnosis through image classification \cite{Anthimopoulos2016LungNetwork, Miki2017ClassificationNetwork}. 

While these works achieve record-breaking diagnostic accuracy, they often rely on supervised learning approaches -- requiring the diagnosis ground-truth to be available during training -- and on the acquisition of large datasets. Furthermore, they focus solely on a single type of data input (EEGs, scans, etc.) -- often acquired by clinicians using specialized, highly accurate equipment -- and do not harvest the potentially rich and complementary sources of information provided by alternative medical modalities. Yet, with the development of at-home diagnostic tests (lateral flow assays, questionnaire data for disease screening, mobile health apps, etc.), this paradigm shifts, and diagnoses have to be established through the combination of multiple sources of cheaper, yet often noisier and more imprecise data. 
In light of the uncertainty exhibited by these various inputs, it also becomes indispensable to pair the provision of a diagnosis with a notion of a confidence interval, thereby thoroughly characterizing our state of knowledge given the available data. 


\xhdr{Multiple modality integration} Interestingly, diagnosis studies fusing different data sources seemed to be more prominent in the first years of CAD. Bayesian networks \cite{Neal1993ProbabilisticMethods} -- also called belief networks -- have been used as a crucial decision tool for automatic diagnosis. Such networks provide a biologically-grounded and interpretable statistical framework that integrates the probabilities of contracting the disease given a patient's clinical history, and the probability of developing symptoms or observing positive test results in the presence or absence of the disease. The advantages of these Bayesian networks are that they are (a) fully unsupervised and do not require ground-truth information, and (b) able to provide meaningful results even in low sample regimes.

Bayesian networks have thus been implemented in a variety of contexts to integrate clinical data and laboratory results, and diagnose conditions ranging from pyloric stenosis to acute appendicitis \cite{Alvarez2006EvaluationStenosis, Gevaert2006PredictingNetworks, Sakai2007AccuracyModel, Gonzalez-Lopez2016DevelopmentUveitis,Sangamuang2018AccuracyDecision}. They have also been used to combine clinical data with medical images, and subsequently applied to assess venous thromboembolism \cite{Kline1984DerivationThromboembolism}, community-acquired pneumonia \cite{Aronsky1998DiagnosingNetwork.}, head injury prognostics \cite{Sakellaropoulos1999DevelopmentTechniques}, or to predict tumors \cite{Wang1999Computer-assistedNetwork,Kahn1997ConstructionCancer, Sneha2013TowardsImages}. As early as 1994, full integration of clinical data, laboratory results and imaging features was performed to diagnose gallbladder disease \cite{Haddawy1994ADisease}. 

However, contrary to our proposed setting, the uncertainty that these Bayesian networks integrate is typically known and controlled. The medical conditions that they study are well characterized by a set of specific questions and physiological exams (e.g. projectile vomiting, potassium levels and ultrasounds in the case of Pyloric stenosis). Not only do these inputs provide a very strong signal for the diagnosis, but the uncertainty arising from the different modalities can often be reliably informed by a test manufacturer, extensive medical research or prior ML studies (for diagnostic inputs obtained through classification algorithms). The uncertainty of a diagnostic method, however, is not always well specified: whether it be (i) a physician's assessment of the symptoms associated with a novel or rare disease, (ii) predictions of an ML algorithm whose accuracy has not been fully characterized outside of a curated research environment or reference datasets (such as MNIST, CIFAR), or simply (iii) biological tests whose sensitivity still has to be determined. 




\xhdr{Multiple noisy modality integration: a COVID-19 case study} To motivate this paper and further understand the limitations of these previous approaches, let us consider a particular use case: COVID-19 antibody testing. 
Lateral flow immunoassay (LFA) antibody tests for COVID-19 are one of the manageable, affordable and easily deployable options to allow at-home testing for large population and provide assessments of our immunity.
Yet, studies have shown that the sensitivity of these tests remains highly variable and highly contingent on the time of testing.
The successful deployment of LFAs thus depends on their augmentation with additional data inputs, such as user-specific risk factors and self-reported symptoms.
The provision of confidence scores is essential to flag potential false negative or positive tests (requiring re-testing or closer scrutiny) and to assess local prevalence levels -- both pivotal for researchers and policymakers in the context of a pandemic.

\xhdr{Contributions and outline} This applications paper is geared towards the practical integration of noisy sources of information for CAD. Our contribution is two-fold. From a methods perspective, we account for the variability of the inputs by devising a two-level Bayesian hierarchical model.
In contrast to existing Bayesian methods for CAD, our model is deeper, and trained using a Stochastic Expectation Maximization (EM) algorithm \cite{celeux1995stochastic,Celeux2001OnAlgorithm,nielsen2000stochastic}. The Stochastic EM allows to overcome the limitations of its non-stochastic counterpart \cite{celeux1995stochastic}, that is (a) its sensitivity to the starting position, (b) its potential slow convergence rate, and (c) its possible convergence to a saddle position instead of a local maximum. From an applications perspective, we gear this algorithm to enhance at-home LFA testing in the context of the COVID-19 pandemic. In particular, we wish to (a) quantify the benefit of multimodal data integration when the diagnostics are uncertain, and (b) show how our method can benefit medical experts or researchers in real life. 
Section~\ref{sec:multimodal} presents the Bayesian model for multimodal integration and the Stochastic Expectation-Maximization algorithm that performs principled and scalable inference. Section~\ref{sec:simu} presents extensive tests of our model on simulated datasets. Section~\ref{sec:real} details the results obtained on the Covid-19 dataset and the impact of our method for affordable and reliable at-home test kits.




\section{Covid-19 Dataset}\label{sec:data}

By way of clarifying the challenges that our algorithm proposes to overcome, we present here the COVID Clearance Remote Recovery \& Antibody Monitoring Platform study 
\footnote{COV-CLEAR, \url{www.cov-clear.com}}, which motivated our approach. The purpose of this study is to track the evolution of the immunity of a cohort of adult participants in the UK with various COVID-19 exposure risks. This paper focuses on a subset of 117 healthcare workers, for which we have both questionnaire information and LFA test results.


\xhdr{Home-testing Immunoassay Data} Participants are issued packs of home-testing kits with written instructions. These kits identify Immunoglobulin M (IgM) and/or Immunoglobulin G (IgG) specific to SARS-CoV2 in blood samples using three gold-standard methods \footnote{Home-sampled (blood) Abbott ARCHITECT® Anti-SARS-CoV-2 IgG CMIA CE Marked; PHE Approved} \footnote{Elecsys® Anti-SARS-CoV-2 Double Antigen, CE Marked, PHE Approved} \footnote{Biopanda Ltd  Product Number: RAPG-COV-019 version 1 In-Vitro LFT IgM IgG Antibody test. CE certified for Professional Use Only.}. Pictures and typical examples of LFA test results are provided in the supplementary materials. For the LFAs, participants self-report a positive result if IgM and/or IgG are detected by the test in addition to a positive control band. Through the collection of this data, the COV-CLEAR study aims to address questions relating to  (a) the quantification of the robustness of the antibody response, and (b) the durability of this response. However, while affordable, we highlight that the LFA tests suffer from vastly varying levels of sensitivity --- registering a sensitivity as low as 70\% on asymptomatic individuals.

\xhdr{Clinical Data: Questionnaire} Additionally, participants are asked to answer a questionnaire, ideally before knowing the result of their LFA test. The form consists of questions related to $K=14$ exhibited symptoms (fever, cough, runny nose, etc.) and $M=2$ subject-specific risk factors (household size and proportion of members with a suspected or confirmed history of COVID-19). The empirical distributions of the symptoms and the risk factors in this cohort are provided in the supplementary materials. 







\section{Multimodal Data Integration: Statistical Model}\label{sec:multimodal}
\vspace{-0.3cm}

In this section, we describe our principled integration of noisy diagnostic test results, with additional clinical data such as symptom data and subject-specific risk factors.
Our approach applies to general heterogeneous medical data where the outputs are binary. The latter could be either self-reported answers to questionnaires, clinician-reported physiological exams, outputs of a diagnostic based on an image, abnormalities of laboratory results, etc., making this a widespread and general setting. Thus, while we implement and showcase the method for the particular purpose of applying it to antibody testing, this could be relevant to any medical diagnostic with binary inputs.

Denote $D$ the true diagnosis of an individual (healthy/sick), $T$ the outcome of the noisy diagnostic test (positive/negative), $S$ the symptomaticity (symptomatic/asymptomatic), $X$ the symptoms exhibited and $Y$ the subject-specific risks factors. The underlying assumption is that given a true diagnosis $D$, the symptoms $X$ and the diagnostic test outcome $T$ are independent. In other words, the probability of the diagnostic test being a false negative $\mathbb{P}[T=0|D=1]$ is independent of the symptoms of a truly infected individual $\mathbb{P}[X|D=1]$. Similarly, given a diagnosis $D$, the test outcome $T$ and the exhibited symptoms $X$ are independent of the risk factors $Y$. We define:
\begin{itemize}[leftmargin=3em]
       \item $\mathbb{P}(D=1|Y) = \pi_\beta(Y)$ the probability of contracting the disease given risk factors $Y$,
       \item $\mathbb{P}[(T=1|D=1) = x$ the sensitivity of the diagnostic test,
       \item $\mathbb{P}[T=0|D=0] = 1-y$ the specificity of the diagnostic test, \textit{i.e.} $y=\mathbb{P}[T=1|D=0]$,
       \item $\mathbb{P}(S=1|D=0) = p_0$ the probability of being symptomatic when whilst not having been infected by that specific disease (the symptoms could be due to another illness for instance),
       \item $\mathbb{P}(S=1|D=1) = p_1$ the probability of having been symptomatic upon infection,
       \item $\mathbb{P}(X_k=1|S=1,D=0) = s_{0k}$ the probability of exhibiting symptom $k$ when not infected,
       \item $\mathbb{P}(X_k=1|S=1,D=1) = s_{1k}$ the probability of exhibiting symptom $k$ upon infection.
\end{itemize}
In the above, the diagnostic test may be any biological test (e.g. LFA antibody tests), or output of a medical imaging classification algorithm provided by a ML research team. We have here splitted our binary variables into two classes according to their variability: (a) the test $T$, for which we have provisional estimates for the sensitivity and specificity (given by the manufacturer) but which we are still uncertain, and (b) the symptoms $X$, which carry additional uncertainty in that neither of them is extremely specific to COVID-19 nor is their prevalence amongst COVID cohorts very well established; Values for $(x, y)$, $p_0, p_1, s_0, s_1$ or $\beta$ may be published, but challenged by complementary research or field experience, or may be completely unavailable.
We thus need to account for these inputs' variability and for the relative importance of different risk factors $\beta$. 
We propose using a hierarchical Bayesian network (see Fig.~\ref{fig:model}) that is consequently deeper than the ones traditionally implemented in the CAD literature \cite{Alvarez2006EvaluationStenosis, Gevaert2006PredictingNetworks, Sakai2007AccuracyModel, Gonzalez-Lopez2016DevelopmentUveitis,Sangamuang2018AccuracyDecision}.
 
The uncertainty in the test sensitivity $x$ and specificity $1-y$ are expressed through Beta priors on $x$ and $y$ with parameters $(\alpha_{x}, \beta_x)$ and $(\alpha_y, \beta_y)$, see Fig.~\ref{fig:model}. Similarly, the uncertainty on the probabilities of symptomaticity and of the different symptoms are expressed via Beta priors with respective parameters $(\alpha_{p_0}, \beta_{p_0}), (\alpha_{p_1}, \beta_{p_1}), \{(\alpha_{s_{0k}}, \beta_{s_{0k}}), (\alpha_{s_{1k}}, \beta_{s_{1k}})\}_k$, see Fig.~\ref{fig:model}. When published estimates exist for $x, y, p_0, p_1, s_0$ and $s_1$ (e.g. as provided by the LFA test manufacturer), we match the moments of the Beta priors with those of the published distributions.
If this information is unavailable, we use the non-informative Beta prior with parameters $(0.5, 0.5)$. In the COVID-19 example, we implement a non-informative prior for the probabilities of symptomaticity and symptoms, as the etiology of the disease remains unknown. These priors are then updated during learning as we aggregate the information from the observed and imputed variables.

Lastly, we model the probability $\pi_\beta(Y)$ of contracting the disease with a logistic regression on the $M$ risk factors $Y$. In the case of the COVID-19 dataset, these include county data, profession, size of household, etc. The coefficients $\beta=\{\beta_m\}_m$ weight the relative importance of each factor $Y_m$, $m=1...M$. We express our uncertainty on the possible impact of a given factor $m$ by introducing a Gaussian prior: $\beta \sim N(0, \sigma_\beta^2)$, see Fig.~\ref{fig:model}.



\begin{figure}[h!]
\centering
\def0.8\columnwidth{0.8\columnwidth}
\input{./total_model_with_legend.pdf_tex}
    \caption{Hierarchial bayesian network integrating the accuracy uncertainty on the data sources, to estimate the true diagnosis $D$. The index $k=1...K$ represents symptom $k$, while $m=1...M$ represents risk factor $m$. The shaded cells represent observed variables or known hyper-parameters. The dashed lines represent switches.}
    \label{fig:model}
\end{figure}
\vspace{-0.3cm}


\section{Inference in the Hierarchical Bayesian Network via Stochastic EM}
\vspace{-0.3cm}

\textit{At the subject level}, we seek to compute the posterior distribution of the true diagnosis $D_i$, informed by the integration of the observed variables $T_i, S_i, X_i, Y_i$. This posterior yields an estimated diagnosis, together with a credible interval that expresses the confidence associated with our prediction. In the case of COVID-19, this provides individual estimates of each citizen's immunity that may inform their social behavior as lockdown is eased.

\textit{At the global level}, we wish to learn: (i) the distributions of the sensitivity $x$ and specificity $y$ of the diagnostic test as observed in the field; (ii) the distributions of the probabilities $p_0, p_1$ of being symptomatic, and $s_0, s_1$ of exhibiting specific symptoms, within a population of interest; and (iii) the distribution of the impact $\beta$ of the risk factors for contracting the disease. In the context of the COVID-19 pandemic, such aggregated figures are pivotal to understand the dynamic of the disease and implement appropriate crisis policy.

Since the true diagnoses are hidden variables, the Expectation-Maximization algorithm \cite{Dempster1977MaximumAlgorithm} is an appealing method to perform inference in our Bayesian model; hence to achieve the aforementioned objectives. However, the EM algorithm requires the computation of an expectation over the posterior of the hidden variables, which may be intractable depending on the probability distributions defined in the model. To allow for flexibility in terms of model's design within our multimodal data integration framework, we offer to rely on the Stochastic EM algorithm (StEM) \cite{Celeux1988AAlgorithm}. StEM effectively estimates the conditional expectation in the EM using the ``Stochastic Imputation Principle", \textit{i.e.} approximating the expectation by sampling once from the underlying distribution. This method allows us to carry out inference with priors that are not necessarily the Beta distributions implemented in our experiments --- thus providing us with additional flexibility in modeling real data. Furthermore, StEM is more robust, being less dependent on the parameters' initialization than the EM -- a definite advantage given our very uncertain framework. Finally, StEM shows better asymptotic behavior: unlike the EM, StEM always leads to maximization of a complete data log-likelihood in the M-step \cite{Nielsen2000TheResults}. 








\vspace{-0.3cm}


\subsection{Stochastic E-step at iteration $(j+1)$}
\vspace{-0.1cm}

Starting iteration $j+1$, the current estimates of the model's parameters are $\theta^{(j)}$. The stochastic E-step computes the posterior distribution $\mathbb{P}(D_i | T_i, S_i, X_i, Y_i, \theta^{(j)})$ of the hidden variable $D_i$. We then sample a diagnosis $\widehat{D_i}$ from this posterior.


\begin{proposition}\label{prop:d}
The odds of the posterior of the hidden variable $D_i$ at iteration $(j+1)$ writes:
\begin{equation}
   \begin{split}
    &\frac{
        \mathbb{P}(D_i=1|T_i, S_i, X_i, Y_i, \theta^{(j)})
        }{
        \mathbb{P}(D_i=0|T_i, S_i, X_i, Y_i, \theta^{(j)})}\\
    &= \frac{
    {x^{(j)}}^{T_i} (1-x^{(j)})^{(1-T_i)}
    \times {s_1^{(j)}}^{S_iX_i} (1-s_1^{(j)})^{S_i(1-X_i)} 
    \times \pi(Y_i, \beta^{(j)})
    \times p_1^{S_i}(1-p_1)^{(1-S_i)}
    }{
    {y^{(j)}}^{T_i}(1-y^{(j)})^{(1-T_i)}\times {s_0^{(j)}}^{S_iX_i}(1-s_0^{(j)})^{S_i(1-X_i)} \times (1-\pi(Y_i, \beta^{(j)}))
        \times p_0^{S_i}(1-p_0)^{(1-S_i)}}.
    \end{split}
\end{equation}
\end{proposition}



\subsection{M-step at iteration $(j+1)$}
\vspace{-0.1cm}

The following proposition shows the updates in the model's parameters at iteration $(j+1)$, performed in the M-step.


\begin{proposition}\label{prop:xy}
The parameters updates write:
\begin{equation}
    \begin{gathered}
       x^{(j+1)}  
       = \frac{
       \sum_{i=1}^nT_i\widehat{D_i}+\alpha_x - 1}{
       \sum_{i=1}^n\widehat{D_i}+(\alpha_x+\beta_x) - 2 }, \quad
       y^{(j+1)}
       = \frac{
       \sum_{i=1}^nT_i(1-\widehat{D_i})+\alpha_y - 1}{
       \sum_{i=1}^n (1-\widehat{D_i})+(\alpha_y+\beta_y) - 2}\\
       p_0^{(j+1)} 
       = \frac{
       \sum_{i=1}^n (1-\widehat{D_i}) S_i + \alpha_{p_0} - 1}{
       \sum_{i=1}^n (1 - \widehat{D_i})+(\alpha_{p_0}+\beta_{p_0}) - 2 }, \quad
       p_1^{(j+1)} 
        = \frac{
       \sum_{i=1}^n \widehat{D_i}S_i + \alpha_{p_1} - 1}{
       \sum_{i=1}^n \widehat{D_i}+(\alpha_{p_1}+\beta_{p_1}) - 2}\\
       s_0^{(j+1)} 
       = \frac{
       \sum_{i=1, \text{s.t. } S_i=1}^n (1-\widehat{D_i}) X_i + \alpha_{s_0} - 1}{
       \sum_{i=1, \text{s.t. } S_i=1}^n (1 - \widehat{D_i})+(\alpha_{s_0}+\beta_{s_0}) - 2 }, \quad
       s_1^{(j+1)} 
        = \frac{
       \sum_{i=1, \text{s.t. } S_i=1}^n \widehat{D_i}X_i + \alpha_{s_1} - 1}{
       \sum_{i=1, \text{s.t. } S_i=1}^n \widehat{D_i}+(\alpha_{s_1}+\beta_{s_1}) - 2}\\
       \beta^{(j+1)}
        = \text{argmin}_\beta \sum_{i=1}^n \frac{||\widehat{D_i} - g(Y_i\beta)||^2}{2\sigma^2} + \frac{||\beta||^2}{2\sigma_\beta^2},
    \end{gathered}
\end{equation}
where the minimization on $\beta$ is performed through Newton-Raphson descent.
\end{proposition}


\xhdr{Algorithm Complexity} Our algorithm only relies on simple sequential updates of the distribution parameters, as highlighted in Prop.~\ref{prop:d} and \ref{prop:xy}. Denoting by $L$ the number of such parameters, and $N$ the number of samples:
\vspace{-0.4cm}

\begin{itemize}[leftmargin=2em]
\setlength\itemsep{0em}
    \item The updates for $X$ and $T$ only involve matrix multiplications of the form $\widehat{D}^TX$ , while those for $\beta$ involve matrix multiplications of size $O(M^2N)$ --- yielding a complexity of $O(LN)$.
    \item Prior updates rely solely on element-wise operations on the log-odds, with $O(N)$ complexity.
\end{itemize} 
Denoting as $B$ the maximum number of steps, the complexity is $O(BLN)$, and thus linear in $N$. Fig~\ref{fig:bench_perf}(A) provides an example of the evolution of the number of steps as the number of samples increases.

\vspace{-0.2cm}

\section{Validation on Synthetic Data}\label{sec:simu}
\vspace{-0.3cm}

Since our approach is unsupervised, we begin by validating it on synthetic datasets where the ground truth is known and controlled --- thus allowing us to characterize the performance of the algorithm, and to showcase the strength of combining multiple noisy modalities. 

 We assume the same generative model as in Fig. \ref{fig:model} and generate synthetic data for various pairs of values of sensitivity and specificity ranging from 60\% to 99\% . In each case, we simulate $N=100$ tuples  of variables ($(Y_i, X_i, T_i)_{i=1}^{n}$, for $n \in \{100, 200, 500, 1000\}$ participants. We also vary the noise level $\sigma \in \{0.1, 0.5, 1.0\}$ for the prior on $D$ in  Fig. \ref{fig:model}. To mimic our COVID data, we simulate $K=14$ symptoms and $M=2$ risk factors, for which we randomly select the parameters\footnote{All of these parameters are provided with the code in the supplementary materials.}. 


\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_SIM/sens_spe_heatmaps.pdf}
    \caption{Performance of the SEM algorithm for n=300 samples, $\sigma=0.5$, and varying levels of specificity and sensitivity. \textbf{(A)} Raw accuracy of the labels imputed via StEM. \textbf{(B)} Difference in accuracy between the StEM imputed diagnostics and the observed test variable $T$.  }
    \label{fig:heatmap_acc}
\end{figure}

\xhdr{Improving upon the sole test $T$} Fig.~\ref{fig:heatmap_acc} (A) and (B) show the Stochastic EM's raw accuracy and accuracy improvement over the sole test result $T$. Note that this difference is always positive, highlighting that our method only improves upon single diagnostic inputs --- even when one input is more reliable than any of the others.
As expected, the most substantial improvements upon $T$ are observed when the test specificity and/or sensitivity are low. For instance, for sensitivity and specificity of 70\%, our method provides an  $86.6 \pm (2.5)$\% accuracy for the diagnosis -- or equivalently, a 16\% gain over $T$.
We further quantify the algorithm's performance in Table 1 of the supplementary materials, for values of the specificity close to those observed on the field.




\xhdr{Benchmarking against other models} We now benchmark our algorithm against other approaches:
\vspace{-0.3cm}

\begin{itemize}[leftmargin=2em]
\setlength\itemsep{0em}
\setlength\topsep{0em}
    \item \textit{Vanilla Classifier:} using  both the context $Y$ and symptoms $X$ as inputs, we fit a logistic regression to the test labels $T$. We choose the regularization parameter using 10-fold cross-validation, and compute confidence intervals for the log-probability of $(D|X,T)$ by bootstrapping.
    \item \textit{Data-Agnostic EM}: we implement a deterministic version of EM, providing uninformative priors for the parameter (thereby reflecting our absence of knowledge of the truth), which are not updated --- i.e, an "uninformed'' equivalent of the belief networks found in the literature.
        \item \textit{Data-Informed EM}: similar to the Data-Agnostic EM, but choosing the priors (which then remain fixed) based on the empirical data.
\end{itemize}
\vspace{-0.3cm}

The results are shown in Fig.~\ref{fig:bench_perf}(B), and further completed in the supplementary materials. We highlight the superiority of our deeper and more adaptive hierarchical model, yielding improvements  (for a reasonable tuple of specificity and sensitivity of 80\%) of up to 4\% over the Data-Informed EM, 8\% over the Data-Agnostic one, and 9\% over the Vanilla classifier.


\xhdr{Assessing the robustness of the Stochastic EM} For a fixed specificity 80\%, Figure \ref{fig:varying_inputs} shows the accuracy of StEM for different values of $\sigma$ (A) or as the number of samples increases (B). These axes of variation seem to yield little impact on the model's performance, providing evidence that our algorithm is fairly robust, especially with respect to low-sample size.

\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_SIM/simulations.pdf}
    \caption{Performance of the Stochastic EM for a fixed specificity of 80\%,  (A) as the number of samples increases (and $\sigma=0.5$) and (B)  as the variance increases  (and $n=300$).}
    \label{fig:varying_inputs}
\end{figure}

\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_SIM/convergence_bis.pdf}
    \caption{(A) Number of EM steps until convergence; and (B) Accuracy comparison against benchmarks: Stochastic EM (StEM), EM with data-informed priors on models' parameters, not updated during learning (EM-2), EM with uninformative priors on models' parameters, not updated during learning (EM-3), Vanilla logistic regression (Log. Reg.).}
    \label{fig:bench_perf}
\end{figure}

\xhdr{Convergence} Finally, we assess the convergence properties of our algorithm. We examine the distribution of the relative difference between recovered coefficients and ground truth (expressed as a percentage of the ground-truth value). The plots, provided in the supplementary materials, show deviations that are within a few percentages of the true value of the coefficients -- thus highlighting the ability of the model to converge to the ground-truth parameters, and making it relevant from a medical perspective to characterize the disease's etiology. Illustrations of the behavior of the number of required EM steps are also provided in the supplementary materials.




\vspace{-0.3cm}

\section{Results on Real Data}\label{sec:real}
\vspace{-0.3cm}

We turn to the real dataset of $n=117$ participants described in Section~\ref{sec:data}. The purpose of this section is to show how our algorithm can be practically applied to process heterogeneous data types and inform participants and researchers alike. At the individual level, we provide each user with (a) the confirmation of the diagnostic, and (b) confidence intervals associated with the uncertainty to shed more light on the uncertainty associated with the diagnostic. At the global level, we provide policymakers with (c) aggregated analysis of the herd immunity, and associated measures of uncertainty. 

\xhdr{At the subject-level} We present two examples, where our algorithm either confirms or infirms the result of the test -- thereby allowing for the potential flagging of false negatives. The first example (Fig.~\ref{fig:real} A) is a user that registers a negative test while being asymptomatic and with a limited number of risk factors. In this case, we expect our model to confirm the result of the test, and provide a narrower confidence interval as per the probability of immunity -- as confirmed by Fig.~\ref{fig:real} (A). The second example showcases an instance where the questionnaire and the test disagree. Subject 108 is a user that registers a negative test, while exhibiting a wide number of known COVID symptoms (dry cough, shortness of breath, fever), but took the test less than 10 days after his illness. While the posterior does not reclassify the subjects' diagnosis, the confidence interval associated with the prediction of immunity reflects the uncertainty associated with this case, and flags it as a potential false negative.





\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{FIG_REAL/examples_short_nm.pdf}
    \caption{Posterior diagnosis distributions on two selected subjects. In each panel (A-B): the left plot represents the posterior of the diagnosis, given the symptoms and risk factors data reported in the questionnaire, while the right plot is the posterior of the diagnosis, given the symptoms, risk factors, and reported result of the diagnostic test. The red dot represents the expectation.}
    \label{fig:real}
\end{figure}


\xhdr{At the global level} At the global level, this framework allows to perform principled inference for the disease and the population of interest.
For this  cohort of $n=117$ healthcare workers in the UK, at-home tests predict $35.8\%$ of immunity. 
The posterior distributions of the model's parameters shed light on the actual accuracy of the LFA tests on our population. Fig.~\ref{fig:global} (A) compares the posteriors of the sensitivity and specificity to the priors built from values reported by manufacturers, for asymptomatic subjects. Furthermore, our results provide information regarding the most prevalent symptoms for COVID-19. For instance, among symptomatic participants, the probability of exhibiting fever is higher for infected subjects (see Fig.~\ref{fig:global} (B)) while cough with sputum does not seem to be associated with COVID-19, being probably the sign of another infection (see Fig.~\ref{fig:global} (C)). Similar plots on additional symptoms are provided in the supplementary materials.

\vspace{-0.2cm}

\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_REAL/global_real.pdf}
    \caption{(A) Posteriors of sensitivity and specificity compared to their priors, for asymptomatic subjects. (B-C): Posteriors of the probability of exhibiting specific symptoms, for symptomatic subjects with estimated negative ($D=0$) or positive ($D=1$) diagnosis.}
    \label{fig:global}
\end{figure}

\vspace{-0.3cm}
\section{Conclusion}
\vspace{-0.3cm}

This applications paper provides a statistical framework for multimodal integration of noisy diagnostic test results, self-reported symptoms, and demographic variables. Compared to previous approaches, our work is more amenable to the handling of the different inputs' uncertainty. While we have focused here on binary symptoms, this adaptive Bayesian framework, together with the flexibility provided by the Stochastic EM algorithm, pave the way for the integration of other continuous-valued inputs (index of the severity of the disease, pain, etc.) as well as interactions -- an extension that we are currently in the process of implementing. The robustness of the Stochastic EM is crucial in ensuring the reliability of the parameters given the medical stakes. Finally, we note that the generative model allows for principled handling of missing data, thus allowing us to make diagnoses even in the presence of incomplete data.






\section{Introduction and Related Work}
\vspace{-0.3cm}

Current medical diagnoses are most often based on the combination of several data inputs by medical experts, typically including (i) clinical history, interviews, and physical exams, (ii) laboratory tests, (iii) electrophysiological signals, and medical images. Advances in the machine learning (ML) community have highlighted the potential of ML to contribute to the field of Computer-Aided Diagnosis (CAD), for which we distinguish two main classes of methods.


\xhdr{Single modality analysis} Most recent efforts in the ML community have focused on analyzing a single medical data source --- called a ``modality". For instance, the parsing of electronic health records (EHR) for clinical history, patients' interviews, and physical exams has shown huge potential for the diagnosis of a broad range of diseases, from coronary artery disease to rheumatoid arthritis \cite{Harrell1984RegressionPrediction, Kurt2008ComparingDisease, Carroll2011NaiveArthritis., Hippisley-Cox2013PredictingScore, Rahimian2018PredictingRecords}. ML algorithms have been developed to process various types of laboratory results, from urine steroid profiles to gene expression data \cite{Huang2018NovelClassification, Wilkes2018UsingProfiles, Demirci2016ArtificialAlgorithms, Agrahari2018ApplicationsMalignancies}. Others have shown success in automatically processing and classifying medical signals such as electrocardiograms (ECG) \cite{Acharya2017AutomatedNetwork, Kiranyaz2016Real-TimeNetworks, Rad2017ECG-BasedAnalysis, Pawiak2018NovelSystem}, or electroencephalograms (EEG) \cite{Richhariya2018EEGMachine}. Finally, a large body of work has focused on medical imaging analysis encompassing a vast number of tasks, such as automatic extraction of diagnostic features \cite{Thomaz2017FeatureImages, Bar2018ChestTraining}, segmentation of anatomical structures \cite{Moeskops2016AutomaticNetwork, Li2015AutomaticNetworks, Shelhamer2017FullySegmentation, Milletari2016V-Net:Segmentation, Dolz20183DStudy, Christ1919AutomaticFields}, or direct diagnosis through image classification \cite{Anthimopoulos2016LungNetwork, Miki2017ClassificationNetwork}. 

While these works achieve record-breaking diagnostic accuracy, they often rely on supervised learning approaches -- requiring the diagnosis ground-truth to be available during training -- and on the acquisition of large datasets. Furthermore, they focus solely on a single type of data input (EEGs, scans, etc.) -- often acquired by clinicians using specialized, highly accurate equipment -- and do not harvest the potentially rich and complementary sources of information provided by alternative medical modalities. Yet, with the development of at-home diagnostic tests (lateral flow assays, questionnaire data for disease screening, mobile health apps, etc.), this paradigm shifts, and diagnoses have to be established through the combination of multiple sources of cheaper, yet often noisier and more imprecise data. 
In light of the uncertainty exhibited by these various inputs, it also becomes indispensable to pair the provision of a diagnosis with a notion of a confidence interval, thereby thoroughly characterizing our state of knowledge given the available data. 


\xhdr{Multiple modality integration} Interestingly, diagnosis studies fusing different data sources seemed to be more prominent in the first years of CAD. Bayesian networks \cite{Neal1993ProbabilisticMethods} -- also called belief networks -- have been used as a crucial decision tool for automatic diagnosis. Such networks provide a biologically-grounded and interpretable statistical framework that integrates the probabilities of contracting the disease given a patient's clinical history, and the probability of developing symptoms or observing positive test results in the presence or absence of the disease. The advantages of these Bayesian networks are that they are (a) fully unsupervised and do not require ground-truth information, and (b) able to provide meaningful results even in low sample regimes.

Bayesian networks have thus been implemented in a variety of contexts to integrate clinical data and laboratory results, and diagnose conditions ranging from pyloric stenosis to acute appendicitis \cite{Alvarez2006EvaluationStenosis, Gevaert2006PredictingNetworks, Sakai2007AccuracyModel, Gonzalez-Lopez2016DevelopmentUveitis,Sangamuang2018AccuracyDecision}. They have also been used to combine clinical data with medical images, and subsequently applied to assess venous thromboembolism \cite{Kline1984DerivationThromboembolism}, community-acquired pneumonia \cite{Aronsky1998DiagnosingNetwork.}, head injury prognostics \cite{Sakellaropoulos1999DevelopmentTechniques}, or to predict tumors \cite{Wang1999Computer-assistedNetwork,Kahn1997ConstructionCancer, Sneha2013TowardsImages}. As early as 1994, full integration of clinical data, laboratory results and imaging features was performed to diagnose gallbladder disease \cite{Haddawy1994ADisease}. 

However, contrary to our proposed setting, the uncertainty that these Bayesian networks integrate is typically known and controlled. The medical conditions that they study are well characterized by a set of specific questions and physiological exams (e.g. projectile vomiting, potassium levels and ultrasounds in the case of Pyloric stenosis). Not only do these inputs provide a very strong signal for the diagnosis, but the uncertainty arising from the different modalities can often be reliably informed by a test manufacturer, extensive medical research or prior ML studies (for diagnostic inputs obtained through classification algorithms). The uncertainty of a diagnostic method, however, is not always well specified: whether it be (i) a physician's assessment of the symptoms associated with a novel or rare disease, (ii) predictions of an ML algorithm whose accuracy has not been fully characterized outside of a curated research environment or reference datasets (such as MNIST, CIFAR), or simply (iii) biological tests whose sensitivity still has to be determined. 




\xhdr{Multiple noisy modality integration: a COVID-19 case study} To motivate this paper and further understand the limitations of these previous approaches, let us consider a particular use case: COVID-19 antibody testing. 
Lateral flow immunoassay (LFA) antibody tests for COVID-19 are one of the manageable, affordable and easily deployable options to allow at-home testing for large population and provide assessments of our immunity.
Yet, studies have shown that the sensitivity of these tests remains highly variable and highly contingent on the time of testing.
The successful deployment of LFAs thus depends on their augmentation with additional data inputs, such as user-specific risk factors and self-reported symptoms.
The provision of confidence scores is essential to flag potential false negative or positive tests (requiring re-testing or closer scrutiny) and to assess local prevalence levels -- both pivotal for researchers and policymakers in the context of a pandemic.

\xhdr{Contributions and outline} This applications paper is geared towards the practical integration of noisy sources of information for CAD. Our contribution is two-fold. From a methods perspective, we account for the variability of the inputs by devising a two-level Bayesian hierarchical model.
In contrast to existing Bayesian methods for CAD, our model is deeper, and trained using a Stochastic Expectation Maximization (EM) algorithm \cite{celeux1995stochastic,Celeux2001OnAlgorithm,nielsen2000stochastic}. The Stochastic EM allows to overcome the limitations of its non-stochastic counterpart \cite{celeux1995stochastic}, that is (a) its sensitivity to the starting position, (b) its potential slow convergence rate, and (c) its possible convergence to a saddle position instead of a local maximum. From an applications perspective, we gear this algorithm to enhance at-home LFA testing in the context of the COVID-19 pandemic. In particular, we wish to (a) quantify the benefit of multimodal data integration when the diagnostics are uncertain, and (b) show how our method can benefit medical experts or researchers in real life. 
Section~\ref{sec:multimodal} presents the Bayesian model for multimodal integration and the Stochastic Expectation-Maximization algorithm that performs principled and scalable inference. Section~\ref{sec:simu} presents extensive tests of our model on simulated datasets. Section~\ref{sec:real} details the results obtained on the Covid-19 dataset and the impact of our method for affordable and reliable at-home test kits.




\section{Covid-19 Dataset}\label{sec:data}

By way of clarifying the challenges that our algorithm proposes to overcome, we present here the COVID Clearance Remote Recovery \& Antibody Monitoring Platform study 
\footnote{COV-CLEAR, \url{www.cov-clear.com}}, which motivated our approach. The purpose of this study is to track the evolution of the immunity of a cohort of adult participants in the UK with various COVID-19 exposure risks. This paper focuses on a subset of 117 healthcare workers, for which we have both questionnaire information and LFA test results.


\xhdr{Home-testing Immunoassay Data} Participants are issued packs of home-testing kits with written instructions. These kits identify Immunoglobulin M (IgM) and/or Immunoglobulin G (IgG) specific to SARS-CoV2 in blood samples using three gold-standard methods \footnote{Home-sampled (blood) Abbott ARCHITECT® Anti-SARS-CoV-2 IgG CMIA CE Marked; PHE Approved} \footnote{Elecsys® Anti-SARS-CoV-2 Double Antigen, CE Marked, PHE Approved} \footnote{Biopanda Ltd  Product Number: RAPG-COV-019 version 1 In-Vitro LFT IgM IgG Antibody test. CE certified for Professional Use Only.}. Pictures and typical examples of LFA test results are provided in the supplementary materials. For the LFAs, participants self-report a positive result if IgM and/or IgG are detected by the test in addition to a positive control band. Through the collection of this data, the COV-CLEAR study aims to address questions relating to  (a) the quantification of the robustness of the antibody response, and (b) the durability of this response. However, while affordable, we highlight that the LFA tests suffer from vastly varying levels of sensitivity --- registering a sensitivity as low as 70\% on asymptomatic individuals.

\xhdr{Clinical Data: Questionnaire} Additionally, participants are asked to answer a questionnaire, ideally before knowing the result of their LFA test. The form consists of questions related to $K=14$ exhibited symptoms (fever, cough, runny nose, etc.) and $M=2$ subject-specific risk factors (household size and proportion of members with a suspected or confirmed history of COVID-19). The empirical distributions of the symptoms and the risk factors in this cohort are provided in the supplementary materials. 







\section{Multimodal Data Integration: Statistical Model}\label{sec:multimodal}
\vspace{-0.3cm}

In this section, we describe our principled integration of noisy diagnostic test results, with additional clinical data such as symptom data and subject-specific risk factors.
Our approach applies to general heterogeneous medical data where the outputs are binary. The latter could be either self-reported answers to questionnaires, clinician-reported physiological exams, outputs of a diagnostic based on an image, abnormalities of laboratory results, etc., making this a widespread and general setting. Thus, while we implement and showcase the method for the particular purpose of applying it to antibody testing, this could be relevant to any medical diagnostic with binary inputs.

Denote $D$ the true diagnosis of an individual (healthy/sick), $T$ the outcome of the noisy diagnostic test (positive/negative), $S$ the symptomaticity (symptomatic/asymptomatic), $X$ the symptoms exhibited and $Y$ the subject-specific risks factors. The underlying assumption is that given a true diagnosis $D$, the symptoms $X$ and the diagnostic test outcome $T$ are independent. In other words, the probability of the diagnostic test being a false negative $\mathbb{P}[T=0|D=1]$ is independent of the symptoms of a truly infected individual $\mathbb{P}[X|D=1]$. Similarly, given a diagnosis $D$, the test outcome $T$ and the exhibited symptoms $X$ are independent of the risk factors $Y$. We define:
\begin{itemize}[leftmargin=3em]
       \item $\mathbb{P}(D=1|Y) = \pi_\beta(Y)$ the probability of contracting the disease given risk factors $Y$,
       \item $\mathbb{P}[(T=1|D=1) = x$ the sensitivity of the diagnostic test,
       \item $\mathbb{P}[T=0|D=0] = 1-y$ the specificity of the diagnostic test, \textit{i.e.} $y=\mathbb{P}[T=1|D=0]$,
       \item $\mathbb{P}(S=1|D=0) = p_0$ the probability of being symptomatic when whilst not having been infected by that specific disease (the symptoms could be due to another illness for instance),
       \item $\mathbb{P}(S=1|D=1) = p_1$ the probability of having been symptomatic upon infection,
       \item $\mathbb{P}(X_k=1|S=1,D=0) = s_{0k}$ the probability of exhibiting symptom $k$ when not infected,
       \item $\mathbb{P}(X_k=1|S=1,D=1) = s_{1k}$ the probability of exhibiting symptom $k$ upon infection.
\end{itemize}
In the above, the diagnostic test may be any biological test (e.g. LFA antibody tests), or output of a medical imaging classification algorithm provided by a ML research team. We have here splitted our binary variables into two classes according to their variability: (a) the test $T$, for which we have provisional estimates for the sensitivity and specificity (given by the manufacturer) but which we are still uncertain, and (b) the symptoms $X$, which carry additional uncertainty in that neither of them is extremely specific to COVID-19 nor is their prevalence amongst COVID cohorts very well established; Values for $(x, y)$, $p_0, p_1, s_0, s_1$ or $\beta$ may be published, but challenged by complementary research or field experience, or may be completely unavailable.
We thus need to account for these inputs' variability and for the relative importance of different risk factors $\beta$. 
We propose using a hierarchical Bayesian network (see Fig.~\ref{fig:model}) that is consequently deeper than the ones traditionally implemented in the CAD literature \cite{Alvarez2006EvaluationStenosis, Gevaert2006PredictingNetworks, Sakai2007AccuracyModel, Gonzalez-Lopez2016DevelopmentUveitis,Sangamuang2018AccuracyDecision}.
 
The uncertainty in the test sensitivity $x$ and specificity $1-y$ are expressed through Beta priors on $x$ and $y$ with parameters $(\alpha_{x}, \beta_x)$ and $(\alpha_y, \beta_y)$, see Fig.~\ref{fig:model}. Similarly, the uncertainty on the probabilities of symptomaticity and of the different symptoms are expressed via Beta priors with respective parameters $(\alpha_{p_0}, \beta_{p_0}), (\alpha_{p_1}, \beta_{p_1}), \{(\alpha_{s_{0k}}, \beta_{s_{0k}}), (\alpha_{s_{1k}}, \beta_{s_{1k}})\}_k$, see Fig.~\ref{fig:model}. When published estimates exist for $x, y, p_0, p_1, s_0$ and $s_1$ (e.g. as provided by the LFA test manufacturer), we match the moments of the Beta priors with those of the published distributions.
If this information is unavailable, we use the non-informative Beta prior with parameters $(0.5, 0.5)$. In the COVID-19 example, we implement a non-informative prior for the probabilities of symptomaticity and symptoms, as the etiology of the disease remains unknown. These priors are then updated during learning as we aggregate the information from the observed and imputed variables.

Lastly, we model the probability $\pi_\beta(Y)$ of contracting the disease with a logistic regression on the $M$ risk factors $Y$. In the case of the COVID-19 dataset, these include county data, profession, size of household, etc. The coefficients $\beta=\{\beta_m\}_m$ weight the relative importance of each factor $Y_m$, $m=1...M$. We express our uncertainty on the possible impact of a given factor $m$ by introducing a Gaussian prior: $\beta \sim N(0, \sigma_\beta^2)$, see Fig.~\ref{fig:model}.



\begin{figure}[h!]
\centering
\def0.8\columnwidth{0.8\columnwidth}
\input{./total_model_with_legend.pdf_tex}
    \caption{Hierarchial bayesian network integrating the accuracy uncertainty on the data sources, to estimate the true diagnosis $D$. The index $k=1...K$ represents symptom $k$, while $m=1...M$ represents risk factor $m$. The shaded cells represent observed variables or known hyper-parameters. The dashed lines represent switches.}
    \label{fig:model}
\end{figure}
\vspace{-0.3cm}


\section{Inference in the Hierarchical Bayesian Network via Stochastic EM}
\vspace{-0.3cm}

\textit{At the subject level}, we seek to compute the posterior distribution of the true diagnosis $D_i$, informed by the integration of the observed variables $T_i, S_i, X_i, Y_i$. This posterior yields an estimated diagnosis, together with a credible interval that expresses the confidence associated with our prediction. In the case of COVID-19, this provides individual estimates of each citizen's immunity that may inform their social behavior as lockdown is eased.

\textit{At the global level}, we wish to learn: (i) the distributions of the sensitivity $x$ and specificity $y$ of the diagnostic test as observed in the field; (ii) the distributions of the probabilities $p_0, p_1$ of being symptomatic, and $s_0, s_1$ of exhibiting specific symptoms, within a population of interest; and (iii) the distribution of the impact $\beta$ of the risk factors for contracting the disease. In the context of the COVID-19 pandemic, such aggregated figures are pivotal to understand the dynamic of the disease and implement appropriate crisis policy.

Since the true diagnoses are hidden variables, the Expectation-Maximization algorithm \cite{Dempster1977MaximumAlgorithm} is an appealing method to perform inference in our Bayesian model; hence to achieve the aforementioned objectives. However, the EM algorithm requires the computation of an expectation over the posterior of the hidden variables, which may be intractable depending on the probability distributions defined in the model. To allow for flexibility in terms of model's design within our multimodal data integration framework, we offer to rely on the Stochastic EM algorithm (StEM) \cite{Celeux1988AAlgorithm}. StEM effectively estimates the conditional expectation in the EM using the ``Stochastic Imputation Principle", \textit{i.e.} approximating the expectation by sampling once from the underlying distribution. This method allows us to carry out inference with priors that are not necessarily the Beta distributions implemented in our experiments --- thus providing us with additional flexibility in modeling real data. Furthermore, StEM is more robust, being less dependent on the parameters' initialization than the EM -- a definite advantage given our very uncertain framework. Finally, StEM shows better asymptotic behavior: unlike the EM, StEM always leads to maximization of a complete data log-likelihood in the M-step \cite{Nielsen2000TheResults}. 








\vspace{-0.3cm}


\subsection{Stochastic E-step at iteration $(j+1)$}
\vspace{-0.1cm}

Starting iteration $j+1$, the current estimates of the model's parameters are $\theta^{(j)}$. The stochastic E-step computes the posterior distribution $\mathbb{P}(D_i | T_i, S_i, X_i, Y_i, \theta^{(j)})$ of the hidden variable $D_i$. We then sample a diagnosis $\widehat{D_i}$ from this posterior.


\begin{proposition}\label{prop:d}
The odds of the posterior of the hidden variable $D_i$ at iteration $(j+1)$ writes:
\begin{equation}
   \begin{split}
    &\frac{
        \mathbb{P}(D_i=1|T_i, S_i, X_i, Y_i, \theta^{(j)})
        }{
        \mathbb{P}(D_i=0|T_i, S_i, X_i, Y_i, \theta^{(j)})}\\
    &= \frac{
    {x^{(j)}}^{T_i} (1-x^{(j)})^{(1-T_i)}
    \times {s_1^{(j)}}^{S_iX_i} (1-s_1^{(j)})^{S_i(1-X_i)} 
    \times \pi(Y_i, \beta^{(j)})
    \times p_1^{S_i}(1-p_1)^{(1-S_i)}
    }{
    {y^{(j)}}^{T_i}(1-y^{(j)})^{(1-T_i)}\times {s_0^{(j)}}^{S_iX_i}(1-s_0^{(j)})^{S_i(1-X_i)} \times (1-\pi(Y_i, \beta^{(j)}))
        \times p_0^{S_i}(1-p_0)^{(1-S_i)}}.
    \end{split}
\end{equation}
\end{proposition}



\subsection{M-step at iteration $(j+1)$}
\vspace{-0.1cm}

The following proposition shows the updates in the model's parameters at iteration $(j+1)$, performed in the M-step.


\begin{proposition}\label{prop:xy}
The parameters updates write:
\begin{equation}
    \begin{gathered}
       x^{(j+1)}  
       = \frac{
       \sum_{i=1}^nT_i\widehat{D_i}+\alpha_x - 1}{
       \sum_{i=1}^n\widehat{D_i}+(\alpha_x+\beta_x) - 2 }, \quad
       y^{(j+1)}
       = \frac{
       \sum_{i=1}^nT_i(1-\widehat{D_i})+\alpha_y - 1}{
       \sum_{i=1}^n (1-\widehat{D_i})+(\alpha_y+\beta_y) - 2}\\
       p_0^{(j+1)} 
       = \frac{
       \sum_{i=1}^n (1-\widehat{D_i}) S_i + \alpha_{p_0} - 1}{
       \sum_{i=1}^n (1 - \widehat{D_i})+(\alpha_{p_0}+\beta_{p_0}) - 2 }, \quad
       p_1^{(j+1)} 
        = \frac{
       \sum_{i=1}^n \widehat{D_i}S_i + \alpha_{p_1} - 1}{
       \sum_{i=1}^n \widehat{D_i}+(\alpha_{p_1}+\beta_{p_1}) - 2}\\
       s_0^{(j+1)} 
       = \frac{
       \sum_{i=1, \text{s.t. } S_i=1}^n (1-\widehat{D_i}) X_i + \alpha_{s_0} - 1}{
       \sum_{i=1, \text{s.t. } S_i=1}^n (1 - \widehat{D_i})+(\alpha_{s_0}+\beta_{s_0}) - 2 }, \quad
       s_1^{(j+1)} 
        = \frac{
       \sum_{i=1, \text{s.t. } S_i=1}^n \widehat{D_i}X_i + \alpha_{s_1} - 1}{
       \sum_{i=1, \text{s.t. } S_i=1}^n \widehat{D_i}+(\alpha_{s_1}+\beta_{s_1}) - 2}\\
       \beta^{(j+1)}
        = \text{argmin}_\beta \sum_{i=1}^n \frac{||\widehat{D_i} - g(Y_i\beta)||^2}{2\sigma^2} + \frac{||\beta||^2}{2\sigma_\beta^2},
    \end{gathered}
\end{equation}
where the minimization on $\beta$ is performed through Newton-Raphson descent.
\end{proposition}


\xhdr{Algorithm Complexity} Our algorithm only relies on simple sequential updates of the distribution parameters, as highlighted in Prop.~\ref{prop:d} and \ref{prop:xy}. Denoting by $L$ the number of such parameters, and $N$ the number of samples:
\vspace{-0.4cm}

\begin{itemize}[leftmargin=2em]
\setlength\itemsep{0em}
    \item The updates for $X$ and $T$ only involve matrix multiplications of the form $\widehat{D}^TX$ , while those for $\beta$ involve matrix multiplications of size $O(M^2N)$ --- yielding a complexity of $O(LN)$.
    \item Prior updates rely solely on element-wise operations on the log-odds, with $O(N)$ complexity.
\end{itemize} 
Denoting as $B$ the maximum number of steps, the complexity is $O(BLN)$, and thus linear in $N$. Fig~\ref{fig:bench_perf}(A) provides an example of the evolution of the number of steps as the number of samples increases.

\vspace{-0.2cm}

\section{Validation on Synthetic Data}\label{sec:simu}
\vspace{-0.3cm}

Since our approach is unsupervised, we begin by validating it on synthetic datasets where the ground truth is known and controlled --- thus allowing us to characterize the performance of the algorithm, and to showcase the strength of combining multiple noisy modalities. 

 We assume the same generative model as in Fig. \ref{fig:model} and generate synthetic data for various pairs of values of sensitivity and specificity ranging from 60\% to 99\% . In each case, we simulate $N=100$ tuples  of variables ($(Y_i, X_i, T_i)_{i=1}^{n}$, for $n \in \{100, 200, 500, 1000\}$ participants. We also vary the noise level $\sigma \in \{0.1, 0.5, 1.0\}$ for the prior on $D$ in  Fig. \ref{fig:model}. To mimic our COVID data, we simulate $K=14$ symptoms and $M=2$ risk factors, for which we randomly select the parameters\footnote{All of these parameters are provided with the code in the supplementary materials.}. 


\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_SIM/sens_spe_heatmaps.pdf}
    \caption{Performance of the SEM algorithm for n=300 samples, $\sigma=0.5$, and varying levels of specificity and sensitivity. \textbf{(A)} Raw accuracy of the labels imputed via StEM. \textbf{(B)} Difference in accuracy between the StEM imputed diagnostics and the observed test variable $T$.  }
    \label{fig:heatmap_acc}
\end{figure}

\xhdr{Improving upon the sole test $T$} Fig.~\ref{fig:heatmap_acc} (A) and (B) show the Stochastic EM's raw accuracy and accuracy improvement over the sole test result $T$. Note that this difference is always positive, highlighting that our method only improves upon single diagnostic inputs --- even when one input is more reliable than any of the others.
As expected, the most substantial improvements upon $T$ are observed when the test specificity and/or sensitivity are low. For instance, for sensitivity and specificity of 70\%, our method provides an  $86.6 \pm (2.5)$\% accuracy for the diagnosis -- or equivalently, a 16\% gain over $T$.
We further quantify the algorithm's performance in Table 1 of the supplementary materials, for values of the specificity close to those observed on the field.




\xhdr{Benchmarking against other models} We now benchmark our algorithm against other approaches:
\vspace{-0.3cm}

\begin{itemize}[leftmargin=2em]
\setlength\itemsep{0em}
\setlength\topsep{0em}
    \item \textit{Vanilla Classifier:} using  both the context $Y$ and symptoms $X$ as inputs, we fit a logistic regression to the test labels $T$. We choose the regularization parameter using 10-fold cross-validation, and compute confidence intervals for the log-probability of $(D|X,T)$ by bootstrapping.
    \item \textit{Data-Agnostic EM}: we implement a deterministic version of EM, providing uninformative priors for the parameter (thereby reflecting our absence of knowledge of the truth), which are not updated --- i.e, an "uninformed'' equivalent of the belief networks found in the literature.
        \item \textit{Data-Informed EM}: similar to the Data-Agnostic EM, but choosing the priors (which then remain fixed) based on the empirical data.
\end{itemize}
\vspace{-0.3cm}

The results are shown in Fig.~\ref{fig:bench_perf}(B), and further completed in the supplementary materials. We highlight the superiority of our deeper and more adaptive hierarchical model, yielding improvements  (for a reasonable tuple of specificity and sensitivity of 80\%) of up to 4\% over the Data-Informed EM, 8\% over the Data-Agnostic one, and 9\% over the Vanilla classifier.


\xhdr{Assessing the robustness of the Stochastic EM} For a fixed specificity 80\%, Figure \ref{fig:varying_inputs} shows the accuracy of StEM for different values of $\sigma$ (A) or as the number of samples increases (B). These axes of variation seem to yield little impact on the model's performance, providing evidence that our algorithm is fairly robust, especially with respect to low-sample size.

\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_SIM/simulations.pdf}
    \caption{Performance of the Stochastic EM for a fixed specificity of 80\%,  (A) as the number of samples increases (and $\sigma=0.5$) and (B)  as the variance increases  (and $n=300$).}
    \label{fig:varying_inputs}
\end{figure}

\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_SIM/convergence_bis.pdf}
    \caption{(A) Number of EM steps until convergence; and (B) Accuracy comparison against benchmarks: Stochastic EM (StEM), EM with data-informed priors on models' parameters, not updated during learning (EM-2), EM with uninformative priors on models' parameters, not updated during learning (EM-3), Vanilla logistic regression (Log. Reg.).}
    \label{fig:bench_perf}
\end{figure}

\xhdr{Convergence} Finally, we assess the convergence properties of our algorithm. We examine the distribution of the relative difference between recovered coefficients and ground truth (expressed as a percentage of the ground-truth value). The plots, provided in the supplementary materials, show deviations that are within a few percentages of the true value of the coefficients -- thus highlighting the ability of the model to converge to the ground-truth parameters, and making it relevant from a medical perspective to characterize the disease's etiology. Illustrations of the behavior of the number of required EM steps are also provided in the supplementary materials.




\vspace{-0.3cm}

\section{Results on Real Data}\label{sec:real}
\vspace{-0.3cm}

We turn to the real dataset of $n=117$ participants described in Section~\ref{sec:data}. The purpose of this section is to show how our algorithm can be practically applied to process heterogeneous data types and inform participants and researchers alike. At the individual level, we provide each user with (a) the confirmation of the diagnostic, and (b) confidence intervals associated with the uncertainty to shed more light on the uncertainty associated with the diagnostic. At the global level, we provide policymakers with (c) aggregated analysis of the herd immunity, and associated measures of uncertainty. 

\xhdr{At the subject-level} We present two examples, where our algorithm either confirms or infirms the result of the test -- thereby allowing for the potential flagging of false negatives. The first example (Fig.~\ref{fig:real} A) is a user that registers a negative test while being asymptomatic and with a limited number of risk factors. In this case, we expect our model to confirm the result of the test, and provide a narrower confidence interval as per the probability of immunity -- as confirmed by Fig.~\ref{fig:real} (A). The second example showcases an instance where the questionnaire and the test disagree. Subject 108 is a user that registers a negative test, while exhibiting a wide number of known COVID symptoms (dry cough, shortness of breath, fever), but took the test less than 10 days after his illness. While the posterior does not reclassify the subjects' diagnosis, the confidence interval associated with the prediction of immunity reflects the uncertainty associated with this case, and flags it as a potential false negative.





\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{FIG_REAL/examples_short_nm.pdf}
    \caption{Posterior diagnosis distributions on two selected subjects. In each panel (A-B): the left plot represents the posterior of the diagnosis, given the symptoms and risk factors data reported in the questionnaire, while the right plot is the posterior of the diagnosis, given the symptoms, risk factors, and reported result of the diagnostic test. The red dot represents the expectation.}
    \label{fig:real}
\end{figure}


\xhdr{At the global level} At the global level, this framework allows to perform principled inference for the disease and the population of interest.
For this  cohort of $n=117$ healthcare workers in the UK, at-home tests predict $35.8\%$ of immunity. 
The posterior distributions of the model's parameters shed light on the actual accuracy of the LFA tests on our population. Fig.~\ref{fig:global} (A) compares the posteriors of the sensitivity and specificity to the priors built from values reported by manufacturers, for asymptomatic subjects. Furthermore, our results provide information regarding the most prevalent symptoms for COVID-19. For instance, among symptomatic participants, the probability of exhibiting fever is higher for infected subjects (see Fig.~\ref{fig:global} (B)) while cough with sputum does not seem to be associated with COVID-19, being probably the sign of another infection (see Fig.~\ref{fig:global} (C)). Similar plots on additional symptoms are provided in the supplementary materials.

\vspace{-0.2cm}

\begin{figure}
\centering
    \includegraphics[width=\textwidth]{FIG_REAL/global_real.pdf}
    \caption{(A) Posteriors of sensitivity and specificity compared to their priors, for asymptomatic subjects. (B-C): Posteriors of the probability of exhibiting specific symptoms, for symptomatic subjects with estimated negative ($D=0$) or positive ($D=1$) diagnosis.}
    \label{fig:global}
\end{figure}

\vspace{-0.3cm}
\section{Conclusion}
\vspace{-0.3cm}

This applications paper provides a statistical framework for multimodal integration of noisy diagnostic test results, self-reported symptoms, and demographic variables. Compared to previous approaches, our work is more amenable to the handling of the different inputs' uncertainty. While we have focused here on binary symptoms, this adaptive Bayesian framework, together with the flexibility provided by the Stochastic EM algorithm, pave the way for the integration of other continuous-valued inputs (index of the severity of the disease, pain, etc.) as well as interactions -- an extension that we are currently in the process of implementing. The robustness of the Stochastic EM is crucial in ensuring the reliability of the parameters given the medical stakes. Finally, we note that the generative model allows for principled handling of missing data, thus allowing us to make diagnoses even in the presence of incomplete data.





