\section{Introduction}

Satellite Earth observation (EO) data can provide critical, relevant, and timely information in support of agricultural monitoring and food security \citep{BeckerReshef2019,NationalAcademies2019}. EO can enable crop yield assessment and forecasting \citep{Becker2010,Franch2015,Franch2017,Lobell2013,Porter2017}, inform agricultural commodities markets, enable early warning and mitigation of impending crop shortages, and inform decisions about subsidies in times of crisis such as extreme drought \citep{Dorward2011} or outbreaks of disease or pests. However, the lack of up-to-date, high-resolution cropland and crop type maps is a critical barrier to agricultural and food security assessments and operational rapid response to humanitarian crises. These data are especially lacking in regions dominated by smallholder agriculture that are most vulnerable to food insecurity (e.g., sub-saharan Africa).

Machine learning techniques have commonly been used to identify cropland in EO data. The majority of studies use tree-based classifiers (primarily random forests or decision trees) or neural network/deep learning methods (primarily recurrent or convolutional neural networks). However, prior methods have typically been applied to small and/or spatially homogeneous areas \citep{RNN_Camargue, RNN_winter, Skakun2019}, require large training datasets (hundreds of thousands of labels) \citep{breizhcrops2020,Wang2020,lstm_crop_identification, RF_france}, or use EO inputs with insufficient spatial resolution for reliably detecting smallholder farms which often have field sizes $<$1 ha \citep{GFSAD,Teluguntla2018,Oliphant2019,CopernicusLandCover,Samasse2020,Waldner2017}. In many regions, particularly developing countries and regions dominated by smallholder farming, the lack of readily accessible ground truth data is a major barrier for training and validating machine learning techniques for cropland classification \citep{Coutu2020}.

In this paper, we demonstrate a novel method for rapidly generating cropland maps for areas with little to no available labeled data over a large heterogeneous area. This is achieved by leveraging global and local datasets. The global dataset contains tens of thousands of crowdsourced labels from diverse geographies \citep{GeoWiki} to learn general features for identifying cropland globally. These global examples supplement a much smaller dataset of curated, local examples (manually labeled using photointerpretation of high-resolution satellite imagery) to learn features specific to cropland in the target mapping region. Combined, these datasets provide sufficient examples to train a Long Short Term Memory network (LSTM) \citep{LSTM} to predict the presence of cropland in each pixel given a monthly time-series of Sentinel-2 multispectral observations.

The COVID-19 global pandemic has severely impacted the global food system, and the threat of food shortages and rising food prices has required countries at risk of food insecurity to rapidly design and enact aid programs in response. In the West African country of Togo, the government announced a program to boost national food production in response to the COVID-19 crisis by distributing aid to farmers, but found that high-resolution spatial information about the distribution of farms across Togo (which primarily consist of smallholder farms under 1 ha) was a critical piece of missing information for designing this program. Using our method, we were able to generate an accurate high-resolution (10 m) map of cropland in Togo for 2019 in under 10 days to respond to the Togolese government's immediate need for information about the distribution of farms across the country. We showed that our model outperforms existing publicly-available cropland maps for Togo \citep{GFSAD, CopernicusLandCover} while also being higher resolution and more recent. We have delivered this product to the Togolese government where it is being used to help inform decision-making in response to the COVID-19 pandemic, thus demonstrating a successful transition of machine learning research to operational rapid response for a real humanitarian crisis. We have made our map as well as all data and code used in this study publicly available to facilitate future research and operational systems in data-sparse regions.
\begin{figure}
  \centering
  \includegraphics[width=0.75\linewidth]{diagrams/label_distribution.pdf}
  \caption{Distribution of labeled examples in training (a) and test (b) datasets.}
  \label{fig:labels}
\end{figure}

\section{Methods}

\subsection{Labeled Data}
\subsubsection{Crowdsourced Labels}
We leveraged a large (35,866 labels) crowd-sourced dataset with labels of ``crop'' or ``non-crop'' from diverse, globally-distributed locations \citep{GeoWiki} to help the model learn features that are useful for detecting a wide range of cropland globally. Details about sample selection, labeling procedures, and quality assessment/control measures can be found in \citet{Bayas2017}. Each location was labelled by multiple labelers; to turn these labels into binary labels, we took the mean of the labels for each location, with a threshold of 0.5 to define a point as ``crop'' or ``non-crop.''

\subsubsection{Active Labeling}
To generate training labels, experts in photointerpretation of agricultural land cover (authors Becker-Reshef, Barker, Hosseini, Kerner, Munshell, and Paliyam) drew polygons over pixels with labels of ``crop'' or ``non-crop.'' We used a combination of basemaps for photointerpretation: a SkySat 72cm composite of January-March 2019 \citep{Planet}, PlanetScope composites of April-July and July-October 2019 \citep{Planet}, and Google Earth Pro basemaps (comprised primarily of high-resolution Maxar images). We used QGIS and Google Earth Pro to draw polygons over basemaps. The locations of labels were chosen to be distributed across different agroecological zones and non-crop land cover types observed across Togo. The central pixel within each polygon was used for training. We obtained 394 crop and 194 non-crop labels for the initial training dataset to supplement the 43 examples within Togo in the \citet{GeoWiki} dataset (631 examples total). After training the model with the hand-labeled and \citet{GeoWiki} dataset, we collected additional labeled examples (294 crop and 394 non-crop) by analyzing the predicted map and identifying areas of model confusion (e.g., transient vegetation on shorelines that the model might confuse with crops). This allowed us to maximize the utility of each training example, as it focused on places where the model was weakest rather than adding redundancy to the dataset. This resulted in a final training dataset of 1,319 hand-labeled examples, visualized spatially in Figure \ref{fig:labels}a.

\subsection{EO Data}
We used Copernicus Sentinel-2 surface reflectance (Level 2A) observations that corresponded to label locations as input to the model. We exported 160 m $\times$ 160 m patches for each labeled pixel using Google Earth Engine (GEE) \citep{GORELICK201718}, then extracted the closest pixel within the patch to the label location. To construct a cloud-free timeseries representation of each pixel, we used the algorithm from \citet{CloudFreeSentinel}, which finds the least cloudy pixel within a defined time period. We used 12 30-day time periods to construct a monthly timeseries spanning 360 days for each labeled pixel. We used all spectral bands except for B1 (coastal aerosol) and B10 (cirrus SWIR), and added the normalized difference vegetation index (NDVI), computed as $\text{NDVI} = \frac{\text{B08}-\text{B04}}{\text{B08}+\text{B04}}$. All bands (which range from 10 to 60m resolution) were upsampled to 10 m during GEE export. We used observations acquired between March of Year N and March of Year N+1 where N is the year the labels were created for. Thus for \citet{GeoWiki}, we used observations acquired March 2017-March 2018 and for our hand-labeled dataset we used observations acquired March 2019-March 2020.

\subsection{Model} \label{section:model}
We trained a one layer multi-headed LSTM model (Figure \ref{fig:lstm}) to predict whether a pixel contained cropland (1) or not (0). The input was expressed as a timeseries $\bm{X} = \{\bm{x}_{1}, \bm{x}_{2}, ... \bm{x}_{12}\}$, where each timestep $\bm{x}_{i} \in \mathbb{R}^{1\times14}$ consisted of the least-cloudy composite over a 30-day period of 11 optical bands (all bands except B01 and B10) plus NDVI. The LSTM hidden layer had 64 units. A classifier and a sigmoid activation were used on the final hidden output to map the hidden layer output (logits) to a value that can be interpreted as the posterior probability of crop presence in the pixel.

We performed a grid search that optimized AUC ROC on the validation set to choose values for two hyperparameters: number of linear layers in the classifier (1 or 2) and whether or not a dropout of 0.2 should be applied between each timestep in the LSTM.

\begin{figure}
  \centering
  \includegraphics[width=\linewidth]{diagrams/model.png}
  \caption{A multi-headed LSTM based model for pixel-wise crop identification. The model has two classification layers: a global linear classifier to classify instances outside Togo, and a local linear classifier to classify instances within Togo.}
  \label{fig:lstm}
\end{figure}

To allow the model to focus on examples within Togo while also learning from the global \citet{GeoWiki} dataset, we employed a multi-headed LSTM model, where one classifier was used for all examples within Togo (local instances) and another was used for all other examples (global instances). All instances were shuffled during training: a batch could therefore contain a combination of global and local data instances. To combine the results from both classifiers, we used the following loss function to train the model:
\begin{equation}
\mathcal{L} = \frac{W}{\alpha}\mathcal{L}_{\textrm{global}} + \mathcal{L}_{\textrm{local}}
\end{equation}
Where $\mathcal{L}_{\textrm{global}}$ and $\mathcal{L}_{\textrm{local}}$ are the binary cross-entropy losses for the global and local data instances respectively, and $\alpha$ is a weighting parameter. $W$ is a batch-specific value, used to weight the ratio of global and local instances in a batch:
\begin{equation}
    W = \frac{\mbox{Number of global instances in batch}}{\mbox{Number of local instances in batch}}
\end{equation}
In our experiments, we used $\alpha = 10$. We used the Adam optimizer with default parameters to tune weights during training \citep{Kingma2014}. To determine when to stop training, we used early stopping with a patience of 10 using a validation set consisting of 20\% of the combined \citet{GeoWiki} and hand-labeled Togo datasets. We implemented the model and experiments using PyTorch \citep{pytorch} in python. All code is publicly available at \url{https://github.com/nasaharvest/togo-crop-mask}.

\begin{table*}
  \caption{Results for each method evaluated on the test set. The best results are in bold.}
  \label{tab:experimental-results}
  \begin{tabular}{cccccc}
    \toprule
    Description &Accuracy &AUC &Precision &Recall &F1\\
    \midrule
    \citet{CopernicusLandCover} (Copernicus 100m) & 0.67 & 0.74 & 0.54 & 0.24 & 0.33\\
    Single-headed LSTM (hand-labeled only) & 0.71 & 0.90 & 0.59 & \textbf{0.93} & 0.72\\
    \citet{GFSAD} (GFSAD 30m) & 0.74 & - & 0.62 & 0.62 & 0.62 \\
    Support Vector Machine (hand-labeled only) & 0.79 & 0.88 & 0.65 & 0.84 & 0.73 \\
    Random Forest (hand-labeled only) & 0.81 & 0.90 & 0.66 & 0.88 & \textbf{0.76} \\
    Multi-headed LSTM (hand-labeled + \citet{GeoWiki}) & \textbf{0.83} & \textbf{0.91} & \textbf{0.81} & 0.68 & 0.74\\
    \bottomrule
  \end{tabular}
  \label{experimental-setup-results}
\end{table*}

\section{Results}
\label{sec:results}

We selected points for a test dataset by randomly sampling 350 points within the Togo country boundaries, constrained by a minimum distance of 50 m between samples (Figure \ref{fig:labels}b). Four experts (authors Barker, Kerner, Munshell, and Nakalembe) labeled each example. Labels were determined by majority vote in order to mitigate label noise in the absence of groundtruth validation, and ties were discarded. The consensus test dataset contained 106 crop and 200 non-crop examples (306 total test examples).

We compared two experimental setups to evaluate the contribution of the \citet{GeoWiki} dataset:
\begin{itemize}

    \item A single-headed LSTM trained on only the hand-labeled Togo data
    \item A multi-headed LSTM trained on both the hand-labeled Togo and \citet{GeoWiki} data
\end{itemize}

In both cases, we used a model with hyperparameters determined by the grid search described in \ref{section:model}. For the single headed model, this consisted of a single classification layer in the classifier \textit{without} dropout in the LSTM. For the multi-headed model, this consisted of two classification layers in the classifier \textit{with} dropout in the LSTM.

\begin{figure}
  \centering
  \includegraphics[width=0.9\linewidth]{diagrams/togo_full.pdf}
  \caption{Cropland maps of Togo generated using our method, \citet{CopernicusLandCover}, and \citet{GFSAD}. White boxes indicate locations of insets shown in Figure \ref{fig:comparisons}.}
  \label{fig:fullmap}
\end{figure}

\begin{figure*}
  \centering
  \includegraphics[width=0.9\linewidth]{diagrams/comparisons.pdf}
  \caption{Qualitative comparison of cropland maps produced using our multi-headed and single-headed LSTM methods, \citet{CopernicusLandCover}, and \citet{GFSAD} for three example regions, with SkySat 72 cm basemap (April-March composite) for reference. Locations of each inset indicated by white boxes in Figure \ref{fig:fullmap}.}
  \label{fig:comparisons}
\end{figure*}

We compared our results to two publicly-available cropland maps of Togo. One map is Global Food Security-support Analysis Data (GFSAD) cropland extent map of Africa \citep{GFSAD}. This map has 30 m resolution and was produced using two pixel-based supervised classifiers, a Random Forest and a Support Vector Machine (SVM), and one object-based classifier, a Recursive Hierarchical Image Segmentation classifier. The GFSAD map does not provide cropland probabilities, but a discrete class assigned to each pixel of 0 (water), 1 (non-cropland), or 2 (cropland). The most recent GFSAD product is from 2015. Prior to our map, the GFSAD map is to our knowledge the highest resolution cropland map publicly available for Togo. Figure \ref{fig:fullmap} (right) shows the GFSAD land cover map for Togo. The second map we compared our results to is the Copernicus Land Cover cropland map \citep{CopernicusLandCover}. The Copernicus map specifies the fraction of each pixel covered by crops (using values 0-100) at 100 m resolution and was produced using a combination of random forests and expert rules with satellite data inputs from the PROBA-V sensor. The most recent Copernicus map is for 2018. To enable quantitative comparison with our results, we treated crop fraction as analogous to crop probability and computed all performance metrics using our test dataset (Table \ref{tab:experimental-results}); however, we emphasize that this may not be how the map was intended to be used and qualitative comparisons may be more faithful to the intended use of the \citet{CopernicusLandCover} map. Figure \ref{fig:fullmap} (center) shows the Copernicus crop fraction map for Togo. A direct comparison between our method and \citet{GFSAD} and \citet{CopernicusLandCover} would require running all three methods on the same data (i.e., 10m Sentinel-2 observations paired with our labeled examples). Since neither the code for these methods nor sufficient information for reproducing them faithfully are available, we additionally performed two baseline comparisons to Random Forest and SVM classifiers which are the basis for the \citet{GFSAD} and \citet{CopernicusLandCover} methods.
For the Random Forest, we used 100 estimators with no maximum depth (nodes were expanded until all leaves were pure, or until the leaves contained fewer than 2 samples). For the SVM, we used a radial basis function kernel. We implemented both methods using Scikit-learn in python \citep{scikit-learn}.

To assess performance, we computed the area under the receiver operating characteristic (ROC) curve (AUC score), accuracy (0.5 threshold on posterior probabilities), precision (user's accuracy), recall (producer's accuracy), and F1-score (harmonic mean between precision and recall). We did not compute AUC for GFSAD, which only provides the discrete predicted class. The results are reported in Table \ref{tab:experimental-results}. Figure \ref{fig:fullmap} (left) shows our predicted cropland map.

To qualitatively compare our method compared with \citet{GFSAD} and \citet{CopernicusLandCover}, we show each map compared to the SkySat 72cm basemap (true color) for several example locations in Figure \ref{fig:comparisons}. Each region in Figure \ref{fig:comparisons} contains primarily cropland, though the level of vegetation in non-crop areas varies in each image since Togo's climate ranges from dry savanna in the north to tropical in the south. While all maps detect most fields in the second row, the \citet{GFSAD} and \citet{CopernicusLandCover} maps failed to detect most of the fields in the first and third rows.

Our resulting map of cropland probabilities at 10 m resolution as well as all hand-labeled training and testing labels are publicly available at: \url{https://doi.org/10.5281/zenodo.3836629}.

\section{Discussion}
 The multi-headed model trained on the combined hand-labeled and \citet{GeoWiki} dataset had the highest score for most performance metrics. While AUC score is a measure of model performance across a range of possible thresholds on posterior probability, the decision boundary when using the sigmoid function as the output activation layer is 0.5. While in practice a different threshold could be chosen, it is not clear how to construct a large and diverse enough validation set that could be used to reliably select an alternate threshold over very large, heterogeneous regions (e.g., the country of Togo). Thus, the improved accuracy and precision of the multi-headed model (using a threshold of 0.5) - in addition to its improved AUC score - made it preferable. Additionally, since the multi-headed model included training data from prior years (2017) and global instances of cropland \citep{GeoWiki}, we expect that this model would generalize better to potential domain shift in future years than the single-headed model trained with only one year of local data. Figure \ref{fig:comparisons} shows that the multi-headed LSTM and single-headed LSTM have similar patterns of detection, but overall the multi-headed LSTM has less extreme probabilities (values closer to 0.5 than 0 or 1). This suggests that the global examples used to train the multi-headed LSTM act like a regularizer for the classifier by moderating the confidence of predictions and preventing overfitting.

 We found that the \citet{CopernicusLandCover} map had the lowest performance for the test dataset, with an F1 score significantly lower than other methods, though we note that these results could underestimate the performance because the map reports the fraction of each pixel covered by crops, not probability, in each pixel. This lower performance may, in part, be due to the coarse spatial resolution of EO data used to produce the map and the growth in cropland area that has likely occurred since the map was produced in 2015. Similarly, the lower performance of the \citet{GFSAD} could also, in part, be due to lower spatial resolution of input data (30 m) and changes in cropland area that might have occurred between 2018 (when the most recent \citet{GFSAD} map was produced) and 2019.

 \subsection{Sensitivity to number of local examples}

\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{diagrams/auc_acc_mean_3.png}
    \caption{A plot of AUC ROC and accuracy scores of the multi-headed model as the size of the hand-labeled dataset size increases. Each score is calculated from a mean of 3 runs with different random seeds.}
    \label{fig:auc_from_size}
\end{figure}

To assess the sensitivity of the model to the frequency and number of local (Togo) examples in the training dataset and provide a guide for future labeling campaigns, we measured the performance of the multi-headed model on the test set as a function of the number of hand-labeled examples used in training (Figure \ref{fig:auc_from_size}). Although the performance of the model increased as the size of the hand-labeled dataset increased, the model performs well even with very little labeled data. This suggests this technique may be applied to other data-sparse regions with a comparable or fewer number of labels as we created for this study.

\subsection{Consensus in photointerpreted labels}

\begin{figure}
  \centering
  \includegraphics[width=0.7\linewidth]{diagrams/consensus.pdf}
  \caption{Tetrahedron showing number of examples (of 350 total) for which each pair of labelers (indicated by initials) agreed.}
  \label{fig:consensus}
\end{figure}

As described in Section \ref{sec:results}, we created a test dataset by randomly sampling 350 points within Togo and having four experts label each point as ``crop'' or ``non-crop'' using photointerpretation of high-resolution satellite imagery. For each point, we chose the label assigned by the majority of labelers (ties were discarded). Figure \ref{fig:consensus} shows the number of examples for which each pair of labelers chose the same label. Agreement between pairs of labelers ranged from 229-278 (65\%-79\%) of the total points and all labelers agreed unanimously on 181 (52\%) of the points. This level of consensus was surprisingly low, and while our photointerpretation-based dataset is a high-quality surrogate for ground-truth labels, it underscores the difficulty for even experts to determine land cover types from high-resolution satellite imagery, particularly in developing countries. This should be considered in the design of future photointerpretation-based labeling efforts, especially crowdsourced efforts in which labelers may not be experts. We plan to conduct a future study to assess agreement between labels determined by experts for the same points from photointerpretation vs. field observation to better characterize the level of accuracy to expect from photointerpreted labels for different regions and land cover types.

\subsection{Lessons learned for rapid response}
The urgent need for a cropland map of Togo that captured smallholder farms (usually $<$1 ha) required us to leverage our research methods to rapidly deliver a result for operational use in less than 10 days. During this process, we gathered several ``lessons learned'' that we share to help prepare other researchers to leverage their methods for rapid humanitarian response in the future, as well as encourage community practices that facilitate rapid response.
\begin{itemize}
\item Good software engineering practices (such as developing modular and well-documented code) is often prioritized in operational rather than research settings. However, developing research methods motivated by their use for rapid humanitarian response in the future requires integration of operational practices into research workflows. We found that developing research methods using operational software engineering practices enabled us to quickly execute experiments and apply the method to large geographic areas when the need for rapid response arose.
\item Use of cloud services (in this case AWS) allowed us to easily scale our model (particularly during country-wide inference) and enable faster computation in a short period of time. Additionally, the use of EC2 instances enabled processes to be run, monitored, and modified by a distributed team rather than a single user.
\item File transfer/download remains a bottleneck for rapid response workflows. Transferring data from storage (e.g., via Google Earth Engine, the Sentinel-2 AWS S3 bucket, or USGS Earth Explorer) to compute (e.g., AWS) constituted a significant part of the time to deliver our cropland map. There is a community need for better tools to enable searching, querying, and subsetting of geospatial data on cloud platforms.
\item Given the complexity of deep learning architectures such as Inception-v3 \citep{Szegedy2016} or ResNet \citep{He2016} designed for image classification, it might seem that classification of land cover or crops in remote sensing data---which often exhibit high intraclass variance and low inter-class variance---would also require complex architectures. However, we found that a relatively simple, shallow architecture achieved good results for this task, and has the benefit of being faster in training/inference, requiring fewer examples, and being less susceptible to overfitting than more complex architectures.
\item The GFSAD \citep{GFSAD} and Copernicus \citep{CopernicusLandCover} cropland maps were easy to find and download and were well-documented, which enabled us to quickly benchmark our map against existing maps. To give others this same benefit with our cropland map, we have made our cropland map, code, and training labels publicly available.
\end{itemize}

\section{Conclusion}
Cropland maps are critical inputs for decision-makers to rapidly design and enact policies during humanitarian crises, e.g., to deliver aid or boost productivity for farmers, yet unavailability of high-resolution, up-to-date cropland maps in most countries globally precludes their use for rapid response. The development of cropland maps, in turn, is limited by the availability of ground data that can be used for training and validating machine learning classifiers.

We present a method for cropland classification in regions with little to no available ground data that uses a multi-headed LSTM network to learn global and local features for identifying cropland based on multispectral time-series Sentinel-2 observations. We used this method to create a cropland map of Togo at 10 m resolution in fewer than 10 days to assist the Togolese government in decision-making about aid distribution to farmers during the COVID-19 crisis. This demonstrated the successful transition of machine learning research to operational use in rapid response for a real humanitarian crisis, and we provided lessons learned from this experience to facilitate this transition for other researchers in the future. Additionally, our approach gives an example of how to leverage existing readily-available labeled datasets and reduce the number of data points to be collected in field campaigns, which can be costly as well as infeasible due to regional insecurity or travel restrictions (e.g., during the COVID-19 pandemic). In future work, we plan to investigate leveraging our multi-headed LSTM approach for crop type classification, subsetting the Geo-Wiki dataset based on geographical or agroecological zones, and including Sentinel-1 synethic aperture radar data in the model input to mitigate the effects of clouds during the growing season. Additionally, we will apply our method for cropland mapping in other data-sparse African countries.

\begin{acks}
We would like to thank Tanya Harrison, Zara Khan, and Charlie Candy from Planet, Inc., for making SkySat and PlanetScope basemaps available for this study. We are also grateful to Noel Gorelick from Google Earth Engine for helping to accelerate our Google Earth Engine processing workflow in order to meet this rapid response timeline. Finally, we would like to thank Steffen Fritz, Linda See, and the rest of the Geo-Wiki team for providing the publicly-available, global, crowdsourced cropland dataset critical to our method.
\end{acks}

\bibliographystyle{ACM-Reference-Format}

\section{Introduction}

Satellite Earth observation (EO) data can provide critical, relevant, and timely information in support of agricultural monitoring and food security \citep{BeckerReshef2019,NationalAcademies2019}. EO can enable crop yield assessment and forecasting \citep{Becker2010,Franch2015,Franch2017,Lobell2013,Porter2017}, inform agricultural commodities markets, enable early warning and mitigation of impending crop shortages, and inform decisions about subsidies in times of crisis such as extreme drought \citep{Dorward2011} or outbreaks of disease or pests. However, the lack of up-to-date, high-resolution cropland and crop type maps is a critical barrier to agricultural and food security assessments and operational rapid response to humanitarian crises. These data are especially lacking in regions dominated by smallholder agriculture that are most vulnerable to food insecurity (e.g., sub-saharan Africa).

Machine learning techniques have commonly been used to identify cropland in EO data. The majority of studies use tree-based classifiers (primarily random forests or decision trees) or neural network/deep learning methods (primarily recurrent or convolutional neural networks). However, prior methods have typically been applied to small and/or spatially homogeneous areas \citep{RNN_Camargue, RNN_winter, Skakun2019}, require large training datasets (hundreds of thousands of labels) \citep{breizhcrops2020,Wang2020,lstm_crop_identification, RF_france}, or use EO inputs with insufficient spatial resolution for reliably detecting smallholder farms which often have field sizes $<$1 ha \citep{GFSAD,Teluguntla2018,Oliphant2019,CopernicusLandCover,Samasse2020,Waldner2017}. In many regions, particularly developing countries and regions dominated by smallholder farming, the lack of readily accessible ground truth data is a major barrier for training and validating machine learning techniques for cropland classification \citep{Coutu2020}.

In this paper, we demonstrate a novel method for rapidly generating cropland maps for areas with little to no available labeled data over a large heterogeneous area. This is achieved by leveraging global and local datasets. The global dataset contains tens of thousands of crowdsourced labels from diverse geographies \citep{GeoWiki} to learn general features for identifying cropland globally. These global examples supplement a much smaller dataset of curated, local examples (manually labeled using photointerpretation of high-resolution satellite imagery) to learn features specific to cropland in the target mapping region. Combined, these datasets provide sufficient examples to train a Long Short Term Memory network (LSTM) \citep{LSTM} to predict the presence of cropland in each pixel given a monthly time-series of Sentinel-2 multispectral observations.

The COVID-19 global pandemic has severely impacted the global food system, and the threat of food shortages and rising food prices has required countries at risk of food insecurity to rapidly design and enact aid programs in response. In the West African country of Togo, the government announced a program to boost national food production in response to the COVID-19 crisis by distributing aid to farmers, but found that high-resolution spatial information about the distribution of farms across Togo (which primarily consist of smallholder farms under 1 ha) was a critical piece of missing information for designing this program. Using our method, we were able to generate an accurate high-resolution (10 m) map of cropland in Togo for 2019 in under 10 days to respond to the Togolese government's immediate need for information about the distribution of farms across the country. We showed that our model outperforms existing publicly-available cropland maps for Togo \citep{GFSAD, CopernicusLandCover} while also being higher resolution and more recent. We have delivered this product to the Togolese government where it is being used to help inform decision-making in response to the COVID-19 pandemic, thus demonstrating a successful transition of machine learning research to operational rapid response for a real humanitarian crisis. We have made our map as well as all data and code used in this study publicly available to facilitate future research and operational systems in data-sparse regions.
\begin{figure}
  \centering
  \includegraphics[width=0.75\linewidth]{diagrams/label_distribution.pdf}
  \caption{Distribution of labeled examples in training (a) and test (b) datasets.}
  \label{fig:labels}
\end{figure}

\section{Methods}

\subsection{Labeled Data}
\subsubsection{Crowdsourced Labels}
We leveraged a large (35,866 labels) crowd-sourced dataset with labels of ``crop'' or ``non-crop'' from diverse, globally-distributed locations \citep{GeoWiki} to help the model learn features that are useful for detecting a wide range of cropland globally. Details about sample selection, labeling procedures, and quality assessment/control measures can be found in \citet{Bayas2017}. Each location was labelled by multiple labelers; to turn these labels into binary labels, we took the mean of the labels for each location, with a threshold of 0.5 to define a point as ``crop'' or ``non-crop.''

\subsubsection{Active Labeling}
To generate training labels, experts in photointerpretation of agricultural land cover (authors Becker-Reshef, Barker, Hosseini, Kerner, Munshell, and Paliyam) drew polygons over pixels with labels of ``crop'' or ``non-crop.'' We used a combination of basemaps for photointerpretation: a SkySat 72cm composite of January-March 2019 \citep{Planet}, PlanetScope composites of April-July and July-October 2019 \citep{Planet}, and Google Earth Pro basemaps (comprised primarily of high-resolution Maxar images). We used QGIS and Google Earth Pro to draw polygons over basemaps. The locations of labels were chosen to be distributed across different agroecological zones and non-crop land cover types observed across Togo. The central pixel within each polygon was used for training. We obtained 394 crop and 194 non-crop labels for the initial training dataset to supplement the 43 examples within Togo in the \citet{GeoWiki} dataset (631 examples total). After training the model with the hand-labeled and \citet{GeoWiki} dataset, we collected additional labeled examples (294 crop and 394 non-crop) by analyzing the predicted map and identifying areas of model confusion (e.g., transient vegetation on shorelines that the model might confuse with crops). This allowed us to maximize the utility of each training example, as it focused on places where the model was weakest rather than adding redundancy to the dataset. This resulted in a final training dataset of 1,319 hand-labeled examples, visualized spatially in Figure \ref{fig:labels}a.

\subsection{EO Data}
We used Copernicus Sentinel-2 surface reflectance (Level 2A) observations that corresponded to label locations as input to the model. We exported 160 m $\times$ 160 m patches for each labeled pixel using Google Earth Engine (GEE) \citep{GORELICK201718}, then extracted the closest pixel within the patch to the label location. To construct a cloud-free timeseries representation of each pixel, we used the algorithm from \citet{CloudFreeSentinel}, which finds the least cloudy pixel within a defined time period. We used 12 30-day time periods to construct a monthly timeseries spanning 360 days for each labeled pixel. We used all spectral bands except for B1 (coastal aerosol) and B10 (cirrus SWIR), and added the normalized difference vegetation index (NDVI), computed as $\text{NDVI} = \frac{\text{B08}-\text{B04}}{\text{B08}+\text{B04}}$. All bands (which range from 10 to 60m resolution) were upsampled to 10 m during GEE export. We used observations acquired between March of Year N and March of Year N+1 where N is the year the labels were created for. Thus for \citet{GeoWiki}, we used observations acquired March 2017-March 2018 and for our hand-labeled dataset we used observations acquired March 2019-March 2020.

\subsection{Model} \label{section:model}
We trained a one layer multi-headed LSTM model (Figure \ref{fig:lstm}) to predict whether a pixel contained cropland (1) or not (0). The input was expressed as a timeseries $\bm{X} = \{\bm{x}_{1}, \bm{x}_{2}, ... \bm{x}_{12}\}$, where each timestep $\bm{x}_{i} \in \mathbb{R}^{1\times14}$ consisted of the least-cloudy composite over a 30-day period of 11 optical bands (all bands except B01 and B10) plus NDVI. The LSTM hidden layer had 64 units. A classifier and a sigmoid activation were used on the final hidden output to map the hidden layer output (logits) to a value that can be interpreted as the posterior probability of crop presence in the pixel.

We performed a grid search that optimized AUC ROC on the validation set to choose values for two hyperparameters: number of linear layers in the classifier (1 or 2) and whether or not a dropout of 0.2 should be applied between each timestep in the LSTM.

\begin{figure}
  \centering
  \includegraphics[width=\linewidth]{diagrams/model.png}
  \caption{A multi-headed LSTM based model for pixel-wise crop identification. The model has two classification layers: a global linear classifier to classify instances outside Togo, and a local linear classifier to classify instances within Togo.}
  \label{fig:lstm}
\end{figure}

To allow the model to focus on examples within Togo while also learning from the global \citet{GeoWiki} dataset, we employed a multi-headed LSTM model, where one classifier was used for all examples within Togo (local instances) and another was used for all other examples (global instances). All instances were shuffled during training: a batch could therefore contain a combination of global and local data instances. To combine the results from both classifiers, we used the following loss function to train the model:
\begin{equation}
\mathcal{L} = \frac{W}{\alpha}\mathcal{L}_{\textrm{global}} + \mathcal{L}_{\textrm{local}}
\end{equation}
Where $\mathcal{L}_{\textrm{global}}$ and $\mathcal{L}_{\textrm{local}}$ are the binary cross-entropy losses for the global and local data instances respectively, and $\alpha$ is a weighting parameter. $W$ is a batch-specific value, used to weight the ratio of global and local instances in a batch:
\begin{equation}
    W = \frac{\mbox{Number of global instances in batch}}{\mbox{Number of local instances in batch}}
\end{equation}
In our experiments, we used $\alpha = 10$. We used the Adam optimizer with default parameters to tune weights during training \citep{Kingma2014}. To determine when to stop training, we used early stopping with a patience of 10 using a validation set consisting of 20\% of the combined \citet{GeoWiki} and hand-labeled Togo datasets. We implemented the model and experiments using PyTorch \citep{pytorch} in python. All code is publicly available at \url{https://github.com/nasaharvest/togo-crop-mask}.

\begin{table*}
  \caption{Results for each method evaluated on the test set. The best results are in bold.}
  \label{tab:experimental-results}
  \begin{tabular}{cccccc}
    \toprule
    Description &Accuracy &AUC &Precision &Recall &F1\\
    \midrule
    \citet{CopernicusLandCover} (Copernicus 100m) & 0.67 & 0.74 & 0.54 & 0.24 & 0.33\\
    Single-headed LSTM (hand-labeled only) & 0.71 & 0.90 & 0.59 & \textbf{0.93} & 0.72\\
    \citet{GFSAD} (GFSAD 30m) & 0.74 & - & 0.62 & 0.62 & 0.62 \\
    Support Vector Machine (hand-labeled only) & 0.79 & 0.88 & 0.65 & 0.84 & 0.73 \\
    Random Forest (hand-labeled only) & 0.81 & 0.90 & 0.66 & 0.88 & \textbf{0.76} \\
    Multi-headed LSTM (hand-labeled + \citet{GeoWiki}) & \textbf{0.83} & \textbf{0.91} & \textbf{0.81} & 0.68 & 0.74\\
    \bottomrule
  \end{tabular}
  \label{experimental-setup-results}
\end{table*}

\section{Results}
\label{sec:results}

We selected points for a test dataset by randomly sampling 350 points within the Togo country boundaries, constrained by a minimum distance of 50 m between samples (Figure \ref{fig:labels}b). Four experts (authors Barker, Kerner, Munshell, and Nakalembe) labeled each example. Labels were determined by majority vote in order to mitigate label noise in the absence of groundtruth validation, and ties were discarded. The consensus test dataset contained 106 crop and 200 non-crop examples (306 total test examples).

We compared two experimental setups to evaluate the contribution of the \citet{GeoWiki} dataset:
\begin{itemize}

    \item A single-headed LSTM trained on only the hand-labeled Togo data
    \item A multi-headed LSTM trained on both the hand-labeled Togo and \citet{GeoWiki} data
\end{itemize}

In both cases, we used a model with hyperparameters determined by the grid search described in \ref{section:model}. For the single headed model, this consisted of a single classification layer in the classifier \textit{without} dropout in the LSTM. For the multi-headed model, this consisted of two classification layers in the classifier \textit{with} dropout in the LSTM.

\begin{figure}
  \centering
  \includegraphics[width=0.9\linewidth]{diagrams/togo_full.pdf}
  \caption{Cropland maps of Togo generated using our method, \citet{CopernicusLandCover}, and \citet{GFSAD}. White boxes indicate locations of insets shown in Figure \ref{fig:comparisons}.}
  \label{fig:fullmap}
\end{figure}

\begin{figure*}
  \centering
  \includegraphics[width=0.9\linewidth]{diagrams/comparisons.pdf}
  \caption{Qualitative comparison of cropland maps produced using our multi-headed and single-headed LSTM methods, \citet{CopernicusLandCover}, and \citet{GFSAD} for three example regions, with SkySat 72 cm basemap (April-March composite) for reference. Locations of each inset indicated by white boxes in Figure \ref{fig:fullmap}.}
  \label{fig:comparisons}
\end{figure*}

We compared our results to two publicly-available cropland maps of Togo. One map is Global Food Security-support Analysis Data (GFSAD) cropland extent map of Africa \citep{GFSAD}. This map has 30 m resolution and was produced using two pixel-based supervised classifiers, a Random Forest and a Support Vector Machine (SVM), and one object-based classifier, a Recursive Hierarchical Image Segmentation classifier. The GFSAD map does not provide cropland probabilities, but a discrete class assigned to each pixel of 0 (water), 1 (non-cropland), or 2 (cropland). The most recent GFSAD product is from 2015. Prior to our map, the GFSAD map is to our knowledge the highest resolution cropland map publicly available for Togo. Figure \ref{fig:fullmap} (right) shows the GFSAD land cover map for Togo. The second map we compared our results to is the Copernicus Land Cover cropland map \citep{CopernicusLandCover}. The Copernicus map specifies the fraction of each pixel covered by crops (using values 0-100) at 100 m resolution and was produced using a combination of random forests and expert rules with satellite data inputs from the PROBA-V sensor. The most recent Copernicus map is for 2018. To enable quantitative comparison with our results, we treated crop fraction as analogous to crop probability and computed all performance metrics using our test dataset (Table \ref{tab:experimental-results}); however, we emphasize that this may not be how the map was intended to be used and qualitative comparisons may be more faithful to the intended use of the \citet{CopernicusLandCover} map. Figure \ref{fig:fullmap} (center) shows the Copernicus crop fraction map for Togo. A direct comparison between our method and \citet{GFSAD} and \citet{CopernicusLandCover} would require running all three methods on the same data (i.e., 10m Sentinel-2 observations paired with our labeled examples). Since neither the code for these methods nor sufficient information for reproducing them faithfully are available, we additionally performed two baseline comparisons to Random Forest and SVM classifiers which are the basis for the \citet{GFSAD} and \citet{CopernicusLandCover} methods.
For the Random Forest, we used 100 estimators with no maximum depth (nodes were expanded until all leaves were pure, or until the leaves contained fewer than 2 samples). For the SVM, we used a radial basis function kernel. We implemented both methods using Scikit-learn in python \citep{scikit-learn}.

To assess performance, we computed the area under the receiver operating characteristic (ROC) curve (AUC score), accuracy (0.5 threshold on posterior probabilities), precision (user's accuracy), recall (producer's accuracy), and F1-score (harmonic mean between precision and recall). We did not compute AUC for GFSAD, which only provides the discrete predicted class. The results are reported in Table \ref{tab:experimental-results}. Figure \ref{fig:fullmap} (left) shows our predicted cropland map.

To qualitatively compare our method compared with \citet{GFSAD} and \citet{CopernicusLandCover}, we show each map compared to the SkySat 72cm basemap (true color) for several example locations in Figure \ref{fig:comparisons}. Each region in Figure \ref{fig:comparisons} contains primarily cropland, though the level of vegetation in non-crop areas varies in each image since Togo's climate ranges from dry savanna in the north to tropical in the south. While all maps detect most fields in the second row, the \citet{GFSAD} and \citet{CopernicusLandCover} maps failed to detect most of the fields in the first and third rows.

Our resulting map of cropland probabilities at 10 m resolution as well as all hand-labeled training and testing labels are publicly available at: \url{https://doi.org/10.5281/zenodo.3836629}.

\section{Discussion}
 The multi-headed model trained on the combined hand-labeled and \citet{GeoWiki} dataset had the highest score for most performance metrics. While AUC score is a measure of model performance across a range of possible thresholds on posterior probability, the decision boundary when using the sigmoid function as the output activation layer is 0.5. While in practice a different threshold could be chosen, it is not clear how to construct a large and diverse enough validation set that could be used to reliably select an alternate threshold over very large, heterogeneous regions (e.g., the country of Togo). Thus, the improved accuracy and precision of the multi-headed model (using a threshold of 0.5) - in addition to its improved AUC score - made it preferable. Additionally, since the multi-headed model included training data from prior years (2017) and global instances of cropland \citep{GeoWiki}, we expect that this model would generalize better to potential domain shift in future years than the single-headed model trained with only one year of local data. Figure \ref{fig:comparisons} shows that the multi-headed LSTM and single-headed LSTM have similar patterns of detection, but overall the multi-headed LSTM has less extreme probabilities (values closer to 0.5 than 0 or 1). This suggests that the global examples used to train the multi-headed LSTM act like a regularizer for the classifier by moderating the confidence of predictions and preventing overfitting.

 We found that the \citet{CopernicusLandCover} map had the lowest performance for the test dataset, with an F1 score significantly lower than other methods, though we note that these results could underestimate the performance because the map reports the fraction of each pixel covered by crops, not probability, in each pixel. This lower performance may, in part, be due to the coarse spatial resolution of EO data used to produce the map and the growth in cropland area that has likely occurred since the map was produced in 2015. Similarly, the lower performance of the \citet{GFSAD} could also, in part, be due to lower spatial resolution of input data (30 m) and changes in cropland area that might have occurred between 2018 (when the most recent \citet{GFSAD} map was produced) and 2019.

 \subsection{Sensitivity to number of local examples}

\begin{figure}
    \centering
    \includegraphics[width=0.9\linewidth]{diagrams/auc_acc_mean_3.png}
    \caption{A plot of AUC ROC and accuracy scores of the multi-headed model as the size of the hand-labeled dataset size increases. Each score is calculated from a mean of 3 runs with different random seeds.}
    \label{fig:auc_from_size}
\end{figure}

To assess the sensitivity of the model to the frequency and number of local (Togo) examples in the training dataset and provide a guide for future labeling campaigns, we measured the performance of the multi-headed model on the test set as a function of the number of hand-labeled examples used in training (Figure \ref{fig:auc_from_size}). Although the performance of the model increased as the size of the hand-labeled dataset increased, the model performs well even with very little labeled data. This suggests this technique may be applied to other data-sparse regions with a comparable or fewer number of labels as we created for this study.

\subsection{Consensus in photointerpreted labels}

\begin{figure}
  \centering
  \includegraphics[width=0.7\linewidth]{diagrams/consensus.pdf}
  \caption{Tetrahedron showing number of examples (of 350 total) for which each pair of labelers (indicated by initials) agreed.}
  \label{fig:consensus}
\end{figure}

As described in Section \ref{sec:results}, we created a test dataset by randomly sampling 350 points within Togo and having four experts label each point as ``crop'' or ``non-crop'' using photointerpretation of high-resolution satellite imagery. For each point, we chose the label assigned by the majority of labelers (ties were discarded). Figure \ref{fig:consensus} shows the number of examples for which each pair of labelers chose the same label. Agreement between pairs of labelers ranged from 229-278 (65\%-79\%) of the total points and all labelers agreed unanimously on 181 (52\%) of the points. This level of consensus was surprisingly low, and while our photointerpretation-based dataset is a high-quality surrogate for ground-truth labels, it underscores the difficulty for even experts to determine land cover types from high-resolution satellite imagery, particularly in developing countries. This should be considered in the design of future photointerpretation-based labeling efforts, especially crowdsourced efforts in which labelers may not be experts. We plan to conduct a future study to assess agreement between labels determined by experts for the same points from photointerpretation vs. field observation to better characterize the level of accuracy to expect from photointerpreted labels for different regions and land cover types.

\subsection{Lessons learned for rapid response}
The urgent need for a cropland map of Togo that captured smallholder farms (usually $<$1 ha) required us to leverage our research methods to rapidly deliver a result for operational use in less than 10 days. During this process, we gathered several ``lessons learned'' that we share to help prepare other researchers to leverage their methods for rapid humanitarian response in the future, as well as encourage community practices that facilitate rapid response.
\begin{itemize}
\item Good software engineering practices (such as developing modular and well-documented code) is often prioritized in operational rather than research settings. However, developing research methods motivated by their use for rapid humanitarian response in the future requires integration of operational practices into research workflows. We found that developing research methods using operational software engineering practices enabled us to quickly execute experiments and apply the method to large geographic areas when the need for rapid response arose.
\item Use of cloud services (in this case AWS) allowed us to easily scale our model (particularly during country-wide inference) and enable faster computation in a short period of time. Additionally, the use of EC2 instances enabled processes to be run, monitored, and modified by a distributed team rather than a single user.
\item File transfer/download remains a bottleneck for rapid response workflows. Transferring data from storage (e.g., via Google Earth Engine, the Sentinel-2 AWS S3 bucket, or USGS Earth Explorer) to compute (e.g., AWS) constituted a significant part of the time to deliver our cropland map. There is a community need for better tools to enable searching, querying, and subsetting of geospatial data on cloud platforms.
\item Given the complexity of deep learning architectures such as Inception-v3 \citep{Szegedy2016} or ResNet \citep{He2016} designed for image classification, it might seem that classification of land cover or crops in remote sensing data---which often exhibit high intraclass variance and low inter-class variance---would also require complex architectures. However, we found that a relatively simple, shallow architecture achieved good results for this task, and has the benefit of being faster in training/inference, requiring fewer examples, and being less susceptible to overfitting than more complex architectures.
\item The GFSAD \citep{GFSAD} and Copernicus \citep{CopernicusLandCover} cropland maps were easy to find and download and were well-documented, which enabled us to quickly benchmark our map against existing maps. To give others this same benefit with our cropland map, we have made our cropland map, code, and training labels publicly available.
\end{itemize}

\section{Conclusion}
Cropland maps are critical inputs for decision-makers to rapidly design and enact policies during humanitarian crises, e.g., to deliver aid or boost productivity for farmers, yet unavailability of high-resolution, up-to-date cropland maps in most countries globally precludes their use for rapid response. The development of cropland maps, in turn, is limited by the availability of ground data that can be used for training and validating machine learning classifiers.

We present a method for cropland classification in regions with little to no available ground data that uses a multi-headed LSTM network to learn global and local features for identifying cropland based on multispectral time-series Sentinel-2 observations. We used this method to create a cropland map of Togo at 10 m resolution in fewer than 10 days to assist the Togolese government in decision-making about aid distribution to farmers during the COVID-19 crisis. This demonstrated the successful transition of machine learning research to operational use in rapid response for a real humanitarian crisis, and we provided lessons learned from this experience to facilitate this transition for other researchers in the future. Additionally, our approach gives an example of how to leverage existing readily-available labeled datasets and reduce the number of data points to be collected in field campaigns, which can be costly as well as infeasible due to regional insecurity or travel restrictions (e.g., during the COVID-19 pandemic). In future work, we plan to investigate leveraging our multi-headed LSTM approach for crop type classification, subsetting the Geo-Wiki dataset based on geographical or agroecological zones, and including Sentinel-1 synethic aperture radar data in the model input to mitigate the effects of clouds during the growing season. Additionally, we will apply our method for cropland mapping in other data-sparse African countries.

\begin{acks}
We would like to thank Tanya Harrison, Zara Khan, and Charlie Candy from Planet, Inc., for making SkySat and PlanetScope basemaps available for this study. We are also grateful to Noel Gorelick from Google Earth Engine for helping to accelerate our Google Earth Engine processing workflow in order to meet this rapid response timeline. Finally, we would like to thank Steffen Fritz, Linda See, and the rest of the Geo-Wiki team for providing the publicly-available, global, crowdsourced cropland dataset critical to our method.
\end{acks}

\bibliographystyle{ACM-Reference-Format}


\section{Introduction}
\input{01_introduction}

\section{Related Work}
\input{02_related_work}

\section{FHIR Transformation}
\input{03_FHIR}

\section{Patient Diagnosis Prediction}
\input{04_prediction_model}

\section{Experiments and Results}
\input{05_experiments}

\section{Conclusion, Limitations, and Future Work}
\input{06_conclusion}

\begin{comment}
\begin{acks}
\end{acks}
\end{comment}

\bibliographystyle{ACM-Reference-Format}

\subsection{FHIR}
The FHIR format is a new format to exchange EHR data. The FHIR format represents data in containers that are consistent, scalable, and hierarchical. This facilitates the exchange of data between participants in health ecosystems. Moreover, the format does not ensure semantic consistency, which allows the format to be dynamically adjusted and implemented by any health system.

Every recorded event along a patient's admission is considered as a resource in the FHIR format. These resources contain multiple attributes. The ``Medication Dispense'' resource for example contains the trade name, the generic name or the medication ingredients as attributes. Each attribute is defined by a specific data structure and type, creating the hierarchical structure of the resources. These structures and requirements are explicitly documented on the FHIR website \cite{hl7_fhir}.

Each FHIR resource type belongs to a specific thematic category, which in turn belongs to a certain level depending on how sensitive and widespread the resource type is. The levels define how consistently the resource types are defined and how tightly they are governed. Resources types from the top levels are the most widespread and support the most common health care transactions.

\subsection{MIMIC-III to FHIR Transformation}
To be able to combine the MIMIC-III dataset with other EHR datasets in future applications, we mapped it to the FHIR format. Consequently, all primary MIMIC-III tables were converted to the appropriate FHIR resource type, if a corresponding FHIR resource type existed \footnote{No correspoding FHIR resource type were found for the 'callout', 'transfers' or 'drg\_codes' tables.}, see Table \ref{table:fhir_mapping} for the corresponding mapping.

\begin{table}[!ht]
	\centering
	\footnotesize
	\begin{tabular}{|c l l|}

		\hline
		 & MIMIC-III table & FHIR Resource Type \\ [0.5ex]
		\hline
		1 & patients &  patient\\
		2& admissions & encounter \\
		3& diagnosis\_icd &  encounter\\
		4& icustay & enounter\\
		5&cptevents & claim\\
		6&noteevents & diagnosticReport\\
		7&inputevents\_cv & medicationDispense\\
		8&inputevents\_mv & medicationDispense\\
		9&prescriptions & medicationRequest\\
		10&chartevents & observation\\
		11&datetimeevents & observation\\
		12&labevents & observation\\
		13&caregivers & practitioner\\
		14&procedures\_icd & procedure\\
		15&procedureevents\_mv & procedure\\
		16&microbiology & specimen\\
		17&outputevents & specimen\\
		18&service & serviceRequest\\
		19&callout & \\
		20&transfers & \\
		21&drgcodes & \\
		\hline
	\end{tabular}
	\label{table:fhir_mapping}
	\caption{MIMIC-III Data Mapping to FHIR Resources}
\end{table}

The MIMIC-III data tables correspond to level 3 and 4 FHIR resource types. Instead of saving the FHIR resources as hierarchical containers, the resources are saved as single level containers that can easily be read into Python Pandas DataFrames.

The FHIR format allows resources to  be saved in multiple file formats such as XML, JSON and Turtle. The proposed Python package saves the FHIR resources as JSON files, as these can be easily read by the Python programming language. Each MIMIC-III table was therefore converted to a collection of FHIR resource objects, saved as a GZIP compressed JSON file.
We provide a visual representation of the admission table entries mapped to the encounter resource type with greater details in the Appendix.

The python package can be downloaded from its online repository~\footnote{https://github.com/leopold-franz/MIMIC-III\_FHIR\_Transformation}. The next steps for its use are as follows: (i) import the transform function from the mimic\_fhir\_transformation.py, and (ii) call the transform function with these inputs: \textit{input\_path} of the original MIMIC-III table CSV file, \textit{output\_path} of where the collection of FHIR resources should be saved as a JSON file. Note that by adding the '.gz' extension the function can read compressed CSV files and save compressed JSON files. The function then saves the FHIR collection as a file and returns a dataframe with the flat hierarchy FHIR resources as rows. This allows the function call to be directly incorporated into any python pre-processing step. The MIMIC-III FHIR collections can then be combined with FHIR collections from other EHR datasets and used to train robust machine learning models.
\subsection{DeepObserver}
DeepObserver is a deep learning model inspired by \cite{kemp2019baselinemimic}. In this model, all numerical observations from the MIMIC-III Chartevents table, that contains vital sign measurements and laboratory results, are used for training instead of using only a few variables.

\subsubsection{Data \& Pre-processing}
The Chartevents table contains more than 330 million measurements. These measurements belong to different observation types, for example some Chartevent entries denote arterial blood pressure measurements and others denote body temperature data. Depending on the observation type, the measurements might not be numerical values, instead  short strings or DateTime objects. As a first step of pre-processing, all values from non-numerical observation types were filtered out and  only remaining 450 observation types  with numerical values were used. Models in  \cite{kemp2019baselinemimic}, which has similar pre-processing steps, also uses categorical string values converted into numeric embeddings. Next, the time difference between each measurement and the corresponding admission discharge time is computed. Subsequently, all observations are grouped together into 4 time interval  \textit{bins} similar to \cite{kemp2019baselinemimic}. The last three bins correspond to the three consecutive 8h intervals before the discharge time point and all values before the last 24h of an admission are grouped together in the first bin, see Figure \ref{figure:deepobserver_binning}.
\begin{figure}[!ht]
	\centering
	\includegraphics[width=0.45\textwidth]{figures/deepObserver_bins.png}
	\caption{Grouping numerical observations into bins.}
	\label{figure:deepobserver_binning}
\end{figure}

Measurements of the same observation type in the same bin are averaged together. The resulting data of a single admission corresponds to  450 numeric observation types each with 4 time bin values. All values of the same type are then normalized to Z-values.

The input for the DeepObserver model corresponds to an array of size (450,4) filled with normalized numerical observations, representing a patient admission.  The model is trained on 29'590 admissions, as not every MIMIC-III patient admission has numerical entries in the Chartevents table.

DeepObserver is a supervised model that predicts multiple disease diagnoses.
Each admission representation is therefore passed to the model with corresponding disease diagnosis labels. There are two locations to record the admission diagnoses within the MIMIC-III dataset.
First location corresponds to a free-text entry in the Admission table and the second corresponds to the Diagnoses ICD table, where each admission is attributed with ICD-9-CM codes that denote the official diagnosis codes for hospital utilization in the United states. There exists 15'073 ICD-9-CM codes, of which the MIMIC-III dataset contains 6'984. Due to their high number, these ICD-9-CM codes are grouped and mapped to 281 meaningful CCS categories. "The Clinical Classifications Software (CCS) is a tool for clustering patient diagnoses and procedures into a manageable number of clinically meaningful categories." \cite{icd9_ccs_mapping}. There are 284 meaningful CCS diagnosis categories.
Each admission representation is labeled with an array of size 281. Every entry in this array is a Boolean item indicating whether the admission is diagnosed with an diagnosis code of the corresponding CCS category. A positive entry signifies that the patient is diagnosed with the corresponding disease category for the admission.

\subsubsection{Sampling} A challenge with working with medical data is that very often we face an unbalanced ratio of positive to negative labels. For example some rare CCS categories appear in less than 10 admissions of the MIMIC-III dataset. This means that at least 29'580 other admissions have a negative binary label for this CCS category. If a model uses all of the training data, it will therefore learn to predict a negative label with a high probability. One solution is to under-sample the majority class, which means that negative samples are excluded from the dataset until an equal number positive and negative samples remains, and then split the dataset into train (80\%), validation (10\%) and test (10\%) partitions. Unfortunately, an equal number of positive and negative samples cannot be reached for each label by under-sampling an extremely unbalanced multi-label dataset. Therefore an iterative train-test-split method \cite{multilabelstrat} is used that maintains the original label distribution in the train, validation and test sets. This method first calculates the label distribution in the original dataset and then calculates the label frequency needed in the desired partitions to maintain the label distribution. Next, in an iterative manner the sample with the rarest label in the original dataset is moved to the partition needing the rare label most to reach the required distribution.

\subsubsection{Model Architecture}
The DeepObserver model consists of a four layer neural network with supplementary dropout layers \cite{srivastava2014dropout} as regularization. The first layer reduces the admission representations dimension by learning patterns across the time dimension. Learning these time behaviours for each observation type is extremely important, and therefore different neural network layer architectures were evaluated for the first layer. The first architecture evaluated was a simple Fully Connected Neural Network (FCNN) layer. The second architecture tested was a Convolutional Neural Network (CNN) layer with filters of size (1,4) spanning across the time dimension. Finally, a simple Recurrent Neural Network (RNN) was evaluated, which is considered effective at learning patterns from sequences. Because this sequence is very short, no Long Short Term Memory (LSTM) or GRU (Gated Recurrent Unit) RNN architectures were considered. The second and third layer of the DeepObserver Network are all fully connected neural network layers of size 512. The last layer is a fully connected neural network layer with an output size equal to the number of CCS diagnosis categories to predict and with a sigmoid activation function that transforms the each output into label probabilities.

\subsubsection{Training}
The DeepObserver model was trained for three epochs with the Adam optimizer \cite{kingma2014adam}, a learning rate of $2e-5$ and a batch size of 32. Increasing the number of epochs did not improve the model. The hardware used for training and evaluating DeepObserver was a machine with one GPU having 11GiB of memory and 2 CPUs nodes with 40GB of RAM.

\subsection{ClinicalBERT\_Multi}
We develop an adaptation over the ClinicalBERT model \cite{clinicalbert} called ClinicalBERT\_Multi, which predicts disease diagnoses instead of readmission. The ClinicalBERT architecture is able to transform clinical notes into embeddings that can provide interpretable clinical insights. Similar to \cite{clinicalbert}, our adapted model is trained only using discharge summaries and clinical notes from the first 3 days of an admission. A model trained on notes from the first 3 days of admission can predict disease diagnoses early on in the admission, when the patient is still in the hospital. This early model is therefore evaluated on clinical notes from the first 3 and first 2 days of admission. Early predictions are exceedingly valuable in a clinical setting, as they allow patient treatments to be adjusted according to the predicted diagnoses while the patients are still in the hospital.

\subsubsection{Data Pre-processing}
\label{section:clinicalbert_preprocessing}
The clinical notes are first pre-processed, lower cased and tokenized. During pre-processing some key abbreviations are replaced (e.g., 'dr.' -> 'doctor'), and some superfluous characters are removed (e.g., the new line character '\textbackslash n').

The ClinicalBERT\_Multi model is not trained with all clinical notes that exist per admission. This differs from the DeepObserver model that does use all available numerical observations per admission. Instead, similar to the original ClinicalBERT \cite{clinicalbert}, the clinical notes dataset is split into three subsets are as follows:
\begin{enumerate}
    \item[$D_{disch}$] contains all MIMIC-III discharge summaries. Discharge summaries contain a summary of the admission, a written form of the diagnosis, a list of prescribed medications and even summaries of other clinical notes, for example a radiology report summary. Each admission has only one discharge summary.
    \item[$D_{3days}$] contains all MIMIC-III clinical notes except discharge summaries from the first three days of admission. Clinical Notes other than discharge summaries address a specific issue or service and can be more comprehensive.
    \item[$D_{2days}$] contains all MIMIC-III clinical notes except discharge summaries from the first two days of admission. This subset is only used as an additional evaluation subset of the model trained on $D_{3days}$.
\end{enumerate}

$D_{disch}$ and $D_{3days}$ are then partitioned into 80\% train, 10\% validation and 10\% test sets using the same iterative train-test-split method \cite{multilabelstrat} as DeepObserver. One ClinicalBERT\_Multi model then trains on the train set of $D_{disch}$ and another ClinicalBERT\_Multi model is trained on the train partition of $D_{3days}$.  A model trained to represent notes from the first three days of admission can also represent notes from the first two days of admission. It is therefore not necessary to train an extra model using only notes from the first two days. Consequently all clinical notes in $D_{2days}$, that also exist in the test partition of $D_{3days}$ are saved as a test partition of $D_{2days}$.The $D_{2days}$ subset is then only used as an additional test set of the model trained on $D_{3days}$. In Table \ref{table:subtasks}, information on the different subsets is shown.

All notes from $D_{3days}$ or $D_{2days}$ belonging to the same admission are concatenated together to form one long text string for each admission.

The ClinicalBERT architecture \cite{clinicalbert} can however only take in 512 inputs at a time. Therefore, long text sequences are split into multiple smaller text chunks. The text chunk strings are then converted to token sequences using the BERT tokenizer. The BERT tokenizer splits a string into word units and additionally splits rare or long word units into subword units. Furthermore, this tokenizer adds an extra token at the beginning of each sequence, which after being embedded by BERT can be used to classify the whole sequence.

Our models have a different task than the original ClinicalBERT model. The binary readmission labels are therefore transformed to diagnosis array labels, that are identical to the DeepObserver labels.

\subsubsection{Model Architecture}
The ClinicalBERT\_Multi model has the same base architecture as the original ClinicalBERT model, which consists of the standard BERT architecture with an extra fully connected neural network layer on top.

To understand BERT we must understand how attention and transformers work. Attention assigns weights to every input feature based on the features importance for the task at hand. By visualizing these weights, the learned importance of the input tokens can be shown. This ability makes the model output interpretable which is highy relevant in the clinical setting where clinicians demand models to be interpretative.  Transformer encoders then use this attention to simultaneously process every element of a sequence into a sequence of deep embedding representations. BERT is based on 12 stacked transformer encoders \cite{vaswani2017attention} with 12 attention heads. The number of attention heads in BERT defines how many attention mechanisms are used per transformer, thereby defining how many filters are learned per transformer.

The last layer after BERT has an output size equal to the number of labels and a sigmoid activation function, which transforms each output into a probability. Consequently, the last layer acts as a classification layer. It is also where the ClinicalBERT\_Multi model differs from the original ClinicalBERT model, as the output size of ClinicalBERT\_Multi's classification layer is not equal to one but to the number of diagnosis categories to predict.

BERT now takes in a text chunk token sequence $\vec{T}$ as input and returns the embeddings of every token in the sequence. As mentioned in the pre-processing step a classification token is inserted at the beginning of every sequence. The last fully connected layer with weights $W$ then uses this embedded classification token  $h_{\text{CLS}}$ output by BERT to compute a diagnosis probability $P$ for each CCS category, see equation \ref{equation:classification layer}.
\begin{equation}
P(\verb|CCS category 1| = 1 | h_{\text{CLS}}) = \sigma(Wh_{CLS})
\label{equation:classification layer}
\end{equation}
Each probability determines whether the admission will be diagnosed with a disease from the corresponding CCS category.
\begin{equation}
	ClinicalBERT(\vec{T}) = \begin{bmatrix} P(\verb|CCS category 1| = 1 | h_{\text{CLS}}) \\ P(\verb|CCS category 2| = 1 | h_{\text{CLS}}) \\ ... \\ P(\verb|CCS category 281| = 1 | h_{\text{CLS}}) \end{bmatrix}
	\label{equation:clinicalbert_input_output}
\end{equation}
Equation \ref{equation:clinicalbert_input_output} shows the ClinicalBERT\_Multi input token sequence $\vec{T}$ and text chunk probability array output.a

Multiple text chunks belong to the same admission, therefore the text chunk probabilities of the same admission and same CCS category are combined to a single admission CCS probability using Equation \ref{equation:chunk_probability_combination} from \cite{clinicalbert}, where $h_{admission}$ represents the admission representation.
\begin{equation}
P(\verb|CCS category x| | h_{admission}) = \frac{P^n_{max} + P^n_{mean}*n/c}{1 + n/c}
\label{equation:chunk_probability_combination}
\end{equation}
This calculated admission level diagnosis category probability improves on the individual disease category probabilities. $P^n_{max}$ is the maximum probability of all text chunks belonging to the same admission and is used to give additional weight to certain text chunks that have a high probability. Analogously $P^n_{mean}$ is the mean probability of all text chunks and is used to attenuate noise from a certain text chunks. $n$ represents the number of text chunks that belong to the corresponding admission. $c$ is a scaling factor that balances the weight given to the maximum probability and the mean probability. $c=2$ was used like \cite{clinicalbert}

Results presented in this paper are based on the admission level diagnosis category probabilities instead of on the text chunk diagnosis probabilities, as the admission level probabilities have shown to "consistently outperform predictions on each subsequence individually by 3-8\%"\cite{clinicalbert}. The resulting admission CCS probabilities can be converted into positive or negative predictions by setting a certain threshold for each CSS category. Any admission CCS probability over the CCS category threshold is converted to a positive prediction otherwise the admission CCS probability is converted to negative prediction. \footnote{Model performance can be evaluated without setting a threshold, which is why such a threshold was not chosen.}

\subsection{Training}
The pre-trained model from \cite{clinicalbert} was used and fine-tuned to the diagnosis prediction task for one epoch with the Adam optimizer \cite{kingma2014adam} using the standard learning rate of $2e-5$ and a batch size of 8.
The training and evaluation of the ClinicalBERT\_Multi model was done on a Google Cloud Platform Virtual Machine Instance was setup with 8 CPUs, each having 30GB of memory, 1 NVIDIA Tesla K80 GPU and a 128GB SSD harddisk.

\subsection{Experiments}
The first experiment aims to evaluate the best performing DeepObserver model. Therefore evaluating which architecture for the first layer works best. The different DeepObserver models are DeepObserver FCNN, DeepObserver CNN and DeepObserver RNN, where the second part of the name defines what architecture is used for the first layer.

The second experiment evaluates the performance of different models predicting a specific CCS category: CCS Category 98, which corresponds to Hypertension, also known as high blood pressure. This particular category was chosen as it is the most frequent CCS category in the MIMIC-III dataset, appearing 21'139 times. Results of this experiment are labelled with 'CCS Cat:98' in the task column of Table \ref{table:results}. The models evaluated in for this subtask are DeepObserver CNN, ClinicalBERT\_Multi and additionally ClinicalBERT\_Binary. ClinicalBERT\_Binary is a replica of the original ClinicalBERT model only the task is changed to predicting CCS category 98, which means that each binary readmission label is changed to a binary label defining whether the admission is diagnosed with the CCS category 98.

The goal of the third experiment is to compare the performance of the ClinicalBERT\_Binary and ClinicalBERT\_Multi when evaluated on the different data subsets $D_{disch}$, $D_{3days}$, $D_{2days}$. This experiment therefore evaluates how well the models perform when predicting the diagnosis categories at different admission time points, i.e. after 2 days of admission, after 3 days or after discharge.

The fourth and main experiment compares the predictive performance of all models predicting all diagnosis categories simultaneously. Therefore comparing ClinicalBERT\_Multi using Discharge Summaries, ClinicalBERT\_Multi using Clinical notes from the first 3 days,  ClinicalBERT\_Multi using Clinical notes from the first 2 days, DeepObserver FCNN, DeepObserver CNN and DeepObserver RNN and the SHiP model on the multi-label prediction task.

\subsection{Results}
\label{section:results}
In Table \ref{table:results}, we present the evaluation of all models using AU-ROC, AU-PR, and Recall at $precision=80\%$ (Recall@Prec80). Note that the AU-ROC score is only representative of the predictive performance of models evaluated on balanced datasets such as the ClinicalBERT\_binary and the original ClinicalBERT. Therefore, to be able to compare all models, the AU-PR scores are compared, as well as the recall@prec80 scores. All multi-label classification models are marked with $ ^{b} $ to signify that the micro averaged scores are shown.

In the case of disease diagnosis support systems for multi-morbid patients, it is important to develop models with low miss rates. High recall scores are therefore required to make sure a few diseases are missed. Low precision scores are not very problematic, as high numbers of false alarms will be disregarded by doctors with a high precision rate.

The AU-PR score of a random classifier is equal to the ratio of positive samples in the evaluation dataset. The ratio of positive CCS Cat.:98 samples in the test partition of the Diagnoses task is 0.071. Therefore, if ClinicalBERT\_Multi and DeepObserver reach AU-PR scores above 0.071 in the CCS Cat.:98 , they outperform a random classifier. Meanwhile, the ClinicalBERT\_Binary model is tested with an equal amount of positive and negative samples, due to the original ClinicalBERT model under-sampling negative samples. Consequently, the AU-PR score of ClinicalBERT\_Binary model should be above 0.5 to outperform the random classifier. Note that the AU-PR scores of all models assessed on the CCS Cat.:98 task can still be compared. The ratio of positive labels in the Diagnoses task is 0.043, which means any Micro-AU-PR score in the Diagnoses task above this value outperforms a random classifier.

\begin{table}[!ht]

	\centering
	\begin{tabular}{l c c c c}
		\hline
		Model & Dataset & \# of Pat. & \# of Adm. & Avg. chunks \# \\ [0.5ex]
		\hline

		ClinicalBERT & $D_{disch}$ & 24'742 & 29'974 & 4.89\\
		\_Binary& $D_{3days}$ & 25'103 & 30'476 & 8.54\\
		& $D_{2days}$ &  &  & 6.28\\
		ClinicalBERT & $D_{disch}$ & 27'238 & 33'684 & 4.64\\
		\_Multi& $D_{3days}$ & 27'564 & 34'152 & 7.31\\
		& $D_{2days}$ &  &  & 5.51\\
		\hline
	\end{tabular}
	\caption{\footnotesize{Details of ClinicalBERT datasets. `\# of Pat.' and `\# of Adm.' show how many patients and admissions are in the train partition of the data subsets. `Avg. chunks \#' indicates the average number of text chunks per admission in the evaluation partition.}}
	\label{table:subtasks}

\end{table}

\begin{table}[!ht]
    \footnotesize
    \centering

	\begin{tabular}{l l l c c c}
		\toprule
		\textbf{Model} & \textbf{Task} & \textbf{$t_{pred}$} & \textbf{AU-ROC}  & \textbf{AU-PR} & \textbf{Recall@} \\
		            &       &           &          &    & \textbf{Prec80} \\ [0.5ex]
		\midrule
		\textbf{SHiP $ ^{a} $} & \textit{Diagnoses} &\textit{Discharge} & \textit{0.897} & \textit{0.352} &\textit{ -} \\
		\midrule
		\textbf{DeepObserver} & Diagnoses & Discharge & 0.897$ ^{b} $ & 0.367$ ^{b} $ & 0.039$ ^{b} $\\
		\hspace{8mm} \textbf{FCNN}         &           &           &               &               &              \\

		\midrule
		\textbf{DeepObserver} & Diagnoses & Discharge & 0.866$ ^{b} $ & 0.231$ ^{b} $ & 0.000$ ^{b} $\\
		\hspace{8mm} \textbf{RNN}         &           &           &               &               &              \\

		\midrule
		\textbf{DeepObserver} & Diagnoses & Discharge & 0.90$ ^{b} $ & 0.372$ ^{b} $ & 0.039$ ^{b} $\\
		\cmidrule(lr{1em}){2-6}
		\hspace{8mm} \textbf{FCNN} & CCS Cat: 98 & Discharge & 0.699 & 0.585 & 0.001\\
		\midrule

		\textbf{ClinicalBERT\_} & CCS Cat.: 98 & Discharge & \textbf{0.873} & \textbf{0.862} & \textbf{0.803}\\
		\cmidrule(lr{1em}){3-6}
		\hspace{8mm} \textbf{Binary} &  & 3 Days & 0.751 & 0.726 & 0.219\\
		\cmidrule(lr{1em}){3-6}
		&  & 2 Days & 0.745 & 0.720 & 0.239\\
		\midrule
		\textbf{ClinicalBERT\_} & Diagnoses & Discharge & 0.919$^{b}$ & 0.408$^{b}$ & 0.08$1^{b}$\\
		\cmidrule(lr{1em}){3-6}
		\hspace{8mm} \textbf{Multi} &  & 3 Days & \textbf{0.921}$ ^{b} $ & \textbf{0.426}$ ^{b} $ & \textbf{0.117}$ ^{b} $\\
		\cmidrule(lr{1em}){3-6}
		&  & 2 Days & 0.708$ ^{b} $ & 0.202$ ^{b} $ & 0.061$ ^{b} $\\
		\midrule
		\textbf{ClinicalBERT\_} & CCS Cat.: 98 & Discharge & 0.739 & 0.626 & 0.032\\
		\cmidrule(lr{1em}){3-6}
		\hspace{8mm} \textbf{Multi} &  & 3 Days & 0.723 & 0.613 & 0.082\\
		\cmidrule(lr{1em}){3-6}
		&  & 2 Days & 0.666 & 0.558 & 0.003\\[1ex]
		\bottomrule
	\end{tabular}
	\caption{\footnotesize{Predictive qualities of deep learning for disease diagnoses. The prediction time point is represented by $t_{pred}$. The \textbf{Bold values} are the values of the best performing model per task. $ ^{a} $ The Sequential, Hierarchical, and Pretrained Model results are retrieved from \cite{kemp2019baselinemimic} paper, which diagnoses multiple ICD-9 codes. $^{b}$ indicates that the values are micro averaged evaluation values. $ ^{c} $ The ClinicalBERT results are retrieved from the \cite{clinicalbert} paper. \textit{Results in italics are retrieved from the original papers}. Any non italicized values are results from models proposed in this paper.}}
	\label{table:results}

\end{table}

Next, we summarize the results of the predictive performance of models and the results of the experiments.

All models proposed in this paper perform better than a random guess, as they all have an AU-PR score well above the AU-PR score of a random classifier. The smallest performance improvement compared to a random classifier for models evaluated on an imbalanced dataset is the ClinicalBERT\_Multi model. It predicts the disease diagnosis codes after \textbf{2} days of admission with a performance improvement of nearly 370\% (0.202/0.043 -1). The ClinicalBERT\_Binary model is the only model trained for this paper, which trains with an equal ratio of positive and negative samples. It shows significant improvement from the random classifier  of 50\%.

The results of the first experiment show that the best performing DeepObserver model is the DeepObserver CNN model with a Micro-AU-PR score of 0.372. This indicates that the CNN filters are best at capturing time dependent patterns in these circumstances. The lower performance of the FCNN model is probably due to the higher number of trainable weights (500\% more weights). The lower performance of the DeepObserver RNN model can be explained by the low number of time steps along the time dimension the RNN can learn patterns from.

For the second experiment, we compare the models performances of models predicting the CCS Cat.:98 task. The best performing model is ClinicalBERT\_Binary predicting the diagnosis of the 98th CCS category at the Discharge time point. It achieves an AU-PR score of 0.862 and an exceptionally high recall@prec80 of 0.803. The ClinicalBERT\_Multi model trained to predict all diagnosis categories at Discharge is evaluated when predicting CCS Cat.:98 and achieves an AU-PR score of 0.626. This shows that fine-tuning the ClinicalBERT architecture on a balanced single diagnosis code dataset compared to on all diagnosis categories achieves a superior performance.

Results from the third experiment show that the prediction time point, which directly influences the amount of input data, has an effect on the predictive performance of models. The AU-PR scores of ClinicalBERT, ClinicalBERT\_Binary and ClinicalBERT\_Multi models improve when evaluating clinical notes of the first 3 days compared to evaluating clinical notes from the first 2 days. The difference in input data amount can be found in the 'Avg. chunks \#' column of Table \ref{table:subtasks}. All previously mentioned models show a higher AU-PR score when making predictions after 3 days of admission compared to when making predictions after 2 days of admission. Furthermore, ClinicalBERT\_Binary and ClinicalBERT\_Multi evaluated on task 'CCS Cat.: 98' have a better AU-PR score when using the discharge diagnoses than when using all clinical notes from the first 3 days of admission. Interestingly, the ClinicalBERT\_Multi sub-model predicting diagnoses after 3 days (AU-PR = 0.426) slightly outperforms ClinicalBERT\_Multi trained and evaluated with discharge summaries(AU-PR = 0.408). This could be due to the higher amount of input data when making predictions using clinical notes of the first 3 days of admission compared to making predictions using discharge summaries. The difference in input data is seen in column 'Avg. chunks \#' of table \ref{table:subtasks}.

The results of the fourth experiment show that the model with the highest predictive performance for the Diagnoses task is the ClinicalBERT\_Multi model predicting disease diagnosis categories using notes of the first 3 days of admission. It has an AU-PR score of 0.426. It slightly outperforms the DeepObserver CNN model making predictions after discharge (AU-PR = 0.372). This performance difference can be explained by the richer information in clinical notes when compared to the information in the numerical observations. Note that the high accuracy of the DeepObserver CNN model remains a significant achievement, as the data from the DeepObserver model is not human generated and does not depend on quality of clinical notes written by caregivers. Furthermore, comparing the results of the ClinicalBERT\_Multi model predicting CCS categories after 3 days with the results of the external SHiP model (AU-PR = 0.352) shows that ClinicalBERT\_Multi has a higher AU-PR score. The SHiP model differs from the CLinicalBERT\_Multi task by the following two points: (1) it uses all notes from an admission to make diagnosis predictions, which should help the SHiP model to achieve a higher predictive score, (2) the SHiP model predicts ICD-9-CM codes, which are harder to predict as there are more of them, instead of CCS categories.

As a summary, the DeepObserver CNN, which uses numerical observations, easily outperforms a random classifier when predicting multiple diagnosis categories simultaneously. Moreover, ClinicalBERT\_Multi, which uses clinical notes, has an even higher predictive performance than DeepObserver CNN when predicting multiple diagnosis categories and is can even provide good predictions when predicting the diagnoses early on in the admission.

\subsection{ClinicalBERT Interpretation}
By visualizing the self-attention heads from BERT, the diagnoses predictions of ClinicalBERT\_Binary and ClinicalBERT\_Multi can be interpreted. The attention function takes as input a key, a query and a value vector. Using these vectors it computes a distribution over all keys for each query and then multiplies these distribution weights with the values. In the case of self-attention the queries, keys and values are each constructed by multiplying the input tokens with a set of learned weights. The attention weight distribution for a query vector $q$ of length $ d $, and a set of keys $K$ each also of the length $ d $ is therefore:

\begin{equation}
\text{Attention}(Q, K, V) = \text{softmax}\left( \dfrac{QK^\top}{\sqrt{d}}\right)V
\end{equation}
The self-attention heads therefore act as filters over the input token sequence.

\begin{figure}[!ht]
	\centering
	\includegraphics[width=0.45\textwidth,keepaspectratio]{figures/one_attention_heads_diagnosis_98_shortInput.png}
	\caption{Self-attention weights aligned to key and query tokens from  the ClinicalBERT\_Binary model trained on discharge summaries.}
	\label{figure:attention_short_input}
\end{figure}

An attention weight with a high value signifies that the interaction between the query and key token is predictive about the task at hand. In Figure \ref{figure:attention_short_input}, the weights of a single self-attention head are visualized. As you can see a the failure token is given a high weight.

Being able to interpret the results from ClinicalBERT\_Multi sub-models is important to gain an understanding of the model predictions and better informing medical practitioner's decision-making.