# Producing Plankton Classifiers that are Robust to Dataset Shift

C. Chen<EMAIL_ADDRESS>Eawag, Überlandstrasse 133, CH-8600, Switzerland
S. Kyathanahally  Eawag, Überlandstrasse 133, CH-8600, Switzerland  M. Reyes
Eawag, Überlandstrasse 133, CH-8600, Switzerland  S. Merkli  Eawag,
Überlandstrasse 133, CH-8600, Switzerland  E. Merz  Eawag, Überlandstrasse
133, CH-8600, Switzerland  E. Francazi  Eawag, Überlandstrasse 133, CH-8600,
Switzerland  M. Hoege  Eawag, Überlandstrasse 133, CH-8600, Switzerland  F.
Pomati  Eawag, Überlandstrasse 133, CH-8600, Switzerland  M. Baity-Jesi
<EMAIL_ADDRESS>Eawag, Überlandstrasse 133, CH-8600, Switzerland

###### Abstract

Modern plankton high-throughput monitoring relies on deep learning classifiers
for species recognition in water ecosystems. Despite satisfactory nominal
performances, a significant challenge arises from Dataset Shift, which causes
performances to drop during deployment. In our study, we integrate the ZooLake
dataset, which consists of dark-field images of lake plankton [1], with
manually-annotated images from 10 independent days of deployment, serving as
test cells to benchmark Out-Of-Dataset (OOD) performances. Our analysis
reveals instances where classifiers, initially performing well in In-Dataset
conditions, encounter notable failures in practical scenarios. For example, a
MobileNet with a 92% nominal test accuracy shows a 77% OOD accuracy. We
systematically investigate conditions leading to OOD performance drops and
propose a preemptive assessment method to identify potential pitfalls when
classifying new data, and pinpoint features in OOD images that adversely
impact classification. We present a three-step pipeline: (i) identifying OOD
degradation compared to nominal test performance, (ii) conducting a diagnostic
analysis of degradation causes, and (iii) providing solutions. We find that
ensembles of BEiT vision transformers, with targeted augmentations addressing
OOD robustness, geometric ensembling, and rotation-based test-time
augmentation, constitute the most robust model, which we call BEsT . It
achieves an 83% OOD accuracy, with errors concentrated on container classes.
Moreover, it exhibits lower sensitivity to dataset shift, and reproduces well
the plankton abundances. Our proposed pipeline is applicable to generic
plankton classifiers, contingent on the availability of suitable test cells.
By identifying critical shortcomings and offering practical procedures to
fortify models against dataset shift, our study contributes to the development
of more reliable plankton classification technologies.

plankton, machine learning, classification, out of domain

††preprint: APS/123-QED

##  I Introduction

Plankton are a fundamental element of water ecosystems. They are responsible
for carbon and nutrient upcycling and are used in water treatment plants and
fisheries [2, 3]; they are a robust indicator of the health of ecosystems [4];
and they also are a paradigm of complex ecosystem which is studied as a model
system for understanding community composition and dynamics [5, 6, 7].
Additionally, algal blooms, which are sudden increases in the abundance of
specific phytoplankton species, can have disruptive effects on water
ecosystems. Consequently, there is an ongoing effort to develop forecasting
methods for predicting the abundances of plankton taxa at future times [8].

For these and various other reasons, high-throughput in-situ monitoring
systems have been increasingly deployed in recent years [9, 10, 11, 12, 13],
generating vast amounts of data [14]. These systems often employ underwater
cameras, producing an overwhelming volume of images that makes manual
annotation unapproachable. This challenge is addressed through the usage of
deep learning classifiers, which use state-of-the-art methods, ranging from
convolutional networks [15, 16, 17, 18, 19, 1] to vision transformers [20, 21,
22]. The results seem promising, exhibiting test accuracies and F1 scores well
above 90% [23]. As a result, ecological studies have commenced leveraging data
from these classifiers [24, 13, 25].

However, already in the early attempts at plankton classification, a mismatch
between the nominal performance and that at deployment time was noticed [26,
27]. In fact, it was pointed out that plankton systems are very likely to be
subject to dataset shift (DS) [28, 29, 12, 30]. In simple terms, DS means that
the train/test dataset can be so different from the data encountered once the
model is deployed, that the declared test performances are not meaningful [31,
32]. When the difference between testing and deployment are in terms of taxon
abundances, we talk of distributional DS. In plankton monitoring, not only the
relatives abundances can change, but also the way images from a same taxon
appear, in which case we talk of compositional DS. Factors which can induce
compositional DS include changes in the species traits, in environmental
factors such as water turbidity, and in instrument conditions like lens
cleanliness. In the presence of DS, it becomes crucial that machine learning
models not only perform well on the test set (which we call in-dataset, ID),
but also that they generalize well on out-of-dataset (OOD) images that may
come from a different data distribution than the original dataset.

One way to address this is to adopt an out-of-domain detection approach, in
which one tries to infer which images come from a different data distribution,
with the intention of abstaining from classifying those images [33]. One
approach with such an intent, applied to plankton, is based on discarding the
examples for which the classifier has a low confidence [34, 35]. Otherwise,
out-of-domain methods are often used with the aim of improving the
generalization to a target domain. There is a large number of existing
strategies for this [36], but they mostly rely on the availability of more
than one training domain,111Out-of-domain generalization usually refers to
training on $M$ distinct domains, and generalizing to a further one [36]. The
specific case of $M=1$ is usually called single-source out-of-domain
distribution [37, 38]. and on a well-defined target domain, which in our
problem we usually do not have access to.222While it is outside the scope of
this work, it is however possible to define multiple domains for plankton
classification. For example, one could take images from different imaging
systems, and each would constitute a different domain. For this reason, one
often talks about DS, where the changes in the data distribution are not
clearly defined a priori.

For this reason, recent literature on plankton classification focuses rather
on addressing DS through open-set recognition. For example, Ref. [39]
addresses the problem of varying class compositions between different
geographic regions, as well as the appearing of novel kinds of external
objects, within an open-set recognition framework.

Since the final objective of these classifiers is to infer abundances, with
the individual images having little importance, several works focus on
quantification [28, 40, 41, 42], that is, methods that optimize the estimation
of the abundances even in the presence of imperfect classifiers. Focusing on
quantification has the additional advantage that, instead of retraining a
classifier, one can simply recalculate the correction to the counts that would
give the best estimates of the real population.

Independently of the used classification methods, one key issue with DS is to
be able to assess performances on a dataset that is constantly varying. This
can be addressed through the creation of several additional annotated
datasets, called test cells, which represent deployment conditions [28, 29].

Most solutions for assessing performance under possible DS are based on the
availability and the production of large amounts of data (see App. A for a
summary of recent work on DS in plankton monitoring), which must be annotated
by expert taxonomists. While this data is available for datasets such as WHOI
[14], this is often not the case with other datasets, which are much smaller,
and there are therefore no studies of DS. Furthermore, while these methods
assess the presence of DS, they do not address its origin: What are the
characteristics of plankton images, that cause DS? While it is clear that
distributional DS appears in the abundances that are present in the ecosystem,
it is not clear what else is changing in the images, giving rise to
compositional DS. The evidence of further changes beyond the abundance
distribution is in fact indirect, by assessing that the test performance drop
at deployment stage. One would therefore like to be able to quantitatively
link model performances with DS, in order to answer questions such as: Are the
model performances directly related to DS? Which are the image features which
most influence the DS robustness? Which classes are harder/easier more robust
under DS? How well will a classifier perform on a new unlabeled dataset?
Moreover, this kind of assessment is model-dependent, and there is to our
knowledge no systematic study of the model dependence of DS. Finally, once one
observes that DS is present, solutions are needed to address it: Which models
and methods are the most efficient to address DS?

In this paper, we address DS on the freshwater plankton camera images from the
Dual Scripps Plankton Camera [12, 13], in a situation that is much more data-
scarce than the one addressed in [28].333In Ref. [28], the ID dataset
contained over 1000 per class, and is the result of a massive annotating
campaign, lasted years, which can hardly be replicated by future studies on
novel monitoring systems. We train models on the ZooLake2.0 dataset, which is
a minor extension of the ZooLake dataset [43]. We perform a typical
train/validation/test splitting, and call the corresponding performances in-
dataset (ID). We then simulate the deployment of the models in the field, by
further manually annotating 10 out-of-dataset (OOD) test cells, each
corresponding to a separate different day of sampling, in a different time of
the year. An alike OOD scheme was already used in the past [28, 29, 30]. With
this setup we conduct our work in three steps as shown in Fig. 1:

Figure 1: Schematic flowchart of our pipeline to address dataset shift.
Identification: Models are usually trained, validated and tested on in-dataset
(ID) data, which, due to dataset shift (DS), differs from the data found at
deployment. By using 10 out-of-dataset (OOD) test cells, we simulate the
performance at deployment, which is lower than the nominal ID test
performance. Characterization: We quantify dataset shift and the sensitivity
of our models to dataset shift. Cure: We study the effect of several modeling
choices on combating DS.

Identify dataset shift

In medium-size lakes, variations in the observed taxa and debris can be argued
to be lower, and events such as the appearing of new taxa are rare. We check
for dataset shift in the most stringent setting (no clear change in conditions
of the instruments, lake conditions, new species, and so on), and still find
that it has a very strong influence on plankton classification. For example,
during regular deployment, the classifier performance drops by between 10% and
20% with respect to the ID test accuracy. We find that while distributional DS
is present, the performance drop must be attributed mainly to compositional
DS. We quantify the two kinds of shifts and propose a way to re-assess the
performance of a classifier on new data, based on the size of these shifts and
how they affect different classes (since some classes suffer strongly, while
others barely do).

Characterize dataset shift

Compositional shift is due to variations between ID and OOD, but it is not
clear what these variations are. We identify these variations by studying
which features of the images mostly affect classification. Inferring whether
there are specific features of the images which are making the OOD
classification harder, allows us to specifically act on those.

Cure dataset shift

We are interested in finding the models that best deal with dataset shift. We
analyze different architectures to see which is more robust with respect to
DS, and use several methods for dealing with DS. We find that RGB histogram
reweighting and adjusted-count quantification do not help improve the OOD
performance. Instead, ensembling, targeted data augmentation, test-time
augmentation, and architecture selection are helpful to address DS.

Our work follows a pipeline that can be fully reproduced for model selection
with any instrument and in any plankton ecosystem (or any generic monitoring
task). We also provide our codes at
https://github.com/cchen07/Plankiformer_OOD.

##  II Materials and Procedures

###  II.1 Dataset shift

The aim of any machine learning classification application is to have a model
$f(\mathbf{x})$ which, given an input $\mathbf{x}$ (in our case $\mathbf{x}$
is an image), can assign it the correct label $y$ with a low error rate. The
error rate is measured through a loss function $\ell(f(\mathbf{x}),y)$, which
compares the output of the model, $f(\mathbf{x})$, with the true label.
Ideally, the model needs to perform well over all possible examples it may see
at deployment stage. That is, the ultimate aim is to minimize the _population
loss_ ,

$\mathcal{L}_{\mathrm{pop}}=\int
P(\mathbf{x},y)\,\ell(f(\mathbf{x}),y)\,d\mathbf{x}\,dy\,,$ (1)

where $P(\mathbf{x},y)$ is the probability that an input $\mathbf{x}$, with
label $y$ will be found at deployment stage.

Since, while training, validating and testing the model, we do not have access
to $P(\mathbf{x},y)$, we usually resort to the _empirical loss_ ,

$\displaystyle\mathcal{L}_{\mathrm{emp}}^{\mathrm{(train)}}$
$\displaystyle=\frac{1}{M_{\mathrm{train}}}\sum_{i=1}^{M_{\mathrm{train}}}\ell(f(\mathbf{x}_{i}),y_{i})\,,$
(2) $\displaystyle\mathcal{L}_{\mathrm{emp}}^{\mathrm{(val)}}$
$\displaystyle=\frac{1}{M_{\mathrm{val}}}\sum_{i=1}^{M_{\mathrm{val}}}\ell(f(\mathbf{x}_{i}),y_{i})\,,$
(3) $\displaystyle\mathcal{L}_{\mathrm{emp}}^{\mathrm{(test)}}$
$\displaystyle=\frac{1}{M_{\mathrm{test}}}\sum_{i=1}^{M_{\mathrm{test}}}\ell(f(\mathbf{x}_{i}),y_{i})\,,$
(4)

where we emphasized that one usually splits the data into a training set of
size $M_{\mathrm{train}}$ (for finding optimal weights), validation set of
size $M_{\mathrm{val}}$ (for model and hyperparameter selection), and test set
of size $M_{\mathrm{test}}$ (for performance reporting).

The examples $\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{M}$ in the empirical loss
are respectively distributed according to
$P_{\mathrm{train}}(\mathbf{x},y),P_{\mathrm{val}}(\mathbf{x},y)$ and
$P_{\mathrm{test}}(\mathbf{x},y)$. If these are representative of the
population data distribution, that is if
$P_{\mathrm{train}}(\mathbf{x},y)=P_{\mathrm{val}}(\mathbf{x},y)=P_{\mathrm{test}}(\mathbf{x},y)=P(\mathbf{x},y)$,
then the reported test performance is representative of the population loss.

While, by enforcing a random splitting of the dataset, stratified by
abundance,444This means we enforce that the relative abundance of each class
stay stable among the splittings. When the nature of the distribution shift is
known, stratification can be used to provide robust performance assessments,
by ensuring that the same shift occurs from train to validation to test set
(see e.g. Ref. [44]).  one can obtain $P_{\mathrm{train}}(\mathbf{x},y)\approx
P_{\mathrm{val}}(\mathbf{x},y)\approx P_{\mathrm{test}}(\mathbf{x},y)$, there
are several reasons why the dataset might not be representative of
$P(\mathbf{x},y)$:

High dimensionality

The inputs $\mathbf{x}$ are very high-dimensional, therefore, and to
appropriately reproduce $P(\mathbf{x},y)$, it needs to be densely sampled
within the test set.

Biased sampling

Often, dataset collection (hence, test set creation) does not faithfully
reproduce the conditions found at deployment stage. For example, there can be
an augmented effort to find examples belonging to minority classes, junk
images can be discarded, or also, images where the human annotator is unsure
may be discarded.

Non-stationarity

The function $P(\mathbf{x},y)$ represents the images that are found at
deployment stage, but these may vary in time for reasons that go from
environmental changes and instrument degradation, to the community dynamics
itself. Instead, the test set is collected once and for all, so it will
eventually not be fully representative of what is found at deployment stage.

Since we are interested in how well the test performances represent those that
one finds at deployment time, we will mostly focus on the similarity between
$P_{\mathrm{test}}(\mathbf{x},y)$ and $P(\mathbf{x},y)$, the first being ID,
and the second OOD.

Since, in plankton monitoring, $P(\mathbf{x},y)$ varies over time, it was
suggested to build $P_{\mathrm{test}}(\mathbf{x},y)$ as an unbiased union of
separate sampling sessions [29]. While this still cannot reproduce the full
variability of plankton populations, it at least reduces distribution-shift
effects due to sampling bias. Likewise, one could suggest that, if one wants
to add a specific taxon to a dataset, they should not only add images related
to the target taxon, but a whole sampling period containing that taxon, in
order not to influence the distribution of labels found in the field. This
however (i) does not address compositional shift and (ii) only partially
addresses distributional shift, since it influences the frequency of
appearance of certain communities with respect to others, and in any case it
is unlikely that any two far-enough sampling days have a similar $P(y)$.
Furthermore, as we will see, the $P(y)$ is only a smaller part of the domain
shift problem.

##### Distributional and compositional dataset shift

We can use the definition of conditional probability (this is related to
Bayes’ formula),

$P(\mathbf{x},y)=P(\mathbf{x}|y)P(y)\,,$ (5)

to elicit the sources of DS in the population loss from Eq. (1),

$\mathcal{L}_{\mathrm{pop}}=\int
P(\mathbf{x}|y)P(y)\,\ell(f(\mathbf{x}),y)\,d\mathbf{x}\,dy\,.$ (6)

In Eq. (6), we see that the value of the population loss changes if either
$P(y)$ or $P(\mathbf{x}|y)$ are different from their equivalents in the test
set, $P^{\mathrm{(test)}}(y)$ and $P^{\mathrm{(test)}}(\mathbf{x}|y)$. So, if
either of these two quantities in the test set does not reflect what is found
in the field, the model performance will be different than expected. If
$P(y)\neq P^{\mathrm{(test)}}(y)$, we talk of distributional DS. If
$P(\mathbf{x}|y)\neq P^{\mathrm{(test)}}(\mathbf{x}|y)$, we have compositional
DS.

In plankton classification, $P(y)$ represents the abundance of each taxonomic
unit. Instead, $P(\mathbf{x}|y)$ is the probability that, given a class $y$,
it has an appearance $\mathbf{x}$. In practice, it represents the fact that a
specific taxon (or auxiliary class) can appear in different ways. For example,
there can be morphological differences within the same taxonomic groups (due
e.g. to nutrient scarcity), more turbidity in the water, or varying light
conditions. In all those cases, the same class will have a different
appearance.

Figure 2: (a): Distribution of class abundances in the ZooLake2.0 and OOD
dataset. (b): Dates of collection of the images in the ZooLake2.0 dataset and
in the test cells.

###  II.2 Data

####  II.2.1 ZooLake2.0 dataset

We use images from the Dual-magnification Scripps Plankton Camera (DSPC)
deployed in lake Greifensee, in Switzerland [13]. These are color images of
varying size with a black background, taken at three meters depth, of objects
which range from 0.1mm to 10mm. These are images from the 0.5x magnification
of the DSPC, which mainly targets zooplankton and large phytoplankton
colonies. The DSPC camera in the lake Greifensee takes one picture every
second for 10 minutes every hour. To avoid repeated imaging of the same
organism, we subset the sampling frequency to one picture every 6 seconds
[13]. Each picture includes several single plankton organisms, small images
are cut out of large picture. The dataset is an extension of the ZooLake
dataset [43], which we integrated with: (i) the addition of a larger number of
junk images, which better approaches the ratio of junk images (e.g. unknown or
dirt categories) that is found in typical monitoring days;555The reason why we
felt the need to extend the need to integrate ZooLake with new images is that
we did not want our results to trivially descend from a distribution which
under-represents junk classes.  (ii) the discard of 4 mislabeled
examples.666Two images of dinobryon, one of asterionella, and one of rotifers
are discarded because of mislabelling or simultaneous occurrence of multiple
species. The new dataset, that we call ZooLake2.0, contains 29499 images,
separated in 35 classes.777Throughout this paper, we use the term class to
indicate a classifier category, and not the taxonomic rank. We show the
abundance distribution in Fig. 2a (blue column). Most of the classes indicate
zooplankton taxa, identified at a genus or species level (e.g. daphnia or
keratella_quadrata). Others, e.g. asterionella, identify phytoplankton
colonies, non-plankton objects such (e.g. dirt or fish), and unidentifiable
objects (e.g. unknown or unknown_plankton). More details on the ZooLake2.0
dataset can be found in App. B. We provide the ZooLake2.0 dataset for free
access at https://doi.org/10.25678/000C6M.

####  II.2.2 OOD test cells

The ZooLake2.0 dataset is treated as a normal dataset for training machine
learning models. We call it the ID dataset, and perform an 80:5:15 splitting
in train, validation and test set, stratified by class abundance. However, if
we want to estimate the OOD performance, we need to go further.

We build 10 OOD test cells, that we will use to validate the OOD performance
of our models. Following previous work on the topic [14, 29], we make sure
that each test cell represents what would actually be found in a generic
deployment context. Each test cell represents a randomly chosen day, with some
mild restrictions. Cells 1-5 are random dates without blooms. For cells 6-10
we did not explicitly impose an absence of blooms, but we restricted the
sampling to days with Chl-a levels under 40 to not have a dataset with a
phytoplankton bloom. The reason for avoiding blooms is that we wanted our test
cells not to be only representing few taxa. Since the plankton populations
vary throughout the day, for each day we evenly sampled from uniformly
distributed hours of the day: 1:00, 6:00, 11:00, 16:00, 21:00. To spare human
annotation effort, when a test cell contains more than 1000 images (9 out of
10 cases), we restrict to 1000 random images. This ensures that human-induced
bias in the image distributions is minimized, and that the test cells
represent what would actually be found in a generic deployment context.

In this work we want to investigate DS in its purest version, i.e. a DS that
is not clearly attributable to any external factor (as could be changing
camera, site, aging of the instrument, and so on). Any of these additional
factors would add up to what we are finding, and be a source of increased DS.
Therefore, all the images (both from ZooLake2.0 and from the OOD test cells)
come from the same camera, deployed at the same depth in the same station,
throughout the whole year. In order to ensure that the aging of the instrument
is not a dominant factor, the ZooLake2.0 images come with a 2-year gap, and
the test cells are randomly sampled in this 2-year gap (Fig. 2-right). With
this construction, we can only impute dataset shift to intrinsic factors of
plankton imaging. In Fig. 3 we show the class distribution, $P(y)$, in each of
the OOD test cells, and of the ID dataset.

Figure 3:  Distribution of classes in the 10 test cells. The column on the
right indicates the day of sampling that the test cell comes from. The ID
dataset is not assigned a date, because it comes from different days (Fig.
2b).

###  II.3 Model training

We train two different families of models, CNNs and Vision Transformers,
including the following deep learning architectures:

Convolutional Networks

MobileNetV3-Large [45], DenseNet-161 [46], EfficientNet-B2 [47],
EfficientNet-B7 [47].

Transformers

DeiT-Base [48], ViT-Base [49], BEiT-Base [50], Swin-Base [51].

The specific choice of the architectures is motivated by previous literature
on plankton classification: DenseNets were found to be well-performing in Ref.
[18], MobileNets and EfficientNets in Ref. [1], DeiTs in Ref. [20], Swin
transformers in Ref. [21], and BEiTs in Ref. [22].

In particular, MobileNets are lightweight models with a performance that is
only slightly lower than that of much bigger models, so they can be used to
systematically test the effectiveness of several techniques, with a relatively
low investment in terms of computing resources.

For all models, we make use of transfer learning [52], using pretrained weight
configurations as an initialization for our models. Each architecture is
trained three times.888The resulting trained network is different because of
the initialization of the final layer, of the randomness in the learning
dynamics induced by stochastic gradient descent (also called gradient noise),
and of data augmentation. This allows us to use the best model based on the
validation dataset, and to create ensemble models. Details on training and
hyperparameter tuning are provided in App. C.

####  II.3.1 Data augmentation

Data augmentation is a common strategy to improve the generalization of
machine learning models. However, the number of possible augmentations is very
large and hard to select. We compare four different kinds of augmentation:

None

This is a baseline to show the model performances when augmentation is not
used.

Basic

These models are trained with the basic augmentations presented in Ref. [1],
which are random rotation and random flipping. Some examples of the augmented
image are shown in Fig. 4. The probability of applying horizontal and vertical
flipping are both 0.5, the degree of rotation is between $0^{\circ}$ and
$180^{\circ}$.

Targeted

This is a set of augmentations which target OOD sensitivity. We calculate the
OOD sensitivities of a model, identify the features to which the model is most
sensitive, and apply augmentations that target those specific features. We
describe it in App. E.

Extra

In addition to applying the basic and targeted augmentations, we add four more
augmentations, which are random posterizing, random solarizing, random
sharpness and random contrast. These are commonly used augmentations that
transform the color characteristics of image.

Figure 4:  Graphical examples of basic augmentations compared to original
images, including random flipping and rotation.

####  II.3.2 Ensemble models

We also train ensemble models, that is, we use $n_{m}=3$ different models to
produce predictions on each image, and determine the class of the image based
on the average confidence vector.

Given $N_{c}$ classes and $n_{m}$ models to be ensembled, we write the
$i^{\mathrm{th}}$ component related to the $j^{\mathrm{th}}$ as $c_{i}^{(j)}$.
The average confidence vector is obtained in two ways:

Arithmetic

For each component $i$, the ensemble prediction confidence is
$c_{i}^{\mathrm{(arith)}}=\frac{1}{n_{m}}\sum_{j}^{n_{m}}c_{i}^{(j)}$.

Geometric

For each component $i$, the ensemble prediction confidence is
$c_{i}^{\mathrm{(geo)}}=\sqrt[n_{m}]{\prod_{j}^{n_{m}}c_{i}^{(j)}}$.

###  II.4 Metrics

We here list the metrics that we use throughout the paper.

####  II.4.1 Per-class descriptors

Let us call the $TP(y),FP(y)$ and $FN(y)$ the true positives, false positives
and false negatives related to class $y$. We can then define the per-class
recall999The recall is sometimes called true-positive rate, or sensitivity
related to class $y$ as

$r(y)=\frac{TP(y)}{TP(y)+FN(y)}\,,$ (7)

the per-class precision as

$p(y)=\frac{TP(y)}{TP(y)+FP(y)}\,,$ (8)

and the per-class $F_{1}$ score as

$F_{1}(y)=2\frac{r(y)\,p(y)}{r(y)+p(y)}\,.$ (9)

Finally, we also define the fall-out (also known as false-positive rate),
which will be useful for adjusted quantification

$f(y)=\frac{FP(y)}{TN(y)+FP(y)}\,.$ (10)

####  II.4.2 Micro-averaged descriptors

Micro-averaged descriptors give the same weight to all the examples in the
testing set, so the errors will be dominated by the largest population.

The accuracy is the number of correct examples out of the total,

$A=\frac{TP}{n}=\frac{\sum_{y}TP(y)}{n}$ (11)

where we defined the total number of true positives as $TP=\sum_{y}TP(y)$, and
the total number of examples $n$, which is connected to the number of examples
per class, $n(y)$, through $n=\sum_{y}n(y)$.

By noticing that $n(y)=TP(y)+FN(y)$, we can also write the accuracy in terms
of the per-class recall,

$A=\sum_{y}\frac{n(y)}{n}\,r(y)\,.$ (12)

From Eq. (12) we can see that the accuracy and the micro-averaged recall are
the same quantity.

####  II.4.3 Expected accuracy

The ratio $\frac{n(y)}{n}$ is trivially an estimator of $P(y)$. Equivalently,
$r(y)$ is related to $P(\mathbf{x}|y)$, since it is a measure of how good the
classifier recognizes elements of a class, given its appearance. If $r(y)$ is
fully representative of an immutable $P(\mathbf{x}|y)$, then future
measurements (also performed on data with the same $P(\mathbf{x}|y)$) will
give the same $r(y)$, regardless of $P(y)$ but not of $P(\mathbf{x}|y)$.

Therefore, in the total absence of dataset shift going from ID to OOD
($P_{\mathrm{ID}}(y)=P_{\mathrm{OOD}}(y)$ and
$P_{\mathrm{ID}}(\mathbf{x}|y)=P_{\mathrm{OOD}}(\mathbf{x}|y)$ do not change),
if we measure the accuracy in a test set, it will remain stable when using
another test set. If the only source of dataset shift is $P(y)$, then $r(y)$
stays stable, and one expects that the accuracy, measured through Eq. (11)
will be equal to

$A_{\mathrm{exp}}=\sum_{y}\frac{n_{\mathrm{OOD}}(y)}{n_{\mathrm{OOD}}}\,r(y)\,.$
(13)

The condition $A=A_{\mathrm{exp}}$ must be met if the dataset shift is purely
distributional.

As for purely compositional DS, it is easier to check that this never happens.
In order to have purely compositional DS, the ID $P(y)$ should stay stable
throughout the OOD cells. Fig. 3 shows that this is not the case.

####  II.4.4 Macro-averaged descriptors

Since $\frac{n(y)}{n}$ in Eq. (12) is an estimator of $P(y)$, this term makes
$A$ dependent of $P(y)$, so it is more influenced by more abundant classes. To
give every class the same influence on $A$, we can replace it with
$\frac{1}{N_{c}}$, that is the same for all classes, and get the macro-
averaged recall, which is not directly affected by the imbalance,

$R=\frac{1}{N_{c}}\sum_{y=1}^{N_{c}}r(y)\,.$ (14)

Equivalently, the macro-averaged precision and $F_{1}$ score are

$\displaystyle P$ $\displaystyle=\frac{1}{N_{c}}\sum_{y=1}^{N_{c}}p(y)\,,$
(15) $\displaystyle F_{1}$
$\displaystyle=\frac{1}{N_{c}}\sum_{y=1}^{N_{c}}F_{1}(y)\,.$ (16)

####  II.4.5 Further macro-averaged metrics for quantification

Metrics for quantification do not check the images one by one, but rather
check how well the abundances of each class are reconstructed. Since in
plankton monitoring one is usually only interested in taxon abundances, it is
normal to resort to such quantities [28, 29].

If we call $\hat{n}(y)$ the abundance estimated by the classifier and $n(y)$
the true abundance, we can define the per-class bias on a dataset as:

$B(y)=\hat{n}(y)-n(y)\,.$ (17)

This quantity only makes sense when defined per-class, since it averages to
zero when taking all the classes ($\sum_{y}[\hat{n}(y)-n(y)]=n-n=0$), but has
the advantage of indicating whether some classes are systematically over- or
under-represented by the classifiers.

The two metrics that we use here are the Normalized Mean Absolute Error (NMAE)

$NMAE=\frac{1}{N_{c}}\sum_{y=1}^{N_{c}}\frac{|\hat{n}(y)-n(y)|}{n(y)}\,,$ (18)

and the Bray-Curtis dissimilarity

$\mathcal{D}=\frac{\sum_{y}|\hat{n}(y)-n(y)|}{\sum_{y}[\hat{n}(y)+n(y)]}\,,$
(19)

which is 0 for two identical distributions, and 1 for two maximally different
ones.

####  II.4.6 Drop in performance

Since we are interested in assessing whether and how much the performance
drops from ID to OOD, we also define some measurements of the drop in
performance. Given a metric $\mathcal{M}$, with a
$\mathcal{M}_{\mathrm{train}}$ value calculated on the training set and a
value $\mathcal{M}_{\mathrm{test}}$ calculated on a generic test set, we
define drop as

$\mathcal{M}^{\mathrm{(drop)}}=\mathcal{M}_{\mathrm{train}}-\mathcal{M}_{\mathrm{test}}\,.$
(20)

By expressing $\mathcal{M}^{\mathrm{(drop)}}$ with respect to the training set
instead of the whole ID dataset, we can compare the DS of the OOD cells with
that of the ID test set, which would serve as a baseline.

####  II.4.7 OOD performances

We have 10 OOD test cells. OOD metrics can be calculated in two ways:

macro-OOD

We calculate the performance separately on each of the OOD sets, and then
average.

micro-OOD

We merge all the OOD sets together, and then calculate the metric on the
resulting aggregated OOD set.

While it is usually desirable to treat the OOD sets one by one (i.e. macro-
OOD), metrics such as the F1 score give artificially low results when
calculated on test sets where some of the classes have a very low number of
examples. Therefore, in such cases micro-OOD performances are more
informative.

###  II.5 Describing the images

We describe the images through 67 standard descriptors, listed in App. D.
Additionally, we define two quantities that we use for a preliminary
description of our data:

##### Crop Index

We want to quantify when and how much plankton images are cropped, due to the
organism being on the border of the camera. To this aim, we define a crop
index, which counts the number of pixels of the organism that overlap with the
image boundary. Since our images have a black background, this quantity is
simply calculated by counting the number of non-black pixels at the image
boundary (pixels with value smaller than 5 are still considered as black).

##### Blurriness

In determining the image’s degree of blurriness, we initially compute its
Laplacian, a two-dimensional isotropic representation of the image’s second
spatial derivative. The Laplacian highlights regions in the image where there
are rapid intensity change. The degree of blurriness of an image is then
calculated by the average of the absolute values of Laplacian.

###  II.6 Quantifying dataset shift

####  II.6.1 Quantifying distributional dataset shift

To calculate the dissimilarity between two distributions of classes, $P(y)$
and $Q(y)$, we use the scalar product,

$d_{y}=1-P\cdot Q\equiv 1-\frac{1}{N_{c}}\sum_{y=1}^{N_{c}}P(y)Q(y)\,.$ (21)

The quantity $d_{y}$ is 0 if $P(y)=Q(y)$, and 1 if the classes appearing in
$P(y)$ do not appear in $Q(y)$ and vice versa.

####  II.6.2 Estimating compositional dataset shift

To estimate compositional dataset shift, we need to have a measure of the
dissimilarity of $P(\mathbf{x}|y)$ between training and deployment.

The quantity $P(\mathbf{x}|y)$ is conditioned on the class, so the total
distance must be a sum over the classes of a per-class distance. The input
$\mathbf{x}$ is the image, and we can describe it in two ways:

$\mathbf{x}$ is a vector of pixels

This is a natural choice, since it is exactly the input that the models take.
However, if we then compute the dissimilarity on a per-pixel basis, we are
neglecting correlations between pixels which can be important in the
definition of the features of a family of images. One way of taking into
account some kind of correlations with pixels is to perform a PCA. However,
the components of this vector are not interpretable, and as we will see this
is a desirable property.

$\mathbf{x}$ is a vector of features

An alternative solution is to extract a series of computer-vision descriptors
of the image (App. D), and then use those as $\mathbf{x}$. The shortcoming of
this solution is that it is not clear which and how many different descriptors
are necessary, that they are potentially collinear, and they act on different
scales.
To solve this issue, we resort to a relatively large number of descriptors,
and apply a principal component analysis (PCA), in order to obtain an
orthonormal space on which to calculate dissimilarities.
When comparing the importance of different features, since we are not summing
them with each other, we do not need to orthogonalize. Therefore, no PCA is
needed, but we still need to standardize the features so that their scale does
not artificially make some distances larger than others.

##### Hellinger Distance

In order to estimate shifts in $P(\mathbf{x}|y)$, we define and compare
several dissimilarity measures. In the main paper, we only report the
definition of the Hellinger distance [53], which we found to best correlate
with the performance drops. In App. G we also show results using the
Wasserstein distance and the Kullback-Leibler divergence.

The Hellinger distance was already used in plankton classification [29, 41].
Given two normalized distributions $p$ and $q$ defined on a discrete one-
dimensional support (with bins $i=1,\ldots,n_{\mathrm{bins}}$), the Hellinger
distance is

$D_{\mathcal{H}}^{f}=\frac{1}{\sqrt{2}}\sqrt{\sum_{i=1}^{n_{\mathrm{bins}}}\left(\sqrt{p_{i}}-\sqrt{q_{i}}\right)^{2}}$
(22)

The superscript $f$ in Eq. (22) represents a single dimension of $\mathbf{x}$
(so a pixel, a principal component, or a feature of the image). Since the
values of dimension $f$ are min-max normalized, the distance ranges from 0 to
1. The full distance between two distributions is obtained by averaging the
Hellinger distance along all the components,

$D_{\mathcal{H}}=\frac{1}{N_{f}}\sum_{f=1}^{N_{f}}D_{\mathcal{H}}^{f}\,.$ (23)

Note that $D_{\mathcal{H}}$ especially in the presence of small abundances per
class, is dependent on the binning. For this reason, unless specified
otherwise, we restrict $D_{\mathcal{H}}$ to those classes that have at least
10 images in the OOD cells.

###  II.7 Using the distances

By checking how much the OOD performance drops, we can use the dissimilarities
defined in Sec. II.6 to address several questions. This procedure expands on
ideas initially presented in Ref. [29], which showed that cross-validation
folds with larger $D_{\mathcal{H}}$ have a lower performance. As described in
the following, each question needs the distances to be calculated in a
slightly different way.

##### Are the model performances directly related to dataset shift? (Sec.
III.2.3)

Here, we want to estimate how much the $P(\mathbf{x}|y)$ changes, and then
average over $y$. We therefore need to calculate the distances on a per-class
basis, and then take the average over the classes. This corresponds to taking
a macro-averaged distance.

Since $\mathbf{x}$ represents several features that could be mutually
correlated, and we want an orthonormal basis, we perform a PCA, and select the
number of principal components such that the amount of explained variance is
greater than 95%. There are large differences between the ranges of initial
features, those features with larger ranges will dominate over those with
smaller range and lead to biased result. Therefore, we standardize the
features by subtracting the mean and dividing by the standard deviation, so
that each feature contributes similarly to PCA. We then calculate the distance
according to each principal component following Eqs. (22) and (23).

##### Which are the image features which most influence the DS robustness?
(Sec. III.2.3)

When comparing different features, we do not need to perform the PCA, but it
is important that we rescale different features in such a way that they are on
the same scale (otherwise, the distances related to some features may
artificially seem larger than others). Since the value of different features
are in different scales, and because of the presence of outliers, we apply
some transformations on the original feature data before distance calculation.
First, the feature values are standardized by subtracting the mean and
dividing by the standard deviation to ensure they are all on a same scale.
Then a logarithmic transformation is applied to reduce the impact of outliers,

$x^{\prime}=\begin{cases}\log(1+x),&\text{if $x\geq 0$}\\\
-\log(1-x),&\text{otherwise}.\end{cases}$ (24)

Since with this construction the components are not orthogonalized, collinear
features will have similar values of $D_{\mathcal{H}}$, but this is not a
problem for the analyses presented in this paper (for example, it is not a
problem if both height and area appear as important).

##### Which classes are harder/easier more robust under dataset shift? (Sec.
III.2.3)

To estimate which classes suffer most from dataset shift we take the distance
on a per-class basis. Also here we perform a PCA on the features.

##### How well will a classifier perform on a new unlabeled dataset? (Sec.
III.2.3)

If presented with a totally new, unlabeled dataset, we do not have knowledge
on the classes which compose it, so we cannot estimate $P(\mathbf{x}|y)$.
Therefore we directly estimate the distances on the distribution
$P(\mathbf{x})$, which corresponds to calculating a micro-averaged distance.
Here we perform a PCA on the input features as well.

###  II.8 Performance Sensitivity

The bigger the performance drop as a function of the dissimilarity, the bigger
we expect that a specific feature or class will be critical in the OOD
robustness of the classification. We can therefore think in terms of a
quantity

$\displaystyle Q_{\mathrm{acc}}$
$\displaystyle=\frac{A_{\mathrm{train}}-A_{\mathrm{test}}}{D(\text{train,
test})}\,,$ (25) $\displaystyle Q_{\mathrm{F_{1}}}$
$\displaystyle=\frac{{F_{1}}_{\mathrm{train}}-{F_{1}}_{\mathrm{test}}}{D(\text{train,
test})}\,,$ (26)

where $D(\text{train, test})$ is any of the distances we defined in Sec. II.6.
The larger $Q$, the larger the influence of dataset shift. In the limit of
small variations, $\frac{dA}{d\mathbf{x}}$ (or $\frac{dF_{1}}{d\mathbf{x}}$),
this quantity is usually called a sensitivity in statistics [54] or a
susceptibility in physics [55].

###  II.9 Improving performances under dataset shift

We test simple methods for countering dataset shift, based on model choice,
model training, and interpretation of the results.

####  II.9.1 RGB reweighting

If we assume that color variations could have a strong effect in the OOD
performances, it can make sense to make sure to impose that the intensity
distributions along the overall color distributions do not vary when the input
data changes. Therefore, we rescale the input data in such a way that it
matches the RGB distributions of the training set (which is different from
that of ImageNet, given that our images have a black background). Operatively,
the mean and standard deviation of RGB values of the whole training set are
calculated. Then during training, validation and testing, the RGB channels of
all images are standardized by the means and standard deviations.

####  II.9.2 Targeted augmentations

Data augmentation can help increase the robustness of models against
distribution shift [56]. Therefore, we test whether we can address OOD
performance drops with data augmentations that are generally not used within
plankton classification.

To identify which kinds of augmentations have the best potential of
improvement, we measure which features have the largest sensitivities $Q$,
defined in Sec. II.8, and perform augmentations based on these sensitivities.

####  II.9.3 Ensembling

Besides choosing the models with the best robustness against dataset shift, we
assess whether and how using ensemble models produces performances that are
more robust under dataset shift. We compare the model performances of single
model, arithmetic averaged ensemble model and geometric averaged ensemble
model, as mentioned in Sec. II.3.2.

####  II.9.4 Architecture

Different architectures may have different robustness to DS. In particular, it
was argued that vision transformers are more robust than CNNs [57]. We compare
the performance of 8 selected architectures listed in Sec. II.3.

####  II.9.5 Test-time augmentation

Test-time augmentation (TTA) consists of performing test predictions on
augmented versions of an image, in addition to the original one, and taking
the average prediction. For example, if one takes an image, and the same image
rotated of $180^{\circ}$, the model prediction should not change. Averaging
the predictions over these two images will produced a more stable confidence
vector where the noise is flattened out.

To our knowledge, TTA was proposed for the first time in Ref. [58], where they
ensembled the prediction vector of five cropped patches (four corner patches
and one center patch) while evaluating the model performance on test set. This
gives the model a more robust prediction. However, most of the images in our
plankton dataset have only one organism in the center surrounded by a dark
background, so cropping does not introduce much diversity. Instead, we apply
rotation during test time. The final prediction of an image is the average of
the predictions on $k$ rotated versions of image, with equally spaced angles.

####  II.9.6 Quantification through adjusted counts

Another solution is to keep the classifiers untouched, and reinterpret the
population counts. This has the advantage of not requiring retraining models.

Here, we use unsupervised quantification algorithms proposed in earlier
plankton literature [29, 41, 42]. The basic quantification algorithm is
Classify and Count (CC), which estimate an abundance by summing up the number
of counts that the classifier associates to that class:

$n_{\mathrm{CC}}(y)=TP(y)+FP(y)\,.$ (27)

Adjusted counts methods correct the CC classification by using knowledge of
the ID test performance of the models. Although these methods are generally
only thought for distributional dataset shift (since the knowledge of the ID
performance could become obsolete), we check whether within lake plankton
classification they can also be helpful in a more generic setting.

The adjusted count method [59], adjusts the CC counts based on the expectation
of how often false positives are found:

$n_{\mathrm{AC}}(y)=n_{\mathrm{test}}\frac{n_{\mathrm{CC}}(y)/n_{\mathrm{test}}-f(y)}{r(y)-f(y)}\,.$
(28)

Another method consists in using the classifier softmax101010The we specify
”softmax” because the confidence vectors should be normalized
($\sum_{y}c(y)=1$), which is a property of softmax. confidences $c_{i}(y)$ as
probabilities for an adjusted count [60]. Following Ref. [41] and call it
probabilistic (adjusted) count. We have both a probabilistic classify and
count (PCC) and a probabilistic adjusted count (PAC), respectively defined as

$\displaystyle n_{\mathrm{PCC}}(y)=$
$\displaystyle\sum_{i=1}^{n_{\mathrm{test}}}c_{i}(y)\,,$ (29) $\displaystyle
n_{\mathrm{PAC}}(y)=$ $\displaystyle
n_{\mathrm{test}}\frac{n_{\mathrm{PCC}}(y)/n_{\mathrm{test}}-f^{\mathrm{pa}}(y)}{r^{\mathrm{pa}}(y)-f^{\mathrm{pa}}(y)}\,.$
(30)

where $r^{\mathrm{pa}}(y)$ and $f^{\mathrm{pa}}(y)$ are the probability-
average recall and fall-out, defined through measurements on the training set,

$\displaystyle r^{\mathrm{pa}}(y)=$
$\displaystyle\frac{\displaystyle\sum_{i\in
y}^{n_{\mathrm{train}}(y)}c_{i}(y)}{n_{\mathrm{train}}(y)}\,,$ (31)
$\displaystyle f^{\mathrm{pa}}(y)=$
$\displaystyle\frac{\displaystyle\sum_{i\notin
y}^{[n_{\mathrm{train}}-n_{\mathrm{train}}(y)]}c_{i}(y)}{\bigg{[}n_{\mathrm{train}}-n_{\mathrm{train}}(y)\bigg{]}}\,,$
(32)

where $n_{\mathrm{train}}(y)$ is the number of training images belonging to
class $y$, and $n_{\mathrm{train}}$ is the total number of training images.

In Ref. [41], they see that quantification algorithms do not add much with
respect to Classify and Count, unless the models have low performances. So we
can think that they are normally not useful, except when the model starts
failing.

We highlight that one can also use semi-supervised quantification methods.
These involve periodic checks by a human expert manually verifying by how much
the counts should be adjusted. These yielded good results within binary
plankton classification [42], but are out of the scope of this work.

####  II.9.7 Abstention

One can hypothesize that the classifier confidence correlates with the
probability of a prediction being correct. In that case, one can decide to
discard the low-confidence predictions, thus increasing the model precision in
exchange for a decrease in recall.

Formally, we can decide to define an abstention threshold $\theta\in[0,1]$,
and only accept predictions where the confidence on the predicted class is
larger than $\theta$. With $\theta=0$ there is no abstention and all images
get classified. With $\theta=1$ all the images get discarded.

##  III Assessment

###  III.1 Identifying Performance Degradation

####  III.1.1 OOD performance degrades

Figure 5:  Four different performance metrics (accuracy, F1 score, Bray-Curtis
dissimilarity, Normalized Mean Absolute Error) in the ID training and test set
(blue), in the aggregated OOD data (red), and in the individual OOD test
cells, obtained with a MobileNet (see Sec. C).

We trained models following the procedures we described in Ref. [1], and
compare their ID test performance with the OOD performance. In Fig. 5 we
report four performance metrics for a MobileNet,111111The same qualitative
conclusions hold for all the architectures we trained, as we will show further
in the appendix (Fig. 37). two of which are based on single-image performance
(accuracy and F1-score), and two of which are based on estimations of the
abundance of each class (see Sec. II.4). Both if we use the individual test
cells (macro-OOD, orange points), or if we merge all the OOD cells into a
single one (micro-OOD, red), the ID test performance is drastically higher
than the actual performance when faced with real-life situations. We stress
that the performance strongly varies among the OOD cells. This gives us an
indication of the degree of fluctuations in the performances which should be
expected: in different days of deployment, the classification performance can
be quite different. Furthermore, there can be specific days where the
performance is particularly different from usual. One example can be the OOD10
cell, where the accuracy is high, but the F1 score is low. This is because
OOD10 comes from a day in which plankton community was dominated by filament.
Since the model classifies filament well, the accuracy is high. The F1 score
is however low, because the other taxonomic units have very few counts, and
the F1 score is typically exaggeratedly low when the number of true positives
is close to zero.

In Fig. 6, we show that different classes suffer differently from OOD: while
nauplius is practically unaffected, the classification of eudiaptomus
deteriorates visibly on OOD data. The unknown class has the greatest OOD drop.
This is expected, since it is a container class which gathers all sorts of
disparate objects.121212We highlight that the unknown class does not represent
instances where the classifier abstained. They instead represent instances
where the taxonomists were not able to assign a ground truth. In App. H we
show the performance drop for all the classes.

Figure 6: (a) Recall and (b) F1 score in the training set, ID test set, and
OOD test sets, obtained by a MobileNet (see Sec. C) for three classes:
nauplius, eudiaptomus and unknown. More classes are shown in App. H.

In Fig. 7 we show a scatter plot of the predicted vs true abundance in the ID
test and in the overall OOD dataset. The ID points are concentrated around the
1:1 line, and the OOD predictions are further away. While the overall OOD
performance is still by many standards acceptable, the abundances of some
classes are wrongly estimated.131313 The reader might be interested whether
this problem can be overcome by aggregating similar classes together. In App.
I we show that, though the overall performances go up, a performance drop also
persists if we group the classes describing coarser taxonomic units.

Figure 7: Scatter plot of Each point represents a single class, either in the
ID test set (grey) or in the aggregated OOD set (pink). A perfect
classification would fall on the 1:1 line (red line). The inset shows the same
data as the main set, but on a linear scale. In App. J, we show this same
scatter plot for the individual OOD cells.

###  III.2 Diagnosing performance degradation

####  III.2.1 Why is performance degrading? A first analysis

Figure 8:  (a): Distribution of the crop index in the ID and in the OOD. The
gray part is the overlap of two distributions. (b): Distribution of the
blurriness in the ID and in the OOD. Crop index and blurriness are defined in
Sec. II.5.

As a start, we can try to explain this phenomenon with the most intuitive
explanation possible: there could be more hard-to-classify images in the OOD
cells than in the ID, due, for example, to sampling bias. The two most
intuitive reasons that we can attribute to image hardness, are that the image
was cropped (only a part of the organism appears in the picture), or it was
blurry (the organism is out of focus).

In Fig. 8, we compare the ID and OOD data according to these two quantities.
In Fig. 8a we see that the crop index (defined in Sec. II.5) is equally
prevalent in the ID and OOD set, so cropping does not seem to be a driver of
the OOD degradation. Similar conclusions are valid for the blurriness (Fig.
8b), which is comparable in the two sets. Therefore, in order to understand
the OOD drop in performance, we need to resort to a more sophisticated
analysis.

Figure 9:  (a): Plot of F1 drop vs $d_{y}$. (b): Correlation between $d_{y}$
and $D_{\mathcal{H}}$. (c): Plot of expected accuracy vs $d_{y}$.

####  III.2.2 Performance in terms of P(y)

We now explore the dependence of the performance drop on distributional DS,
i.e. variations in the distribution of classes, $P(y)$, from ID to OOD. From
Fig. 3, we already know that distributional DS is present. We quantify it
through the distance $d_{y}$ [Eq. (21)]. As we show in Fig. 9a, the
performance drop has some correlation with the amount of shift. However, we
show in Fig. 9b that distributional DS correlates with compositional DS, so
the two sources of DS cannot be fully untangled. This is reasonable, since for
example, the occurring of a bloom will change the relative abundances of taxa,
and bring more diversity in plankton morphology.

We can however discard that the main driver of the performance drop is
distributional shift. If in fact the performance degradation is driven by
dataset shift alone, the OOD expected accuracy [Eq. (13)] and accuracy should
match. Instead, the measured accuracy is systematically lower than the
expected accuracy (Fig. 9c), and the size of this drop is similar to the
overall drop from ID to OOD.

####  III.2.3 Performance in terms of $P(\mathbf{x}|y)$

##### Overall performance

Since $P(y)$ cannot hold the whole responsibility of the degradation, we need
to inspect $P(\mathbf{x}|y)$. In particular, we want to estimate how much the
performance changes when varying $P(\mathbf{x}|y)$ from the training set to
each of the (ID or OOD) test sets. We estimate the change in $P(\mathbf{x}|y)$
from training to tests sets through the distances defined in Sec. II.6 and
App. G. In the main text, we show results for the Hellinger distance
$D_{\mathcal{H}}$, and show analogous results for other kinds of distances in
App. G.

In Fig. 10 (squares), we show the F1 and accuracy drop in each test cell, as a
function $D_{\mathcal{H}}$. We find a visible correlation between performance
drop and distance from the training set, with the F1 score correlating better
than the accuracy (this is also found with other metrics, App. G, as well as
in Sec. III.2.3).

Figure 10: (a): F1 drop [Eq. (20)] as a function of the Hellinger distance
from the training set [Eq. (23)]. Each dot indicates a specific class in an
OOD cell. The dashed lines are the percentile 5, 50 and 95 regressions. The
solid line is the mean regression. Squares indicate a single OOD cell (error
bars are the standard error among the different classes in each test cell).
(b): Same, for the recall. Similar plots for other distances and architectures
are shown in App. G and App. K.2.

The dots represent instead each single class of each dataset. We show the mean
regression, and the quantiles 5, 50 and 95 of the distribution. All increase
steadily with $D_{\mathcal{H}}$.141414In particular, if $D_{\mathcal{H}}$ is
not a good indicator of the performance drop (as for example the KL
divergence, App. G), one could have that the highest quantiles even decrease
with $D_{\mathcal{H}}$, i.e. a fraction of the classes become easier, but this
is not the case. Note that regressions lines on the OOD cells (on the squares
instead of the dots) would give an even higher slope and better correlation.

##### Distinguishing features

One major question we wish to address about DS degradation is: what is
changing in the images, that makes the performance go down?. We already saw in
Sec. III.2.1 that it is not the cropping nor the blurriness. We cannot connect
the answer of this question to features such as water temperature or organism
health, which cannot be easily inferred from the images. Therefore, we address
the explanation of the OOD performance in terms of computer vision features,
which can be extracted from the images (App. D). To do so, we compute the
Hellinger distance from the ID training set for each single feature, using Eq.
(22).151515These features are not used for model training, but rather only to
quantify the DS.

We perform a simple linear regression of the F1 drop as a function of
$D_{\mathcal{H}}^{f}$. This tells us, to first order, how the performance
varies while increasing the distance related to $f$. The slope of these
regressions is the sensitivity $Q_{F_{1}}$ defined in Eq. (26). We plot those
regression curves in Fig. 11a. For clarity reasons, we only plot the curves
with $Q_{F_{1}}>0.7$, with an additional line representing the average of the
remaining features. We summarize the values of $Q_{F_{1}}$ in Fig. 11b, where
we see that no specific feature stands out. The most relevant features are
(definitions are listed in App. D) saturation_std, intensity_G_std,
eccentricity, compactness, formfactor, convexity, hull_area, height, ESD,
blurriness, angle_rot, image moments $M_{ij}$, $\mu_{ij}$, $\eta_{ij}$ and
$I_{i}$. These are related to color, shape, object size, blurriness and
orientation, indicating that no single computer vision feature (or family of
features) is dominating DS.

Figure 11:  (a): Performance drop as a function of $D_{\mathcal{H}}^{f}$. For
every feature $f$, we calculate regression lines, and show those related to
the features with slopes $Q_{F_{1}}>0.7$. (b): Measured value of $Q_{F_{1}}$,
for each feature. All results in this figure are for a MobileNet model,
trained with basic augmentations.

##### Distinguishing classes

Similarly as done for the single features, we can also check which classes are
more sensitive to DS. However, we can only do this for those classes which
occur at least 10 times in at least 4 OOD cells, to ensure reasonable
correlations. We show them in Fig. 12a. We summarize the per-class
sensitivities in Fig. 12b, where we can see that unknown, dinobryon and
keratella_cochlearis are the classes that is most sensitive to DS. Note,
however, that this result is model-dependent. In fact, we will show that we
can produce a model where the sensitivity of dinobryon and
keratella_cochlearis are strongly reduced (App. K.2).

Figure 12: (a): Per-class F1 drop as a function of the Hellinger distance.
(b): Sensitivity $Q_{F_{1}}$ related to each class of those for which enough
data is available. In App. K.2 we show the same figure for the BEsT model.

##### Unlabeled datasets

For new data, we can give a rough estimate of the performance degradation. We
do this by estimating the distance in the $P(\mathbf{x})$ instead of the
$P(\mathbf{x}|y)$, as described in Sec. III.2.3. We find that this procedure
seems to work reasonably well for the F1 score drop, but not for the accuracy.
Details are in App. K.1.

###  III.3 Cure

With the objective of improving the OOD performance, we test seven different
methods:

  1. 1.

RGB channel standardization (Sec. III.3.1)

  2. 2.

Targeted augmentations (Sec. III.3.2)

  3. 3.

Architecture selection (Sec. III.3.3)

  4. 4.

Ensemble learning (Sec. III.3.4)

  5. 5.

Test-time augmentation (Sec. III.3.5)

  6. 6.

Adjusted counts (Sec. III.3.6)

  7. 7.

Abstention (Sec. III.3.7)

In the following, we describe the effectiveness of each of those methods for
our application case.

####  III.3.1 RGB channel standardization

Since, as noted in Sec. III.2.3, saturation_std appears as one of the most
sensitive features to DS. Therefore, we standardize the RGB color
distributions, in such a way that those of the OOD cells match those of the
training set, as described in Sec. II.9.1. However, as we show in App. F, this
procedure does not improve the performances.

####  III.3.2 Targeted augmentations

We compare the four strategies of image augmentation in Fig. 13.

Figure 13: Comparison between the different augmentation schemes (Sec. II.3.1,
App, E), for a MobileNet model. Note that, although the ID performances are
similar for any of the three applied augmentations, their OOD performances
spread.

The targeted augmentation (the choosing of targeted augmentations is discussed
in detail in App. E) emerges as a clear winner (the green points have the best
performance across all OOD cells). Note that, while the difference between
different augmentation schemes is similar in the ID test set, it becomes
larger in the OOD cells. We also highlight that training a model with targeted
takes only slightly longer than with basic augmentation, but much less than
the extra augmentation (Sec. II.3.1).

####  III.3.3 Architecture selection

We now compare the performance of different architectures: MobileNetV3 [45],
EfficientNet-B2 [47], EfficientNet-B7 [47], DenseNet [46], ViT [49], DeiT
[61], Swin-T [51] and BEiT [50]. The specific choice of the architectures is
motivated by previous literature on plankton classification: DenseNets were
found to be well-performing in Ref. [18], EfficientNets in Ref. [1], DeiTs in
Ref. [20], Swin transformers in Ref. [21], and BEiTs in Ref. [22].

We compare the OOD performance of those architectures in Fig. 14, and show the
performances on the single OOD cells in App. L. The best CNN is
EfficientNet-B7, which has high F1 score and BC, and low Accuracy and NMAE.
The transformers perform generally better (except for the BC), and in
particular the BEiT has the best performances in terms of population metrics
(BC and NMAE). Therefore we select the BEiT architecture.161616We show in App.
M that similar conclusions hold when using ensemble models.

Figure 14:  The classification results of basic models using eight
architectures, tested on micro and macro OOD.

####  III.3.4 Ensemble learning

As shown several times before, using several models to increase the quality of
the predictions can strongly increase the quality of the classification [18,
1, 20, 22].

Here, we focus on ensembling over several replicates of the same architecture:
we train the same model several times, with different initial weight
configurations. This is more convenient than ensembling over different
architectures, since the targeted augmentations are chosen on a per-model
basis, and the overall performance is similar (App. M).

In Fig. 15 we show the result of ensembling on the OOD cells through 3 single
BEiT learners, trained with targeted augmentation.

Figure 15: OOD performance without ensembling, and with two different kinds of
ensembling, as defined in Sec. II.3.2. We report two OOD metrics: micro OOD is
the performance over the aggregated OOD cell. Macro OOD is the average of the
performances across test cells.

The improvement is sizable, regardless of the chosen performance metric, and
geometric averaging seems slightly better. A reason can be that, with
geometric averaging, if one of the learners gives a low score to a specific
class, this will be more relevant than the possibly high scores given by the
other learners.

####  III.3.5 Test-time augmentation

We now test the efficacy of test-time augmentation (TTA), consisting in
producing predictions on augmented versions of each test image, in addition to
the images themselves. With TTA, it is a good practice that the augmented
images are true positives, and not deformed versions of the image that may
never occur in reality. This reduces the number of possible augmentations
available. As we show in Fig. 13, despite the targeted data augmentation,
models still suffer OOD degradation, and apparently one of the strongest
factors is rotation (as shown in Fig. 11b). Therefore, we use rotated images
for the TTA (see Sec. II.9.5).

Figure 16: Performance on aggregated OOD dataset under various types of TTA,
of a 3-BEiT geometric ensemble trained with targeted data augmentation.

From Fig. 16, we see that increasing the number of angles generally increases
the performance.171717In App. N we also try flipping, but we find no advantage
with respect to only doing rotations. We must however consider that the
prediction time increases proportionally to the number of TTAs, so a lower
number of angles is better in terms of wallclock time. Therefore, we choose to
limit the number of angles $k$ to 4, since the gain is small after that.

####  III.3.6 Adjusted counts

We also test the efficacy of methods that adjust the counts of the populations
based on the ID metrics, since these methods seemed to have an effect in
improving the classification performances in previous work [28, 41, 42].
However, as we see from Fig. 17, none of the adjusted counts methods
outperforms the vanilla classify and count procedure. In fact, although it is
never the best, it is the second best both in terms of Bray-Curtis
dissimilarity and of Normalized Mean Absolute Error. Therefore, we do not
adjust counts in our final model.

Figure 17: Comparison of different ways of performing adjusted counts. Methods
of calculating these adjusted count are described in Sec. II.9.6.

Since these methods are thought to balance the effects of distributional (not
compositional) DS, the failure of count adjusting is consistent with our
previous observations that DS is a mainly compositional (not distributional)
problem.

####  III.3.7 Confidences and Abstention

We show in Fig. 18a, that the model confidences are high when ID, and lower
when OOD. We see that the average confidences correlate well with
$D_{\mathcal{H}}$ (Fig. 18b), indicating that both quantities are good
indicators of how out-of-distribution the OOD test cells are.

This can be exploited, as already previously done [34], by utilizing
abstention, which consists in asking the models to abstain from classifying
when the confidence is too low. This increases the precision, and decreases
the recall. Depending on the use case, it might be beneficial, when higher
confidences correlate with better guesses, as it happens in our case (Fig.
18c). Since the performance under abstention depends on the chosen abstention
threshold $\theta$, instead of providing performances at a given threshold, we
show in Table. 1 how precision and recall vary when changing this threshold.
For the performances reported in the rest of the paper, we do not use
abstention (this is equivalent to using a threshold $\theta=0$).

Figure 18:  (a): Prediction confidence of the BEsT model. (b): Average confidence on each OOD cell, as a function of their Hellinger distance from the ID training set. (c): Average accuracy in each OOD cell, as a function of the prediction confidence.  Table 1: Abstention result, macro-averaged precision and recall with applying different confidence threshold during prediction on OOD. Threshold | %classified | Precision | Recall
---|---|---|---
0 | 100 | 0.78 | 0.82
0.1 | 99.72 | 0.78 | 0.81
0.2 | 98.19 | 0.80 | 0.81
0.3 | 95.49 | 0.81 | 0.79
0.4 | 91.20 | 0.83 | 0.78
0.5 | 86.44 | 0.85 | 0.76
0.6 | 81.84 | 0.87 | 0.74
0.7 | 77.06 | 0.87 | 0.72
0.8 | 71.10 | 0.90 | 0.68
0.9 | 63.41 | 0.92 | 0.64
0.95 | 56.95 | 0.93 | 0.61
0.99 | 45.87 | 0.98 | 0.56
0.999 | 33.41 | 0.98 | 0.48

####  III.3.8 Ablation study

We show several methods that can be used together, and show that they improve
the OOD performance. The reader can still be interested in knowing whether
each of these methods is useful, with and without the presence of the others.
Therefore, we perform an ablation study on BEiT models, comparing the
performances that are obtained with each combination of model improvements
(Tab. 2). We can see that though the most crucial one is ensembling, best
performances are obtained by mixing all the approaches. On the right part of
the table, we compare how each methodology slows down the training and the
inference times, since this needs to be taken into account as an additional
factor, on top of sheer performance.

Table 2: BEiT model performance metrics on different ablation settings. Ensemble | TTA | Targeted aug. |  Accuracy $\uparrow$ |  F1-score $\uparrow$ |  BC $\downarrow$ |  NMAE $\downarrow$ | Training cost | Inference cost
---|---|---|---|---|---|---|---|---
(Micro-OOD/Macro-OOD)
✗ | ✗ | ✗ | 0.782/0.784 | 0.673/0.516 | 0.117/0.132 | 0.441/0.506 | $\mathcal{O}(1)$ | $\mathcal{O}(1)$
✓ | ✗ | ✗ | 0.821/0.822 | 0.704/0.569 |  0.085/0.100 |  0.255/0.413 | $\mathcal{O}(n)$ | $\mathcal{O}(n)$
✗ | ✓ | ✗ | 0.795/0.797 | 0.702/0.557 | 0.111/0.127 | 0.363/0.474 | $\mathcal{O}(1)$ | $\mathcal{O}(m)$
✗ | ✗ | ✓ | 0.793/0.796 | 0.716/0.556 | 0.100/0.113 |  0.256/0.464 | $\mathcal{O}(1)$ | $\mathcal{O}(1)$
✓ | ✓ | ✗ | 0.824/0.825 | 0.714/0.588 |  0.085/0.102 |  0.258/0.399 | $\mathcal{O}(n)$ | $\mathcal{O}(nm)$
✓ | ✗ | ✓ | 0.818/0.821 | 0.751/0.581 | 0.088/0.102 | 0.310/0.453 | $\mathcal{O}(n)$ | $\mathcal{O}(n)$
✗ | ✓ | ✓ | 0.812/0.815 | 0.746/0.581 | 0.096/0.109 | 0.269/0.412 | $\mathcal{O}(1)$ | $\mathcal{O}(m)$
✓ | ✓ | ✓ |  0.829/0.832 |  0.761/0.610 |  0.086/0.100 |  0.261/0.398 | $\mathcal{O}(n)$ | $\mathcal{O}(nm)$

###  III.4 The final model

In Fig. 19a we show a scatter plot of the true vs predicted OOD counts, for
each population. We show two models: the baseline model and the final selected
model, which we call BEsT model. The improvement is steady across more or less
all the classes. Fig. 19c displays how well the final model performs on OOD
data.

The only two classes with a large number of misclassified images are rotifers
and unknown, which get mixed one with the other. These are container classes,
which gather in themselves a large number of different kinds of objects. One
way of increasing this performance would be to distinguish different kinds of
subclasses among them. Since, however, one of these classes expresses objects
that were not identifiable by the taxonomists, there is a possibility that the
unknown images classified as rotifers are actually rotifers.

In Fig. 19b, we show the performances of our final model, both ID and OOD. The
reader might argue that, while it is nice that the new performances are better
than the baseline ones, the baseline performances are not that bad either.
However, what we are testing here, is the solidity of our models to unknown
changes in the future. The new models are also more robust to DS, as is
demonstrated by the smaller sensitivities depicted in Fig. 19d and in App.
K.2. In fact, the performances that we are reporting are obtained on the same
the OOD cells that were used to perform the model selection. Therefore, these
are validation and not test performances. True performances should be given on
new OOD data, so models with a lower sensitivity give us higher trust that
they will be robust to varying conditions.

Figure 19:  (a): Scatter plot of true vs predicted abundance, for the baseline
model, and the BEsT model. (b): Train, test, micro-OOD and macro-OOD
performances of the BEsT model. (c): Bar plot showing the true and the
predicted plankton abundances, for the BEsT model, on the aggregated OOD data.
(d): Sensitivities $Q_{F_{1}}$ for the BEsT model.

##  IV Discussion

##### Summary

Deploying classifiers without accounting for DS has the risk of providing
highly deceptive results. Here, we show that DS is indeed a problem with
plankton classification, we display a formalism to describe it, and a
procedure to address it.

In particular, we address classification of lake plankton from the DSPC
installed in lake Greifensee in Switzerland, focusing on dataset shift (DS).
Our main aim is to assess whether DS can be a hindrance, and to find a
classifier which is robust to DS. We see that the dataset shift is present,
that it negatively influences model performances, and that compositional
dataset shift is driving the performance decrease. We then study the
sensitivity of our models to DS. This depends on the specific architectures
and how they are trained. The model performance under DS is mostly affected by
color, rotation and shape features and less, for example, by noise in the
image. Due to the imbalance of class abundances, we cannot safely assess the
sensitivity of all classes. Among those with enough available data, the most
sensitive to DS are dinobryon, keratella_cochlearis and unknown. Finally, we
evaluate a series of methods which can improve the model robustness to DS.
Some of those, such as adjusted counting and RGB channel standardization, do
not work, while others, such as ensembling and test-time augmentation,
increase robustness. In Fig. 20, we summarize the steps that we took to get
our best model, in terms of OOD performance.

Figure 20: OOD performances after each step of model improvement. The reader
interested in these steps applied in a different order can consult Tab. 2.

##### Other possible strategies to address dataset shift

One way to increase OOD performance is to increase ID performance. In fact,
there is evidence that models that have better in-domain performance will also
have a better OOD performance [62, 63]. However, it is important to keep in
mind that the ID performance is not representative of the true performance of
the model, and one should always have an estimate of the OOD performance.

Beyond those that we already tested, there are many other ways that one could
address DS in our plankton classification problem. For example, one could use
Domain Adaptation [64, 65] to infer the final $P(\mathbf{x},y)$ and
consequently reweigh the loss function.181818Software for this already exists,
e.g. https://domainadaptation.org This is doable in our case, because
$P(\mathbf{x},y)=P(y|\mathbf{x})P(\mathbf{x})$ and, since one potentially has
access to a large number of unlabeled images (all the images that need to be
classified), we have access to the $P(\mathbf{x})$.

More related to the results of this specific work, we see that rotations still
play an important role, while translations are marginal, since our images are
usually centered. Rotation-invariant CNNs might therefore be of help [66].
Also, we find that the main source of OOD error in our tests is that the
rotifer and the unknown classes get mixed. Since these classes include many
kinds of different objects, robustness would likely increase if one was able
to create reasonable sub-classes out of them. For example, unknown with some
specific shapes and unknown with some specific colors.

##### How serious is DS?

Although we find that DS is present, and it strongly affects accuracy and F1
scores, even with the worse models it does not seem that the overall
populations are so strikingly different. One could therefore be tempted, in
future studies, from dismissing an analysis of model performance under DS. We
highlight, however, that we are working in a setting where the sources of DS
are minimized (no temporal dependence, no varying external conditions, stable
lake environment), so we expect that in generic plankton classification
settings OOD will be larger than the one we find. Ultimately, in order to
really see whether DS shift is present, and to assess the utility of different
architectures addressing it, one should compare the effect of different time
series on the data that is generated through these classifiers. For example,
one should compare the time series created by the two different models, and
check whether they are quantitatively similar, through the measurement of
autocorrelation functions and related quantities.

##### On the size of the OOD cells

The size of our OOD cells meets a trade-off between data availability, and
expert work. OOD cells of size 1000, with 35 classes, do not ensure a fully
representative distribution of the data, but creating extremely large and
numerous OOD cells is beyond the scope of this work. This reflects on an
increased measurement uncertainty on Hellinger distances, performances, and
sensitivities, which is exacerbated with rare classes. However, when engaging
in a large plankton classification campaign, one must likely regularly
manually check that the classification is stable and to perform active
learning [67]. This ensures the availability of increasing amounts of data
that can be used as OOD cells, thus providing progressively more interpretable
results.

##### On data leakage

The procedure we propose performs model selection based on the OOD cells.
Since the OOD cells are used for model selection, they should not be used to
report performances [68]. In future work, we will retrain our selected model
on a combination of all ID and OOD datasets, and deploy it on a very large
amount of unlabeled data, to construct plankton abundance time series. We will
test the classifications on random days, and these will be used as an
assessment of the performance.

##  Appendix A Related Work

Here, we provide a brief summary of recent literature related to DS in the
context of plankton monitoring.

Ref. [28] studies DS with an AlexNet on the WHOI dataset [14], which contains
millions of annotated black & white images from a Flow Cytobot deployed in a
sea environment. They build a training set plus 21 test cells, which are small
datasets used to assess the OODG, from where they clearly show the presence of
class-distribution shift. They conclude that domain adaptation methods for
quantification outperform random sampling.

Ref. [29] poses emphasis on how to test classifiers in the presence of DS,
since with DS the test performance uses significance. Instead of testing on
hold-out data as done in Ref. [28], Ref. [29] proposes to use a cross-
validation (CV) scheme, where each of 60 CV folds corresponds to a day of
sampling. However, CV is viable because they work with Support Vector Machines
and Random Forests, but it is generally not viable with deep neural networks.
Further, the features that we use for calculating Hellinger distances are also
used as input features, whereas in our case we use pixel values as an input
for the models, and features for calculating Hellinger distances. Similarly to
us, they assess dissimilarity among datasets through the Hellinger Distance
$D_{\mathcal{H}}$, and find a weak correlation between accuracy and
$D_{\mathcal{H}}$.

Ref. [42] studies corrections to quantification in the context of binary
classification of diatom chains, proposing to manually assess the performance
on the target domain, and to use this information to correct the
quantification. However, while viable for binary classification, it is not
clear whether this would also be efficient in the multi-class case, since
reconstructing the probabilities on the target domain would require the manual
annotation of large volumes of images.

One alternative is to only focus on identifying rare classes [30], and
treating the rest as a background. Treating most of the images as a background
allows to use a hard-negative mining scheme, which consists of retraining the
model with images that would belong to the background and were wrongly
classified. However, while useful for specific target studies, this does not
help when one is interested in the entire taxonomic distribution of the target
ecosystem.

Ref. [35] studies out-of-domain generalization with particular emphasis on the
appearing of new classes, and on a radical domain shift such as sampling from
a different monitoring station. For example, they train on data coming from
one station, and test on data coming from another. It does however not address
the DS appearing within a same station, i.e. the difference in performance
from the same station, just because the data are sampled at a different point
in time. In particular, the focus is on differing stations, without specifying
whether the OOD cells belong to a specific day of sampling, thus representing
the data distribution in that particular point in time.

##  Appendix B Data

In Fig. 21, we show a collection of images from each class. For a detailed
description of all the classes, we refer the reader to the appendix of Ref.
[1].

Figure 21: Sample images from each class in the ZooLake2.0 dataset.

##  Appendix C Models

Transfer learning [69] is used to save time and computing resources. We import
the models and weights pre-trained on the ImageNet dataset [70] from the TIMM
Python package [71], and replace the final layer with a dense layer which
connects to $N_{c}$ outputs. The whole model training scheme is organized into
two stages. In the first stage, all layers of the pre-trained model but the
last layer are frozen, only the parameters in the last layer are trained. The
training optimizer is AdamW [72]. The two hyperparameters in AdamW, learning
rate ($\eta$) and weight decay are tuned by Bayesian Optimization search [73],
using the Ray.Tune package [74]. All images of ZooLake2.0 are resized to a
size of $224\times 224$, then 25 trials with different hyperparameter
combinations are trained for 50 epochs, to find the best combination of
hyperparameters by Bayesian Optimization search. The best hyperparameter set
is then saved, and used for the second stage of training. In the second stage,
by unfreezing all layers, all parameters of the model are trained for 100
epochs, with a low learning rate of $\eta=10^{-5}$. During the training, the
real-time training curve is saved after every epoch. The model is saved as a
checkpoint whenever the F1 score on the validation set improves after a new
epoch. Once the entire training is completed, the best model and the last
epoch model are saved. The training of all models are performed in the PyTorch
framework [75], with two Nvidia GTX 2080Ti GPUs.

##  Appendix D Feature description

In order to describe the morphological and color feature of image, 67
descriptors are extracted for each image. The explanations of these feature
descriptors are given below.

width, height

the width and height of bounding rectangle of the object in the image.

w_rot, h_rot

the width and height of rotated bounding rectangle of the minimum area.

angle_rot

the angle between rotated and original bounding boxes, ranging from 0 to 90.

aspect_ratio

the ratio of width to height.

rect_area

the area of bounding box, i.e. the product of width and height.

contour_area, contour_perimeter

the area and perimeter of the contour of the object.

extent

the ratio of contour area to bounding rectangle area.

compactness

the squared contour perimeter, divided by $4\pi\times$ contour area. The
circle has a compactness of 1.

formfactor

the multiplicative inverse of compactness.

hull_area, hull_perimeter

the area and perimeter of the convex hull. Convex hull is the smallest convex
polygon that can fit in the object.

solidity

the ratio of the object area to the convex hull area. This measures the
density of an object.

ESD

equivalent spherical diameter, the diameter of the circle whose area is same
as the contour area.

major_axis, minor_axis

the length of major and minor axis of the fitted ellipse.

angle

the angle between major axis and vertical axis, ranging from 0 to 180.

eccentricity

the ratio of minor axis to major axis.

convexity

the ratio of hull perimeter to contour perimeter.

roundness

the squared hull perimeter, divided by $4\pi\times$ contour area. The circle
has a roundness of 1.

intensity_R_mean, intensity_G_mean, intensity_B_mean

the mean pixel value of 3 color channels respectively.

intensity_R_std, intensity_G_std, intensity_B_std

the standard deviation of pixel values of 3 color channels respectively.

hue_mean, saturation_mean, brightness_mean

the mean value of hue, saturation and brightness respectively.

hue_std, saturation_std, brightness_std

the standard deviation of hue, saturation and brightness respectively.

blurriness

the mean of absolute values of Laplacian.

noise

the standard deviation of all pixel values of an image.

$M_{ij},\mu_{ij},\eta_{ij},I_{ij}$

these are the moments of object contour, up to the third order, as described
in Ref. [76].191919A simple description of these moments can be found at
https://en.wikipedia.org/wiki/Image_moment (visited on Jan. $1^{\mathrm{st}}$,
2024). These moments are weighted sums of the pixels, which can be used to
characterize images. The $M_{ij}$ are raw moments - they do not have any
specific property beyond characterizing the image; the $\mu_{ij}$ are the
central moments - these are translation invariant and can provide information
about the orientation and elongation of the object in the image; the
$\eta_{ij}$ describe the distribution of pixel intensities with respect to the
horizontal and vertical axes, which is invariant to both translation and
scale, they so can be related to the orientation of the image; the $I_{i}$ are
invariant with respect to translation, scale, and rotation, and are therefore
useful for characterizing and recognizing shape in image.

##  Appendix E Targeted Augmentations

Figure 22:  Graphical examples of targeted augmentations compared to original
images.

According to the sensitivity analysis in Sec. III.2.3 the most relevant
features are saturation_std, intensity_G_std, eccentricity, compactness,
formfactor, convexity, hull_area, height, ESD, blurriness, angle_rot and some
image moments, including raw moments, translation invariant moments, scale
invariant moments and rotation invariant moments. These are related with
color, shape, object size, blurriness and orientation. These results do not
highlight a single distinguishable feature, possibly due to the small size and
number of OOD cells, which give us uncertain estimates of the sensitivities.
However, it still indicates that robustness towards changes in shape, color,
size, orientation and blurriness should be addressed. Therefore, we deploy the
following augmentations:

  * •

Color jittering includes the tuning of brightness, contrast, saturation and
hue. This is the targeted augmentation for saturation_std, intensity_G_std.

  * •

Gaussian blur blurs image with randomly chosen Gaussian blur. This is the
targeted augmentation for blurriness.

  * •

Random affine includes the translation and shearing of original images. This
is the targeted augmentation for shape-dependent features, e.g. compactness,
form factor, eccentricity, convexity and moments.

  * •

Random resized crop crops a random portion of image and resizes it to a given
size. This is useful to handle the fluctuation in size-dependent features,
e.g. height, hull_area, ESD, raw moments and translation invariant moments.

  * •

Random perspective performs a random perspective transformation of the given
image with a given probability. This effectively adjusts the shape of the
object and helps to enhance the model robustness against the variation of
shape-dependent features, including compactness, formfactor, eccentricity,
convexity and all moments.

Examples of the augmented image are shown in Fig. 22.

##  Appendix F RGB channel standardization

As we report in Fig. 23, standardizing the RGB channels does not help to
improve the OOD performance of our models.

Figure 23: ID and OOD Performances of a MobileNet model, with and without the
RGB channel standardization.

##  Appendix G Other ways of estimating the distribution shift

While the Hellinger distance is the one that best correlates with the
performance drop in our data, we also try several other dissimilarity measures
as estimator of dataset compositional shift. We report on the following:

  * •

Hellinger distance (described in the main text)

  * •

Wasserstein distance

  * •

Kullback-Leibler divergence

##### Wasserstein Distance

The Wasserstein distance, also known as Earth Mover’s Distance (EMD), is a
mathematical measure of the distance between two probability distributions
[77]. It measures the minimum amount of work required to transform one
distribution into the other, where work is defined as the amount of ”mass”
moved multiplied by the distance it is moved. Given two normalized
distributions $p$ and $q$ defined on a discrete one-dimensional support, the
Wasserstein distance is

$D_{\mathcal{W}}^{f}=\frac{1}{n_{\mathrm{bins}}}\sum_{i=1}^{n_{\mathrm{bins}}}\left|p_{i}-q_{i}\right|.$
(33)

##### Kullback-Leibler Divergence

Kullback-Leibler (KL) divergence, also known as relative entropy, is a measure
of how one probability distribution differs from a second, reference
distribution [78]. Different from the Hellinger and Wasserstein distance, the
Kullback-Leibler divergence is an asymmetric measure.202020So, strictly
speaking, it is not a distance. The divergence from $q$ to $p$
$D_{\mathcal{KL}}^{f}(p\parallel q)$ and the divergence from $p$ to $q$
$D_{\mathcal{KL}}^{f}(q\parallel p)$ are different. For discrete probability
distributions $p$ and $q$ defined on the same sample space $\mathcal{X}$, the
KL divergence from $q$ to $p$ is

$D_{\mathcal{KL}}^{f}(p\parallel
q)=\sum_{x\in\mathcal{X}}p(x)\log{\frac{p(x)}{q(x)}}\,.$ (34)

In our study, we keep the in-domain feature distribution always as reference
distribution, calculate the divergence of out-of-domain distribution to the
reference.

In Fig. 24 and Fig. 25, we show the equivalent of Fig. 10, but using
$D_{\mathcal{W}}$ and $D_{\mathcal{KL}}$. The correlation between distance and
performance drop is lower than using the Hellinger distance.

Figure 24: (a): F1 drop [Eq. (20)] as a function of the Wasserstein distance
from the training set [Eq. (23)]. The error bars are calculated by taking the
fluctuations among the different classes in each test cell. (b): Same, for the
recall.  Figure 25: (a): F1 drop [Eq. (20)] as a function of the Kullback-
Leibler divergence from the training set [Eq. (23)]. The error bars are
calculated by taking the fluctuations among the different classes in each test
cell. (b): Same, for the recall.

##  Appendix H Per-class performance drop

Fig. 26 shows the F1-score of MobileNet model on ID_train, ID_test and OOD
(micro and macro) sets, for all 35 classes. Some classes such as asterionella,
bosmina, leptodora and polyarthra have worse OOD performance compared to other
classes. Fig. 27 displays the similar results, but using the BEsT model. In
general, the OOD performance of each single class improves. Note, also, that
most of the low performances are not associated with living organisms.

Figure 26: Performance drop for each class, with baseline model.  Figure 27:
Performance drop for each class, with BEsT model.

##  Appendix I Performance drop with aggregated classes

The classes in the ZooLake2.0 dataset can be coarsened in superclasses. As a
result, the classes hydra and kellicottia are merged into a new superclass
hydra, classes keratella_cochlearis and keratella_quadrata are combined as
keratella, classes nauplius and paradileptus become a new nauplius, classes
collotheca, conochilus, polyarthra, rotifers, synchaeta and trichocerca are
merged into a single superclass rotifer, classes dirt, unknown and
unknown_plankton go to the superclass unknown. In Fig. 28, we show the OOD
performance drop, when using superclasses.

Figure 28:  Same figure as Fig. 5, but using superclasses.

In Fig. 29, we show the classification performance on the BEsT model using
superclasses. Although all performance metrics show improvement over the
baseline model, the OOD performance still lags far behind that of the ID.

Figure 29:  Same figure as Fig. 28, but using BEsT model.

In Fig. 30, we show the equivalent of Fig. 19c, updating the categories with
superclasses.

Figure 30:  Bar plot showing the true and the predicted plankton abundances,
for the BEsT model, on the aggregated OOD data on superclass level.

##  Appendix J OOD degradation on the single OOD cells

In Fig. 31, we show the population scatter plot for 10 individual OOD cells,
comparing the OOD performance of the baseline MobileNet model with the ID test
performance. The diagonal lines represent the 1:1 relationship between
predicted and true abundances in each OOD cell. In Fig. 32 we show the same
figure, but comparing the MobileNet with the BEsT model.

Figure 31: Same figures as Fig. 7, but for the individual OOD cells. Figure
32: Same figures as Fig. 19a, but for the individual OOD cells.

##  Appendix K Sensitivities for all models and classes

###  K.1 Estimating performance drop on new data

We show the F1 drop ratio as a function of $D_{\mathcal{H}}$ in Fig. 33a, for
a MobileNet model (we also show this in Fig. 34, for the BEsT model). This
could allow us to estimate the F1 drop in a generic OOD dataset. For example,
if a new unlabeled dataset has $D_{\mathcal{H}}\simeq 0.12$, we can expect an
F1 drop of around 10%. This estimate must however take into account the
reduced size of the OOD cells: as highlighted in Sec. II.4, the F1 score tends
to be artificially lower in the presence of small abundances, so we expect the
drop to be lower in the presence of bigger OOD datasets.212121This is, for
example, confirmed by the fact that micro OOD F1 scores (calculated on the
aggregated OOD data) are larger than macro OOD F1 scores (calculated on the
OOD cells one by one, and only later averaged). We see this e.g. in Fig. 14.
Moreover, in this case, since neither the F1 score nor the accuracy drop
correlate well with the Hellinger distance, so more work/data is needed to
confidently estimate the performance degradation with completely new unlabeled
data.

Figure 33: (a) F1 and (b) accuracy drop as a function of the Hellinger
distance between the ID and OOD $P(\mathbf{x})$, for a MobileNet model. An
equivalent figure for the BEsT architecture is provided in Sec. K.1.  Figure
34: (a) F1 drop and (b) accuracy drop as a function of the Hellinger distance
between the ID and OOD $P(\mathbf{x})$, for the BEsT model. Same figure as
Fig. 33, which is for the baseline model.

###  K.2 Sensitivities for all models and classes

In Fig. 35, we show the equivalent of Fig. 10, but for the BEsT model. The
correlation between performance drop and $D_{\mathcal{H}}$ persists also with
a different model.

Figure 35: (a): F1 drop [Eq. (20)] as a function of the Hellinger distance
from the training set [Eq. (23)]. The error bars are calculated by taking the
fluctuations among the different classes in each test cell. (b): Same, for the
recall. These figures are similar to Fig. 10, but using the BEsT model.

In Fig. 36, we show the per-class sensitivities for the BEsT model. This
figure is the equivalent of Fig. 12. Note that, as a result of our pipeline
addressing OOD robustness, the values of the sensitivities are lower, and that
dinobryon is not the most sensitive plankton class anymore.

Figure 36: Similar figures to Fig. 12, but using the BEsT model. (a): per-
class F1 drop as a function of the Hellinger distance. (b): sensitivity
$Q_{F_{1}}$ related to each class of those for which enough data is available.

##  Appendix L Comparison of different architectures on the OOD cells

In Fig. 37 we compare the performance of all the simulated architectures in
each of the ID and OOD sets. There is no single architecture that performs
best in each of the OOD cells.

Figure 37: Comparison of different architectures in each of the ID and OOD
sets.

##  Appendix M Different ways of ensembling

Ensembling can be carried out both by comparing different architectures, and
by comparing different instances of the same architecture [1]. We train 3
models for each architecture, with basic augmentations, and compare both
methodologies.

For the single-architecture ensembling, we join the three trained instances.
For the multi-architecture ensembling, we do not mix CNNs with Vision
Transformers, because they have very different confidence vectors [20]. Hence,
we obtain a single CNN ensemble and a single Transformer ensemble (we call
them E-CNN and E-Trans). In both cases, we take a single model for each
architecture, the one with the best validation F1 score. Note that multi- are
more expensive than single-architecture ensembles, since they use 4 models,
and require the training of 12.

In Fig. 38a, we compare the OOD performance of all the trained ensembles. We
have several remarks:

  * •

Multi-architecture CNN ensembles outperform those with transformers.

  * •

While, depending on the chosen metrics, the best-performing ensemble changes,
we see that the multi-architecture ensembles are not outperforming the single-
architecture ones. However, while E-Trans are not better than their single-
model counterparts, E-CNN offers a big improvement with respect to single-
architecture CNN ensembles. This is consistent with the ID results reported in
Ref. [1].

  * •

The best single-architecture CNN is EfficientNet-B7. This is consistent with
the (in-dataset) results of Ref. [1].

  * •

Single-model transformer ensembles, in particular BEiTs, have very good
performances. The good performance of BEiTs was already reported in Ref. [22],
while that of single-model Transformer ensembles was noted in Ref. [20].

  * •

The best ensembles are EfficientNet-B7, and BEiT. EfficientNets have the best
macro-averaged metrics, while BEiTs have the best micro-averaged metrics. This
essentially indicates that the former perform better with minority classes,
and the latter with majority classes.

In Fig. 38b, we show the performances of all models on each single OOD cell.

Figure 38: Performances of the ensemble models, trained with basic data
augmentation. (a): Micro- and Macro-OOD performances, for each ensemble. (b):
Performances on each OOD cell. Error bars on the macro-OOD are the standard
errors among the single OOD cells.

##  Appendix N More test time augmentations

In Tab. 3, we show a comparison between flipping and rotating. Flipping does
not increase the performance with respect to rotating. The time complexity of
the TTA operation in the last row is doubled, but the model performance is not
improved. This is expected, since the only difference stays in the chirality
of the images, but plankton images are generally not chiral.

Table 3: Comparison of TTA result of with/without additional flipping. The  |  Accuracy $\uparrow$ |  F1-score $\uparrow$ |  BC $\downarrow$ |  NMAE $\downarrow$
---|---|---|---|---
| (Micro-OOD/Macro-OOD)
TTA with 4 rotation angles | 0.829/0.832 | 0.761/0.610 | 0.086/0.100 | 0.261/0.398
TTA with 4 flipped states | 0.827/0.829 | 0.758/0.595 | 0.086/0.102 | 0.274/0.382
TTA with 4 rotation angles and corresponding vertical flipping | 0.829/0.832 | 0.761/0.616 | 0.086/0.100 | 0.262/0.386

## References

  * Kyathanahally  _et al._ [2021a] S. P. Kyathanahally, T. Hardeman, E. Merz, T. Bulas, M. Reyes, P. Isles, F. Pomati, and M. Baity-Jesi, Deep learning classification of lake zooplankton, Frontiers in microbiology , 3226 (2021a), arXiv:2108.05258 .
  * Behrenfeld  _et al._ [2001] M. J. Behrenfeld, J. T. Randerson, C. R. McClain, G. C. Feldman, S. O. Los, C. J. Tucker, P. G. Falkowski, C. B. Field, R. Frouin, W. E. Esaias, D. D. Kolber, and N. H. Pollack, Biospheric primary production during an enso transition, Science 291, 2594 (2001), https://science.sciencemag.org/content/291/5513/2594.full.pdf .
  * Falkowski [2012] P. Falkowski, Ocean science: the power of plankton, Nature 483, S17 (2012).
  * Xu  _et al._ [2001] F.-L. Xu, S. Tao, R. W. Dawson, P.-g. Li, and J. Cao, Lake ecosystem health assessment: indicators and methods, Water Research 35, 3157 (2001).
  * Franks [2002] P. J. Franks, Npz models of plankton dynamics: their construction, coupling to physics, and application, Journal of Oceanography 58, 379 (2002).
  * Ser-Giacomi  _et al._ [2018] E. Ser-Giacomi, L. Zinger, S. Malviya, C. De Vargas, E. Karsenti, C. Bowler, and S. De Monte, Ubiquitous abundance distribution of non-dominant plankton across the global ocean, Nature ecology & evolution 2, 1243 (2018).
  * Hu  _et al._ [2022] J. Hu, D. R. Amor, M. Barbier, G. Bunin, and J. Gore, Emergent phases of ecological diversity and dynamics mapped in microcosms, Science 378, 85 (2022).
  * Zohdi and Abbaspour [2019] E. Zohdi and M. Abbaspour, Harmful algal blooms (red tide): a review of causes, impacts and approaches to monitoring and prediction, International Journal of Environmental Science and Technology 16, 1789 (2019).
  * Benfield  _et al._ [2007] M. C. Benfield, P. Grosjean, P. F. Culverhouse, X. Irigoien, M. E. Sieracki, A. Lopez-Urrutia, H. G. Dam, Q. Hu, C. S. Davis, A. Hansen, C. H. Pilskaln, E. M. Riseman, H. Schultz, P. E. Utgoff, and G. Gorsky, Rapid: Research on automated plankton identification, Oceanography 20, 172 (2007).
  * Cowen and Guigand [2008] R. K. Cowen and C. M. Guigand, In situ ichthyoplankton imaging system (isiis): system design and preliminary results, Limnology and Oceanography: Methods 6, 126 (2008).
  * Campbell  _et al._ [2013] L. Campbell, D. W. Henrichs, R. J. Olson, and H. M. Sosik, Continuous automated imaging-in-flow cytometry for detection and early warning of karenia brevis blooms in the gulf of mexico, Environmental Science and Pollution Research 20, 6896 (2013).
  * Orenstein  _et al._ [2020a] E. C. Orenstein, D. Ratelle, C. Briseño-Avena, M. L. Carter, P. J. S. Franks, J. S. Jaffe, and P. L. D. Roberts, The scripps plankton camera system: A framework and platform for in situ microscopy, Limnology and Oceanography: Methods 18, 681 (2020a), https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.1002/lom3.10394 .
  * Merz  _et al._ [2021] E. Merz, T. Kozakiewicz, M. Reyes, C. Ebi, P. Isles, M. Baity-Jesi, P. Roberts, J. S. Jaffe, S. R. Dennis, T. Hardeman, _et al._ , Underwater dual-magnification imaging for automated lake plankton monitoring, Water Research 203, 117524 (2021), biorxiv:2021.04.14.439767 .
  * Orenstein  _et al._ [2015] E. C. Orenstein, O. Beijbom, E. E. Peacock, and H. M. Sosik, Whoi-plankton- a large scale fine grained visual recognition benchmark dataset for plankton classification (2015), arXiv:1510.00745 [cs.CV] .
  * Py  _et al._ [2016] O. Py, H. Hong, and S. Zhongzhi, Plankton classification with deep convolutional neural networks, in  _2016 IEEE Information Technology, Networking, Electronic and Automation Control Conference_ (2016) pp. 132–136.
  * Ellen  _et al._ [2019] J. S. Ellen, C. A. Graff, and M. D. Ohman, Improving plankton image classification using context metadata, Limnology and Oceanography: Methods 17, 439 (2019).
  * Orenstein and Beijbom [2017] E. C. Orenstein and O. Beijbom, Transfer learning and deep feature extraction for planktonic image data sets, in  _2017 IEEE Winter Conference on Applications of Computer Vision (WACV)_ (2017) pp. 1082–1088.
  * Lumini and Nanni [2019] A. Lumini and L. Nanni, Deep learning and transfer learning features for plankton classification, Ecological Informatics 51, 33 (2019).
  * Eerola  _et al._ [2020] T. Eerola, K. Kraft, O. Grönberg, L. Lensu, S. Suikkanen, J. Seppälä, T. Tamminen, H. Kälviäinen, and H. Haario, Towards operational phytoplankton recognition with automated high-throughput imaging and compact convolutional neural networks, Ocean Science Discussions , 1 (2020).
  * Kyathanahally  _et al._ [2022] S. P. Kyathanahally, T. Hardeman, M. Reyes, E. Merz, T. Bulas, F. Pomati, and M. Baity-Jesi, Ensembles of vision transformers as a new paradigm for automated classification in ecology, Scientific Reports 12, https://doi.org/10.1038/s41598-022-21910-0 (2022), arXiv:2203.01726 .
  * Yue  _et al._ [2023] J. Yue, Z. Chen, Y. Long, K. Cheng, H. Bi, and X. Cheng, Toward efficient deep learning system for in-situ plankton image recognition, Frontiers in Marine Science (2023).
  * Maracani  _et al._ [2023] A. Maracani, V. Pastore, L. Natale, L. Rosasco, and F. Odone, In‑domain versus out‑of‑domain transfer learning in plankton image classification, Scientific Reports 13, 10443 (2023).
  * Eerola  _et al._ [2023] T. Eerola, D. Batrakhanov, N. V. Barazandeh, K. Kraft, L. Haraguchi, L. Lensu, S. Suikkanen, J. Seppälä, T. Tamminen, and H. Kälviäinen, Survey of automatic plankton image recognition: Challenges, existing solutions and future perspectives, arXiv preprint arXiv:2305.11739 (2023).
  * Kenitz  _et al._ [2020] K. M. Kenitz, E. C. Orenstein, P. L. D. Roberts, P. J. S. Franks, J. S. Jaffe, M. L. Carter, and A. D. Barton, Environmental drivers of population variability in colony-forming marine diatoms, Limnology and Oceanography 65, 2515 (2020), https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.1002/lno.11468 .
  * Merkli  _et al._ [2024] S. Merkli, E. Merz, M. Reyes, P. Ntetsika, S. Dennis, P. Isles, S. Kyathanahally, M. Baity-Jesi, S. To, A. Merico, and F. Pomati, Automated plankton monitoring suggests a key role of microzooplankton and temperature for predicting dynamics of phytoplankton size classes (2024).
  * Solow  _et al._ [2001] A. Solow, C. Davis, and Q. Hu, Estimating the taxonomic composition of a sample when individuals are classified with error, Marine Ecology Progress Series 216, 309 (2001).
  * Bell and Hopcroft [2008] J. L. Bell and R. R. Hopcroft, Assessment of zooimage as a tool for the classification of zooplankton, Journal of Plankton Research 30, 1351 (2008).
  * Beijbom  _et al._ [2015] O. Beijbom, J. Hoffman, E. Yao, T. Darrell, A. Rodriguez-Ramirez, M. Gonzalez-Rivero, and O. H. Guldberg, Quantification in-the-wild: data-sets and baselines, in  _NIPS 2015 Workshop on Transfer and Multi-Task Learning, Montreal, CA_ (2015) arXiv:1510.04811 .
  * González  _et al._ [2017a] P. González, E. Álvarez, J. Díez, Á. López-Urrutia, and J. J. del Coz, Validation methods for plankton image classification systems, Limnology and Oceanography: Methods 15, 221 (2017a).
  * Walker and Orenstein [2021] J. L. Walker and E. C. Orenstein, Improving rare-class recognition of marine plankton with hard negative mining, in  _Proceedings of the IEEE/CVF International Conference on Computer Vision_ (2021) pp. 3672–3682.
  * Quiñonero-Candela  _et al._ [2009] J. Quiñonero-Candela, M. Sugiyama, N. D. Lawrence, and A. Schwaighofer,  _Dataset shift in machine learning_ (Mit Press, 2009).
  * Moreno-Torres  _et al._ [2012] J. G. Moreno-Torres, T. Raeder, R. Alaiz-Rodríguez, N. V. Chawla, and F. Herrera, A unifying view on dataset shift in classification, Pattern Recognition 45, 521 (2012).
  * Saadati  _et al._ [2023] M. Saadati, A. Balu, S. Chiranjeevi, T. Z. Jubery, A. K. Singh, S. Sarkar, A. Singh, and B. Ganapathysubramanian, Out-of-distribution detection algorithms for robust insect classification, arXiv preprint arXiv:2305.01823 (2023).
  * Ellen  _et al._ [2015] J. Ellen, H. Li, and M. D. Ohman, Quantifying california current plankton samples with efficient machine learning techniques, in  _OCEANS 2015 - MTS/IEEE Washington_ (2015) pp. 1–9.
  * Yang  _et al._ [2022] Z. Yang, J. Li, T. Chen, Y. Pu, and Z. Feng, Contrastive learning-based image retrieval for automatic recognition of in situ marine plankton images, ICES Journal of Marine Science 79, 2643 (2022).
  * Wang  _et al._ [2022] J. Wang, C. Lan, C. Liu, Y. Ouyang, T. Qin, W. Lu, Y. Chen, W. Zeng, and P. Yu, Generalizing to unseen domains: A survey on domain generalization, IEEE Transactions on Knowledge and Data Engineering (2022).
  * Peng  _et al._ [2022] X. Peng, F. Qiao, and L. Zhao, Out-of-domain generalization from a single source: An uncertainty quantification approach, IEEE Transactions on Pattern Analysis and Machine Intelligence (2022).
  * Ouyang  _et al._ [2022] C. Ouyang, C. Chen, S. Li, Z. Li, C. Qin, W. Bai, and D. Rueckert, Causality-inspired single-source domain generalization for medical image segmentation, IEEE Transactions on Medical Imaging 42, 1095 (2022).
  * Badreldeen Bdawy Mohamed  _et al._ [2022] O. Badreldeen Bdawy Mohamed, T. Eerola, K. Kraft, L. Lensu, and H. Kälviäinen, Open-set plankton recognition using similarity learning, in  _International Symposium on Visual Computing_ (Springer, 2022) pp. 174–183.
  * González  _et al._ [2017b] P. González, A. Castaño, N. V. Chawla, and J. J. D. Coz, A review on quantification learning, ACM Computing Surveys (CSUR) 50, 1 (2017b).
  * González  _et al._ [2019] P. González, A. Castaño, E. E. Peacock, J. Díez, J. J. Del Coz, and H. M. Sosik, Automatic plankton quantification using deep features, Journal of Plankton Research 41, 449 (2019).
  * Orenstein  _et al._ [2020b] E. C. Orenstein, K. M. Kenitz, P. L. Roberts, P. J. Franks, J. S. Jaffe, and A. D. Barton, Semi-and fully supervised quantification techniques to improve population estimates from machine classifiers, Limnology and Oceanography: Methods 18, 739 (2020b).
  * Kyathanahally  _et al._ [2021b] S. Kyathanahally, T. Hardeman, E. Merz, T. Kozakiewicz, M. Reyes, P. Isles, F. Pomati, and M. Baity-Jesi, Data for: Deep learning classification of lake zooplankton (2021b).
  * Schuer  _et al._ [2023] C. Schuer, L. Gasser, F. Perez-Cruz, K. Schirmer, and M. Baity-Jesi, A benchmark dataset for machine learning in ecotoxicology, Scientific Data 10, 718 (2023), bioRxiv:2023.05.27.542160 .
  * Howard  _et al._ [2019] A. Howard, M. Sandler, G. Chu, L.-C. Chen, B. Chen, M. Tan, W. Wang, Y. Zhu, R. Pang, V. Vasudevan, _et al._ , Searching for mobilenetv3, in  _Proceedings of the IEEE/CVF international conference on computer vision_ (2019) pp. 1314–1324.
  * Huang  _et al._ [2017] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, Densely connected convolutional networks, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2017) pp. 4700–4708.
  * Tan and Le [2019] M. Tan and Q. Le, Efficientnet: Rethinking model scaling for convolutional neural networks, in  _International conference on machine learning_ (PMLR, 2019) pp. 6105–6114.
  * Touvron  _et al._ [2020] H. Touvron, M. Cord, M. Douze, F. Massa, A. Sablayrolles, and H. Jégou, Training data-efficient image transformers & distillation through attention, CoRR abs/2012.12877 (2020), 2012.12877 .
  * Dosovitskiy  _et al._ [2020] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, _et al._ , An image is worth 16x16 words: Transformers for image recognition at scale, arXiv preprint arXiv:2010.11929 (2020).
  * Bao  _et al._ [2021] H. Bao, L. Dong, S. Piao, and F. Wei, Beit: Bert pre-training of image transformers, arXiv preprint arXiv:2106.08254 (2021).
  * Liu  _et al._ [2021] Z. Liu, Y. Lin, Y. Cao, H. Hu, Y. Wei, Z. Zhang, S. Lin, and B. Guo, Swin transformer: Hierarchical vision transformer using shifted windows, in  _Proceedings of the IEEE/CVF international conference on computer vision_ (2021) pp. 10012–10022.
  * Tan  _et al._ [2018a] C. Tan, F. Sun, T. Kong, W. Zhang, C. Yang, and C. Liu, A survey on deep transfer learning, in  _International conference on artificial neural networks_ (Springer, 2018) pp. 270–279.
  * González-Castro  _et al._ [2013] V. González-Castro, R. Alaiz-Rodríguez, and E. Alegre, Class distribution estimation based on the hellinger distance, Information Sciences 218, 146 (2013).
  * Tortorelli and Michaleris [1994] D. A. Tortorelli and P. Michaleris, Design sensitivity analysis: overview and review, Inverse problems in Engineering 1, 71 (1994).
  * Huang [1987] K. Huang,  _Statistical Mechanics_ , 2nd ed. (John Wiley and Sons, Hoboken, NJ, 1987).
  * Hendrycks  _et al._ [2021] D. Hendrycks, S. Basart, N. Mu, S. Kadavath, F. Wang, E. Dorundo, R. Desai, T. Zhu, S. Parajuli, M. Guo, D. Song, J. Steinhardt, and J. Gilmer, The many faces of robustness: A critical analysis of out-of-distribution generalization, in  _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ (2021) pp. 8340–8349.
  * Paul and Chen [2022] S. Paul and P.-Y. Chen, Vision transformers are robust learners, in  _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 36 (2022) pp. 2071–2081.
  * Krizhevsky  _et al._ [2012] A. Krizhevsky, I. Sutskever, and G. E. Hinton, Imagenet classification with deep convolutional neural networks, in  _Advances in neural information processing systems_ (2012) pp. 1097–1105.
  * Forman [2008] G. Forman, Quantifying counts and costs via classification, Data Mining and Knowledge Discovery 17, 164 (2008).
  * Bella  _et al._ [2010] A. Bella, C. Ferri, J. Hernández-Orallo, and M. J. Ramirez-Quintana, Quantification via probability estimators, in  _2010 IEEE International Conference on Data Mining_ (IEEE, 2010) pp. 737–742.
  * Touvron  _et al._ [2021] H. Touvron, M. Cord, M. Douze, F. Massa, A. Sablayrolles, and H. Jégou, Training data-efficient image transformers & distillation through attention, in  _International conference on machine learning_ (PMLR, 2021) pp. 10347–10357.
  * Recht  _et al._ [2018] B. Recht, R. Roelofs, L. Schmidt, and V. Shankar, Do cifar-10 classifiers generalize to cifar-10?, arXiv:1806.00451 (2018).
  * Recht  _et al._ [2019] B. Recht, R. Roelofs, L. Schmidt, and V. Shankar, Do imagenet classifiers generalize to imagenet?, in  _International Conference on Machine Learning_ (PMLR, 2019) pp. 5389–5400.
  * Wang and Deng [2018] M. Wang and W. Deng, Deep visual domain adaptation: A survey, Neurocomputing 312, 135 (2018).
  * Patel  _et al._ [2015] V. M. Patel, R. Gopalan, R. Li, and R. Chellappa, Visual domain adaptation: A survey of recent advances, IEEE signal processing magazine 32, 53 (2015).
  * Mo and Zhao [2024] H. Mo and G. Zhao, Ric-cnn: rotation-invariant coordinate convolutional neural network, Pattern Recognition 146, 109994 (2024).
  * Bochinski  _et al._ [2019] E. Bochinski, G. Bacha, V. Eiselein, T. J. W. Walles, J. C. Nejstgaard, and T. Sikora, Deep active learning for in situ plankton classification, in  _Pattern Recognition and Information Forensics_, edited by Z. Zhang, D. Suter, Y. Tian, A. Branzan Albu, N. Sidère, and H. Jair Escalante (Springer International Publishing, Cham, 2019) pp. 5–15.
  * Kapoor and Narayanan [2022] S. Kapoor and A. Narayanan, Leakage and the reproducibility crisis in ml-based science (2022).
  * Tan  _et al._ [2018b] C. Tan, F. Sun, T. Kong, W. Zhang, C. Yang, and C. Liu, A survey on deep transfer learning, in  _Artificial Neural Networks and Machine Learning–ICANN 2018: 27th International Conference on Artificial Neural Networks, Rhodes, Greece, October 4-7, 2018, Proceedings, Part III 27_ (Springer, 2018) pp. 270–279.
  * Deng  _et al._ [2009] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, Imagenet: A large-scale hierarchical image database, in  _2009 IEEE conference on computer vision and pattern recognition_ (Ieee, 2009) pp. 248–255.
  * Wightman [2019] R. Wightman, Pytorch image models, https://github.com/rwightman/pytorch-image-models (2019).
  * Loshchilov and Hutter [2017] I. Loshchilov and F. Hutter, Decoupled weight decay regularization, arXiv preprint arXiv:1711.05101 (2017).
  * Snoek  _et al._ [2012] J. Snoek, H. Larochelle, and R. P. Adams, Practical bayesian optimization of machine learning algorithms, Advances in neural information processing systems 25 (2012).
  * Liaw  _et al._ [2018] R. Liaw, E. Liang, R. Nishihara, P. Moritz, J. E. Gonzalez, and I. Stoica, Tune: A research platform for distributed model selection and training, arXiv preprint arXiv:1807.05118 (2018).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems 32 (2019).
  * Hu [1962] M.-K. Hu, Visual pattern recognition by moment invariants, IRE transactions on information theory 8, 179 (1962).
  * Kantorovich [1960] L. V. Kantorovich, Mathematical methods of organizing and planning production, Management science 6, 366 (1960).
  * Kullback and Leibler [1951] S. Kullback and R. A. Leibler, On information and sufficiency, The annals of mathematical statistics 22, 79 (1951).

## Acknowledgments

This project was funded by the Eawag DF project Big-Data Workflow
(#5221.00492.999.01), the Swiss Federal Office for the Environment (contract
Nr Q392-1149), and the Swiss National Science Foundation (project 182124).