# Physics-assisted machine learning for THz spectroscopy: sensing moisture on
plant leaves

Milan Koumans  Daan Meulendijks  Haiko Middeljans  Djero Peeters  Jacob C.
Douma  Centre for Crop System Analysis, Wageningen University, 6700 AK
Wageningen, the Netherlands  Dook van Mechelen  correspondence to:
<EMAIL_ADDRESS>
###### Abstract

Signal processing techniques are of vital importance to bring THz spectroscopy
to a maturity level to reach practical applications. In this work, we
illustrate the use of machine learning techniques for THz time-domain
spectroscopy assisted by domain knowledge based on light-matter interactions.
We aim at the potential agriculture application to determine the amount of
free water on plant leaves, so-called leaf wetness. This quantity is important
for understanding and predicting plant diseases that need leaf wetness for
disease development. The overall transmission of a moist plant leaf for 12,000
distinct water patterns was experimentally acquired using THz time-domain
spectroscopy. We report on key insights of applying decision trees and
convolutional neural networks to the data using physics-motivated choices.
Eventually, we discuss the generalizability of these models to determine leaf
wetness after testing them on cases with increasing deviations from the
training set.

## Introduction

The quest to implement societally and industrially relevant applications of
THz technology is impeded by aspects such as costs and performance as compared
to alternatives. Despite the wide range of exploratory studies of THz
technology since the late 1980s, and the numerous suggested applications,[1]
very few products exist that use THz-based techniques. Often demonstrators
underperform compared to a cheaper alternative that already exists or that is
conveniently adapted to a new application. However, since THz technology is
rather new on the market, costs will remain high at least for a while. The
only promising way out is to find a novel application for the technology,
which it can almost uniquely serve, and for which the business case is strong
enough to support the high costs.[2] In this case, its overall performance
needs to justify its usage. On the hardware side, technological maturity has
strongly improved over the last decades, although the progress at an
integrated level is lagging.[3] On the software side, despite crucial
advancements such as model-based analysis that is now widely employed,[4] the
limitations are generality, robustness, and speed, which are essential for
realistic application cases.

The agriculture sector is an area where a multitude of sensing technologies
are employed to aid management decisions. Certainly due to this reason, there
is a vivid interest in searching for applications in this field where THz
spectroscopy can make a difference.[5] Globally, pests and pathogens are a big
threat to crop production, with yield losses reported in the range of $9-21$
%.[6] For some pathogens, such as water molds and some fungi, the presence or
absence of free water on the surface of leaves, so-called leaf wetness, is key
for infection and/or sporulation and is therefore an important parameter in
disease epidemiology. A famous example is Phytophthora infestans, the causal
agent of potato late blight, responsible for the Irish Potato Famine in the
mid-19th century. The development of P. infestans depends on the presence of
leaf wetness and the surrounding temperature.[7, 8] In ideal circumstances, P.
infestans can decimate a potato crop in less than 10 days. Control of late
blight, as well as that of other pests and pathogens, is nowadays mostly done
using crop protection products.[9] There is, however, a strong push from
policymakers to reduce this. Early detection and improved predictions of when
and where diseases may be expected can help targeted (preventive) measures.
Instead of directly detecting the molds, which in the field are difficult to
observe, predicting leaf wetness during the growth season is an important
input for decision support systems to advise on spraying crop protection
products.

Terahertz spectroscopy is particularly suited to accurately sense little
amounts of water and is proven to probe leaf properties.[10, 11, 12, 13] Other
technologies that can sense water, such as visible imaging, will have
difficulties estimating very small amounts of water due to a lack of contrast.
Currently used electric measurements are very local, not directly performed on
a plant leaf, and lack interpretability.[14] Although several strong water
absorptions are situated in the infrared spectral range, its frequency domain
operation will make it difficult to distinguish between surface and content
water of a leaf, and its shorter wavelength also causes more local sensing.

To accurately extract parameters from a canopy of plant leaves, which may
move, are curved, have an anisotropic structure, and are in an environment
that is determined by weather conditions, a suitable signal processing method
is needed. Model-based signal processing, e.g., based on the transfer matrix
method, can be powerful in well-behaved conditions, such as paint layer
inspection in an automotive paint shop and wafer metrology in a cleanroom,[15]
it will not describe well the complex situation at hand. Signal processing
using a data-based learning method could be the solution in this case. Many
other studies have used machine learning on THz spectroscopic data.[16, 17,
18, 19] However, most of these studies suffer from a lack of transparency of
the used method such that the quality of the result is unclear in aspects of
generality and reproducibility.

Here, we study the application of machine learning models to THz spectroscopic
data of moist leaves to predict leaf wetness. The primary focus of this paper
is to present a clear and concise approach for applying decision trees and
convolutional neural networks to THz data. Based on the light-matter
interaction, we motivate the architecture of the used models and give
particular emphasis to feature engineering. Both methods accurately predict
the amount of water on a leaf, independent of the droplet pattern. Eventually,
we discuss the influence of model variability and deduce their
generalizability.

Figure 1: (a) Experimental setup showing a leaf on a scale in a THz time-
domain transmission configuration. (b) Selected images from a measurement
series showing the formation of a droplet pattern with weight $g=2.5$, $12.5$
and $22.5$ mg, respectively. (c) Distribution of the experimental data
$\mathbf{E}(t)$ categorized by the total gravimetric weight comprised of the
plasticized leaf and water pattern weight $g$. (d) Distribution of the same
data as in panel (a) but here categorized by $g$.

## Results

The transmitted electric field $E(t)$ of water patterns on plasticized plant
leaves has been experimentally recorded in transmission geometry at THz
frequencies ($0.1-3$ THz) in the time-domain for about 12,000 distinct
patterns as described in the methods (Fig. 1a,b). We recorded two data sets,
$\mathbf{E}(t)$ for which the water patterns are deposited on the top side of
the leaf, and $\mathbf{E^{\prime}}(t)$ for which the water patterns are on the
bottom side of the leaf. For $\mathbf{E}(t)$, the distribution of the total
sample weight, composed of the leaf and the water pattern on top, manifests a
roughly normal distribution (Fig. 1c), whereas the distribution of only the
water weight $g$ is predominant for low weights (Fig. 1d). The maximum value
of $g$ is given by the condition where droplets run down the leaf. Fig. 2a
shows a given measurement series $E(t)\in\mathbf{E}(t)$ for selected values of
$g$, which visually drastically changes with $g$. Besides, each $E(t)$ is also
related to a given water pattern with its droplet size distribution.
Therefore, spectra $E(t)$ with identical value $g$ may mutually differ. The
standard deviation $\sigma(t)$ of $\mathbf{E}(t)$ demonstrates that the
influence of water on $E(t)$ is predominantly present in the range (4, 5) ps,
related to a single passage of THz radiation through the droplet pattern and
leaf (Fig. 2b). However, maybe of more importance for sensing water is the
range of $E(t)$ where the radiation internally reflects within the water
droplets. From simulating $E(t)$ using a transfer matrix model, we conclude
that the first internal reflection inside droplets shows up in $E(t)$ as a
shoulder right after the largest positive peak, that is, around 5.5 ps. This
can be better visualized through the quantity
$\xi=E_{i}(t)/\max(E(t))-E_{0}(t)/\max(E_{0}(t))$ with $i$ the number of
acquisitions within a measurement series, and where $i=0$ stands for a leaf
with $g=0$ (see Fig. 2c). When plotting $\xi$ for the entire data set
$\mathbf{E}(t)$ as a function of $g$ (Fig. 2d), this effect shows up a fork
for low $g$ around 5.5 ps. The ray trajectory in the droplet-plastic-leaf-
plastic system is, however, sufficiently complex that more subtle features in
$E(t)$ cannot be uniquely attributed to a specific radiation path. To obtain
an accurate prediction of the leaf wetness $g$, we employ two different data-
driven methods, decision tree regression and convolutional neural networks.
For each algorithm, we determine the mean absolute error and median percentage
difference, defined as the median of $|g_{p}-g_{b}|/(g_{b}+\epsilon)$ with
predicted weight value $g_{p}$, benchmark weight value $g_{b}$, and small
$\epsilon$ for stability when $g_{b}\approx 0$.

Figure 2: (a) Experimentally determined $E(t)$ for selected values of $g$ in
ambient conditions (see methods). (b) Standard deviation $\sigma(t)$ of
$\mathbf{E}(t)$. (c) $\xi=E(t)/\max(E(t))-E_{0}(t)/\max(E_{0}(t))$ with
$E_{0}(t)$ the transmission for $g=0$. This quantity provides detailed insight
into temporal ranges of large variation. (d) $\mathbf{\xi}$ vs. $g$. The inset
highlights the temporal shift around 5 ps as a function of $g$.

### Decision Trees

The above-described problem to predict $g_{p}$ from the experimental input
matrix $\mathbf{E}(t)$ is equivalent to determining a function $f$ for which
holds that $f:\mathbf{E}(t)\rightarrow g$. We have ascertained that the
relation between distinct features of $\mathbf{E}(t)$, which will be discussed
below, and the target variable $g$ is not linear. To evaluate this problem, we
choose a decision tree method which is proven to work with non-linear mappings
$f$ for mid-sized data sets such as that of $\mathbf{E}(t)$.

A decision tree is a supervised machine learning method that is shaped as a
flowchart in which iterative decisions lead to a piece-wise approximation of
the target variable (see Fig. 3a for a single class and two features $x_{1}$
and $x_{2}$). Decisions within regression trees use quantitative split
criteria such as the absolute error (L1-norm) or squared error (L2-norm) to
determine the best split at a decision point D. A single decision tree is
deterministic in its prediction, making it prone to overfitting. To overcome
this shortcoming, an ensemble of simplified trees is used. Each simplified
tree independently predicts the target value, after which the predictions are
averaged. This reduces the variance of the model while maintaining the
predictive power of the decision trees. Simplifying the trees by limiting the
seen input data is called bagging and is a well-proven method that is capable
of learning nonlinear relations between input features and the target
variable. Due to the simplistic nature of decision trees, inference is
interpretative and most importantly, ensemble modeling leads to improved
robustness compared to single estimators. As common for machine learning
algorithms, decision trees fit on features that are reminiscent of the data.
Many reports employing machine learning methods on THz time-domain data either
first convert to the frequency domain, thereby losing crucial spatial
information, or simply use the time trace $E(t)$.[20, 17, 18] However, points
within $E(t)$ are time-correlated and when directly used as features, the
correlation may hamper determination of the independent effect of each feature
on the target variable. We therefore use feature engineering to extract the
most reminiscent features from $E(t)$, which are not temporally correlated. We
hereto fit a polynomial function of degree $n$ to relevant temporal ranges of
$E(t)$ (Fig. 3b), as further detailed in the methods. $n$ is chosen as small
as possible, but large enough to capture reminiscent features of $E(t)$ to
best predict $g$. Although visually the fit may not seem optimal, for the
algorithm it is. Besides the coefficients of the polynomial terms, also the
beginning of the time window $t_{\text{start}}$ is added as a feature as well
as the absolute air humidity $a$, which has a strong influence on $E(t)$.
Optimization for the displayed range in Fig. 3b leads to $n=11$ as described
in the methods. The feature vector thus reads [$t^{0}$, $t^{1}$, …, $t^{11}$,
$t_{\text{start}}$, $a$], where $t^{0}$ is the bias term and will be further
indicated like that. Among these features there can be, however, ones that are
mutually correlated and ones that are only a little related to the target
variable. We hereto employ recurrent feature elimination to reduce the
dimensionality of the feature matrix and enhance the performance of the model,
as further detailed in the methods. Eventually, hyperparameters responsible
for the regularization of the final ensemble model are determined using 5-fold
cross-validation on a training set consisting of 85 % of the total data set
$\mathbf{E}(t)$.

The predictive performance of a bagged decision tree as described here and
further detailed in the methods was evaluated on an unseen test set consisting
of 15 % of the total data set $\mathbf{E}(t)$ and is shown in Fig. 5a. The
mean absolute error on $g_{p}$ is $0.35^{+0.17}_{-0.1}$ mg and the median
percentage difference is $3.4^{+1.5}_{-1}\%$. The mean inference time is 139
ms ($\sigma=22$ ms) per sample using the hardware as mentioned in the methods.
The indicated error bars are motivated in the Discussion.

Figure 3: (a) Schematic view of the mechanism of a decision tree. The dots
represent data containing features $x_{1}$ and $x_{2}$, and the black lines
show the division of the parameter space according to split criteria Ds with
$s=1..4$. The predicted value of each section is the average of the benchmark
values in the corresponding division. (b) Experimentally determined $E(t)$ for
selected values of $g$ (circles) together with polynomial fits with $n=11$
(solid lines). Figure 4: Activation of $E(t)$ per convolutional layer $\ell$
for (a-f) $g=4$ mg and (g-l) $g=15$ mg averaged over all feature maps (panels
(a,g) correspond to $\ell=1$ etc.). (m) Architecture of the convolutional
neural network, indicating the output dimensions for each layer. The
dimensions after each max-pooling layer are given as $1\times d_{\ell}\times
k_{\ell}$ where $d_{\ell}$ is the length of layer $\ell$ and $k_{\ell}$ the
number of filters of layer $\ell$.

### Convolutional Neural Network

An artificial neural network is a convenient tool to autonomously discover
intricate patterns and representations of signals. The simplest neural network
predicts an output variable $y_{i}$ given input data $X_{train}$ and known
target variable $y_{train}$ as in linear regression by learning a mapping
function $f:X_{train}\rightarrow y_{train}$. The architecture of such a neural
network consists of an input layer with input data $x_{i}\in X_{train}$ and an
output layer with output data $y_{i}$. The real strength of a neural network
is, however, to find patterns and representations in high-dimensional input
data. Hereto, cross-correlations between the inputs are learned using so-
called perceptrons (artificial neurons) that output a latent variable $h$
depending on input $x_{i}$. Subsequently, an activation function, for instance
$\varphi(h)=\max(0,h)$, is applied to $h$. This has the effect that only
neuron output of sufficient importance is fed deeper into the network and
correlated to the output of neurons in the subsequent layer. In this way, the
network can learn non-linear relationships beyond the simple perceptron model.
In case the input can be represented as an image, patterns are more
effectively learned by using a so-called convolutional neural network (CNN).
In such a network a small kernel matrix is scanned over the image matrix to
learn local relations between data points. For pictures, these kernel matrices
can represent lines and circles, but also more complicated patterns, which are
learned by the network. In addition, the method also works for correlated
1-dimensional signals.

The aim here is to construct a CNN, train it on the experimentally acquired
time-domain data $\mathbf{E}(t)$, and compare its performance to predict $g$
to prediction using decision trees as discussed in the previous section. For
the few studies where THz time-domain data is used for CNNs, rarely the full
functional $E(t)$ is used despite its spectral richness and possibility to
effectively augment the number of samples as mentioned in the Discussion.[21,
22, 19] A typical CNN consists of many layers which can be grouped into a
feature extraction part and a regression part. The feature extraction part
consists of multiple convolutional layers $\ell$ that utilize kernel
operations to convolve over the input vector, which enables the detection of
pertinent features in signal $E(t)$. Early layers rather reveal the local
context of signal $E(t)$, whereas deeper layers combine activations of
different temporal regions to extract the global context of signal $E(t)$. The
patterns of the kernels are determined through backpropagation as described in
the methods. For each layer, the number of distinct kernels $k$, also called
filters, has been empirically chosen as a function of the size of $E(t)$, the
complexity of detected patterns in $E(t)$, and the condition to keep the
network size as small as possible. Hence, to capture basic patterns, layer
$\ell=1$ of our CNN has four filters, each of dimension $(1\times 3)$ applied
to each $E(t)$ of dimension $(1\times 760)$ (see Fig. 4m). This results in a
feature map of dimension $(1\times 760\times 4)$. Fig. 4a,g displays the
activation of $E(t)$ for $g=4$ and 15 mg, respectively, that is, the average
of the four filters projected onto $t$. Similar to model-based signal
processing of $E(t)$,[23] CNN shows the largest activation in the range around
the absolute minimum of $E(t)$ and hence demonstrates the local character of
this first layer. Batch normalization has been applied after each
convolutional layer, which normalizes the features during forward propagation.
This ensures activation throughout a deep network and results in improved
convergence while it simultaneously works as a regularizer.[24] A so-called
max-pooling layer is inserted behind each normalization layer to reduce the
network size by downsampling the resulting feature maps. Conceptually, this
layer structure increases the receptive field of each neuron. $\ell$ should be
chosen such that the receptive field of neurons in the last layer covers the
complete input signal $E(t)$. In our case, we empirically evaluated that
$\ell=6$. We can verify the feature extraction performance of this
architecture by inspecting the activations for two cases of $g=4$ and 15 mg
(see Fig. 4a-l). With increasing $\ell$, the activation loses its local
character and spreads throughout the entire temporal range. The concomitant
widening and increasing height of the block functions is the effect of the
max-pooling operations, where every iteration halves the time window and
eliminates the smallest values. The activation shows that for earlier layers,
the range (4, 5) ps is of main importance for the network. This is in
agreement with the earlier observation that the largest amplitude of $E(t)$
varies most with $g$ (cf. Fig. 2a), and as such is the most basic pattern of
$E(t)$. Interestingly, for $\ell=3,4$ and 5, the activations show increased
values also in the range (5,7) ps. As mentioned before, for $g<7$ mg first
internal reflections within the droplets occur around 5.5 ps. In addition,
from modeling $E(t)$ using the transfer matrix method, we find that the first
internal reflection within the leaf material, although strongly damped, occurs
around 7 ps. It thus turns out that these regions which have an increased
importance from a light-matter perspective are likewise important for a CNN.
Moreover, the increased activity in other regions makes our CNN sensitive to
details of $E(t)$ that can contain aspects that are difficult to incorporate
into a physical model.

Figure 5: Predicted weight $g_{p}$ using (a) bagged decision trees and (b) a
convolutional neural network, both versus the benchmark weight $g_{b}$. The
insets display the distribution of the absolute error of $g_{p}$ for both
models.

After feature extraction, regression is performed by a fully connected
artificial neural network. Hereto, all previous activation layers are
flattened to a 1D vector consisting of a concatenation of the 64 $1\times 11$
activations as well as the absolute humidity $\mathbf{a}$. The fully connected
network consisting of two layers subsequently learns the nonlinear relation
between determined input features and the target variable by minimizing its
loss function given by the squared error $|\mathbf{g}_{b}-\mathbf{g}_{p}|^{2}$
using gradient descent optimization.

The CNN is trained on $\mathbf{E}(t)$ for 300 epochs with a batch size of 128
using a train-test ratio of 85:15, where an epoch is a single iteration for
which the complete training set has been (forward and backward) propagated
through the network. During training, a validation set of $10\%$ of the
training set is used to validate the performance after each epoch. Fig. 5b
graphically shows the performance of the model, having a mean absolute error
of $0.38\,^{+0.17}_{-0.1}$mg and a median percentage difference of
$4.1^{+1.5}_{-1}\%$. The mean inference time is 70 ms ($\sigma=9$ ms) per
sample using the hardware as mentioned in the methods.

Figure 6: Predicted weight $g_{p}$ vs. benchmark weight $g_{b}$ using (a,b,e)
bagged decision trees and (c,d,g) CNN, for three test cases. For each case,
the test data set is mentioned in the figure title. Case I tests
generalizability towards unseen water patterns using two unseen measurement
series from $\mathbf{E}(t)$. Case II tests generalizability towards water
patterns on top and below a leaf using the enlarged data set
$\mathbf{E}(t)\cup\mathbf{E^{\prime}}(t)$. Case III tests generalizability
towards water patterns on an unseen leaf surface by training on
$\mathbf{E}(t)$ and testing on $\mathbf{E^{\prime}}(t)$). The inset displays
$E(t)$ and $E^{\prime}(t)$, both at $g=0$, indicating the difference in the
transmitted electric field between irradiating a leaf from the top side and
the bottom side, respectively. $\Delta=g_{b}-g_{p}$ for the 5 data sets
(indicated by the color scheme) of $\mathbf{E}^{\prime}(t)$ vs. the number of
acquisitions $i$ for (f) bagged decision trees and (h) CNN.

## Discussion

The presented results demonstrate that using two very different data-driven
approaches, leaf wetness can be determined using THz spectroscopy within a
confidence interval of about 3 to 4%, independent of the method. This is well
beyond the performance of model-based analysis to determine leaf wetness using
THz spectroscopy.[23] Upon comparison to methodologies reported for the
closely related property of leaf water content,[25, 12, 10] we estimate that
our study, with its large data set and comprehensible approach, gives a
realistic view of the chances of performing quality control on plants.

Although most studies on machine learning conclude at this point, we are here
interested in the generalizability of the models. Therefore, we consider the
following further test cases. For case I, we use 37 out of the 39 measurement
series of dataset $\mathbf{E}(t)$, train the learning methods according to the
same procedure as described before, and test on the two remaining series. Both
methods underperform as compared to the results shown in Fig. 5, with decision
trees having a median percentage difference of 8.8 % and CNN 6.7 % (see Fig.
6(a,c)). The question thus poses about the origin of this decreased
performance. The shape of $E(t)$ is both determined by $g$ and the water
patterns, as mentioned before. However, the temporal range (4, 5) ps with the
largest amplitudes is of leading importance for predicting $g$, and
$\min_{t}\mathbf{E}(t)\approx 5$ ps manifests a monotonous variation with $g$
(see the inset of Fig. 2d), independent of the water patterns, both suggesting
that $g$ predominantly determines $E(t)$. Yet, the values of $g$ of the test
set are a subset of the range of $g$ on which the models are trained. The
effect of water patterns on $E(t)$, on the other hand, is more subtle and
beyond the range of the largest amplitudes, as discussed before. We conjecture
that the test set contains water patterns, rather than different values of
$g$, that are unseen to both methods. The relative underperformance of
decision trees is then probably directly linked to its mechanism for which a
slight variation in $E(t)$, here due to a different pattern, alters the
polynomial fits and thus features on which it has been trained. CNN on the
other hand directly learns regions in $E(t)$ which are of relevance for
learning $g$ for all trained water patterns. Case II considers $\mathbf{E}(t)$
as well as an additional 1501 time traces $E(t)$, referred to as
$\mathbf{E^{\prime}}(t)\notin\mathbf{E}(t)$ where the moisture pattern is
created on the bottom side of the leaf. Instead of shining on the smooth and
reflective top surface as for $\mathbf{E}(t)$, for $\mathbf{E^{\prime}}(t)$
the beam now enters the leaf through the dull and rough bottom surface,
thereby significantly changing $E(t)$ (see inset of Fig. 6). Please note that
the water pattern is always on the emitter side. The methods are trained on a
subset of $\mathbf{E}(t)\cup\mathbf{E^{\prime}}(t)$ and tested on an unseen
subset of the same data set. Fig. 6(b,d) shows a similar performance in
predicting $g_{p}$ as when considering only $\mathbf{E}(t)$ (cf. Fig. 5), also
indicated by the median percentage difference of 4.4 % for decision trees and
3.8 % for CNN. This indicates that the models can extract relevant features
within $E(t)$ no matter the underlying leaf material. In case III, we train
the models on $\mathbf{E}(t)$ and test them on the unseen data set
$\mathbf{E^{\prime}}(t)$. This would correspond to the practical situation
where the method has learned based on water patterns on leaves which
physically are not the same as those on which the model is inferred. The
performance of the methods is rather lousy (Fig. 6e,g), as indicated by the
median percentage difference of 39 % for decision trees and 50 % for CNN. Fig.
6f,h displays the absolute deviation $\Delta=g_{b}-g_{p}$ grouped according to
the five measurement series of $\mathbf{E^{\prime}}(t)$ as a function of the
number of acquisition $i$. Each series starts at low $i$ with $g=0$ after
which both $i$ and $g$ increase concomitantly. It turns out that $\Delta$ is
small at low $g$ and deviates with increasing $g$. Signatures of water
patterns with small droplets are thus overall well recognized by the models,
despite the unseen patterns of $E(t)$ due to the flipped leaf. We conjecture
that this sensitivity may be related to the presence of reminiscent features
of $E(t)$ for $g<7$ mg, due to internal reflections inside the droplets and
leaf material (see Results), on which the network can train. These three test
cases indicate that the learning model needs input training data that are
quite close to the data set from which will be inferred. Although it is
promising that $g$ is well predicted from data corresponding to water patterns
that have not been seen before, unseen variations of the leaf material cause
changes in $E(t)$ to which the methods are not robust. This could be resolved
by training the models on data that besides a large variation of $g$ and water
patterns also include a large diversity of leaves.

The confidence and related error bars of each machine learning model also
depend on specific parameters that are chosen within the architecture and
implementation. For decision trees, for instance, we made use of domain
knowledge to select four temporal regions within $E(t)$. However, a different
choice of the temporal intervals will alter the optimal polynomial order $n$,
modify the feature vector, and thus the performance of the method. For the
convolutional neural network, the specific network architecture expressed by
variables such as the number of layers, filters, and the number of epochs
significantly varies the performance of the network. For the current data set,
the quoted performance due to the described variations has an estimated error
bar of the order of around $1.5\%$. However, also the data set size is an
important parameter. Whereas for decision trees 9000 time traces with 11
features each gives a significant data set of $10^{5}$ features to train on,
for CNNs the final feature size is more difficult to estimate. Although all
temporal points of $\mathbf{E}(t)$ with size (9000, 760) are used as input,
feature extraction leads to much fewer features than the number of elements of
this matrix. By visual inspection of the activations (see Fig. 4a-l), we
estimate that each $E(t)$ provides some $10^{2}$ patterns reminiscent of $g$,
summed over all feature maps. On the other hand, for each $E(t)$ many patterns
together are needed for an accurate prediction of $g$, which reduces the
feature space from the number estimated before. To clarify this point we draw
the analogy with having a picture of a cat, where a cat is defined by its
eyes, ears, tail, fur, etc. Having only one of these features will not lead to
an accurate prediction of the picture showing a cat. We therefore estimate
that the total data size on which the CNN trains should be also of the order
of $10^{5}$, which is on the low side of what is common for CNNs.

The performance of the studied method is also determined by the experimental
setup. The nebulized water does not always entirely end up in the area covered
by the THz beam, thereby contributing to the measured gravimetric weight that
is not seen by the THz beam. We expect that these errors are the cause of the
horizontal lines in Fig. 2d. Additionally, the accuracy of the gravimetric
measurements of $g_{b}$ is rather low, mostly due to natural air convection
and air streams caused by the nebulizers. However, the leaf cannot be placed
in an enclosed box, as then the relative humidity gets spatial and temporal
fluctuations with drastic consequences on $E(t)$. We estimate the absolute
error on $g_{b}$ between 0.1 and 0.2 mg.

Although these results indicate that our data-driven approaches provide a
performant model to determine leaf wetness in a lab setting, the real baptism
of fire for the methodology should be a test in the application environment. A
leaf wetness sensor for the agriculture sector would need to continuously
sense leaf wetness during the growth season on a representative size of the
crop canopy. A reflection geometry would be preferred over a transmission one,
and instead of a single leaf, many leaves will need to be probed
simultaneously. In addition, leaves will not be clean and flat but occur with
a large distribution of appearances. Moreover, the air can be highly humid or
contain rain, and dust, and the beam path may be (temporally) disturbed by
passing insects. We estimate that despite the reported performance of the
algorithms presented in this work, this real setting, which may vary depending
on the canopy, will be challenging for the described methods.

## Conclusions

We have studied the application of conventional machine learning methods to
THz time-domain data to determine leaf wetness. Hereto, we experimentally
acquired a large data set of 12,000 distinct time domain traces $E(t)$
corresponding to a distribution of water patterns on plant leaves. Using
domain knowledge related to the light-matter interaction, we designed and
trained a decision tree model and a convolutional neural network and gave
special attention to feature extraction. Both models predict leaf wetness with
an accuracy of about 4 %. The generalizability of the methods was evaluated on
unseen datasets with increasing deviations from the training set. We observe a
similarity between features important for CNN and those having a physical
interpretation. In conclusion, conventional machine learning models can be of
additional value when compared to model-based signal processing, especially in
cases where the sample configuration is complex and well beyond a multilayer
structure, although the variability with slightly different input data is an
obstacle for practical application.

## Methods

### Samples

The plant material used in this work is Alliaria petiolata, also known as
garlic mustard. The plant is widely spread in the Netherlands, in the wild,
and in private gardens. The plant has been grown in the open ground, and
leaves are harvested in May to June. Of the selected leaf, a circle of 30 mm
diameter is cut and immediately embedded in between two plastic sheets each
0.08 mm thick to prevent it from fast drying out. The plastic sheets have not
been sealed at their edges. During the measurement campaign, the leaf
naturally dried out on the order of $18\%$ (that is, about 11 mg) in 5 days.
This range of relative leaf water content corresponds to variations due to
well-watered vs. severe-drought growth conditions.

To create moisture patterns on the plasticized leaf sample, two nebulizers
(Medisana, Inhalator IN 500) are used. Two specially designed nozzles, which
are symmetrically positioned perpendicular to the plane of incidence and
facing the sample, ensure a homogeneous droplet pattern on the sample. The
nebulizers are filled with distilled water.

### Experimental setup

The optical properties of a moistened plant leaf have been obtained by
performing THz time-domain transmission spectroscopy (Toptica Teraflash Pro).
The sample has been put out-of-focus such that the THz beam covers the surface
area on which water is nebulized. The two nebulizers have been switched on at
random to create a wide range of different droplet patterns. The gravimetric
weight $g_{b}$ of the moist sample is recorded simultaneously with the THz
data. Hereto, the plasticized leaf sample was mounted on a dedicated holder
attached to a precision scale (Sartorius WZA224-L). In addition, an RGB
picture is taken from each water pattern using an optical camera (Logitech
Brio 4K). Simultaneously, also the absolute air humidity $a$ has been
recorded. During the entire measurement campaign, which lasted 7 days, the
temperature was 25.0 ∘C $\pm$ 0.3 ∘C, and the relative humidity was $42\%\pm
6\%$.

Training and inference of the described models are performed on a CPU (Intel
i7-7700HQ @ 2.8 GHz), using 4/4 cores and multi-threading enabled.

### Data acquisition

Using the experimental setup described above, the transmitted electric field
$E(t)$ has experimentally been determined of a leaf sample that has been
moistened until a first droplet ran down the sample. A measurement sequence is
divided into separate time slots: 0.5 s for nebulization, 1.0 s waiting time
where the nebulization cloud deposits itself onto the sample and diffuses away
from the beam path, and 1.5 s for acquiring 50 averages of $E(t)$, recording
$g_{b}$ and $a$, and obtaining an RGB picture of the moisture pattern. For
data set $\mathbf{E}(t)$ the THz beam shines onto the upper side of the leaf,
whereas for data set $\mathbf{E^{\prime}}(t)$ the THz beam shines onto the
bottom side of the leaf. In both cases, the water pattern is directed towards
the emitter side. $\mathbf{E}(t)$ contains a total of 10,609 different
moisture patterns measured in 39 measurement series within 5 days.
$\mathbf{E^{\prime}}(t)$ contains 1,501 patterns recorded in 5 measurement
series within 2 days. Each $E(t)$ contains 760 data points and spans a delay
of 38 ps.

### Data processing

#### Decision trees

Feature extraction is performed by fitting a polynomial function of order $n$
to a selected temporal range $m$ of $E(t)$. As $E(t)$ has a high temporal
variation, and to keep $n$ low for optimal performance of the method, we
choose several ranges to capture all relevant patterns of $E(t)$. The standard
deviation $\sigma$ of $\mathbf{E}(t)$ is enhanced in the range (4, 22) ps (see
Fig. 2b) which we split into $m=4$ equal temporal regions (4, 7); (7, 10);
(15.5, 18.5) and (18.5, 21.5) ps. As $E(t)$ for $g=0$, $E_{0}(t)$, manifests a
slight temporal shift due to factors such as leaf water content, air humidity,
and drift, the onset of all measured traces $E(t)$ occurs at a slightly
different value of $t$. To correct for this, the onset of region $m=1$ is
determined by the slope $d/dt$ of $E(t)/\max(E(t))$. Hereto, $E(t)$ has been
interpolated to enhance the temporal precision. Determination of the order $n$
of each range has been done using a grid search algorithm, which is a brute-
force search of a predefined set of options. Hereto, we calculate the L2 loss
of the validation set after training on $n\in[0,20]$. For the four different
temporal regions of interest we obtain optimal polynomial fits with order
$n_{1}=11,n_{2}=2,n_{3}=4$, and $n_{4}=8$, respectively. Subsequently, we
eliminate recursive features from the input features using permutation
importance as a ranking criterion. By shuffling a feature column in the input
matrix, the correlation between the feature and target value is evaluated. As
such, features with little correlation to the target value are iteratively
eliminated, as well as features that are correlated to others that already
show dependence on the target value. This reduces the number of input features
and in turn, increases the predictive power of the validation set. In this
way, the feature vectors obtained by applying the decision tree algorithm
simultaneously to all regions $m$ read [$t^{2}$, $t^{4}$, $t_{\text{start}}$,
$a$], [bias, $t^{2}$], [$t^{1}$, $t^{2}$], [bias, $t^{2}$, $t^{6}$], for
$m=1..4$ respectively. One can see that although both $a$ and
$t_{\text{start}}$ are present within the features for $m=1$, they are
eliminated from the features for $m>1$ due to mutual correlations.

To train a generalizable model, we use an ensemble of randomized decision
trees, opposite to for instance a single tree which is prone to overfitting.
The randomization is effectuated by the well-established principle of
bootstrap aggregation (or bagging), where different subsets $D_{i}$ are
sampled from training set $D$ with replacement, meaning that the same sample
can be sampled multiple times. We further optimize the trees by tuning
hyperparameters that control regularization, including the number of samples
per tree n_samples, the number of features per tree n_features, and the
maximum depth of branches of a tree max_depth. After optimization of these
hyperparameters using the mentioned grid search algorithm, we obtain n_samples
= 2000, and max_features turns out to be equal to n_features such that the
decision tree becomes a bagged decision tree. The minimum number of samples
after a split is set to 5 to smooth the predicted value of a tree. [26]

### Convolutional Neural network

The architecture of the CNN is shown in Fig. 4m. The input data is inserted in
a concatenation of six groups each containing a convolutional layer. This
structure is motivated by a receptive field which increases with each deeper
layer. Since deeper layers learn more complex patterns, the number of filters
is increased by a factor of two for each layer. Moreover, the convolutional
layers are (after batch normalization, see below) succeeded by a max-pooling
layer that decreases the length by a factor of two. Given the kernel dimension
$(1\times 3)$ and that of the input vector $E(t)$ $(1\times 760)$, we have
chosen for zero padding the extremes of $E(t)$ with a single entry to have the
dimension unaltered after convolution. After this feature extraction part, the
output is flattened and feature $\mathbf{a}$ is appended before inserting it
in a two-layer fully connected network. The entire network has 72,385
trainable parameters: 26,144 are responsible for feature detection and 46,241
for regression.

Input data needs to be normalized before feeding it into the network to ensure
proper convergence during gradient descent, as steepest gradient descent
algorithms do not possess the property of scale invariance. Normalization of
the absolute humidity is performed as
$a_{i}^{\prime}=(a_{i}-\mu_{a})/\sigma_{a}$ where $\mu_{a}$ is the mean of $a$
with respect to the training population and $\sigma_{a}$ its standard
deviation. For all $E(t)$, a uniform division factor of 4 is applied to ensure
all data points are within the range $(-1,1)$. For activation, we make use of
the Rectified Linear Unit (ReLU) function, defined as $\varphi(h)=\max(0,h)$.

Convolutional networks require supervised learning to learn the intrinsic
patterns to predict the target value. Supervised deep learning models use
backpropagation to autonomously shape function $f$ to find relation
$f:x\rightarrow y$, where $x$ is the input and $y$ the target. We have used
the Adam optimizer to find the optimal model parameters. Adam is a simple and
computationally efficient algorithm that combines the advantages of AdaGrad
and RMSProp, resulting in an optimizer that is robust and well-suited to a
wide range of non-convex optimization problems in the field of machine
learning.[27]

## References

  * [1] Koch, M., Mittleman, D. M., Ornik, J. & Castro-Camus, E.  Terahertz time-domain spectroscopy.  _Nature Reviews Methods Primers_ 3, 48, DOI: 10.1038/s43586-023-00232-z (2023).
  * [2] van Mechelen, D.  An industrial THz killer application?  _Optics & Photonics News_ 26, 16–18 (2015).
  * [3] Leitenstorfer, A. _et al._ The 2023 terahertz science and technology roadmap.  _Journal of Physics D: Applied Physics_ 56, 223001, DOI: 10.1088/1361-6463/acbe4c (2023).
  * [4] van Mechelen, J. L. M., Kuzmenko, A. B. & Merbold, H.  Stratified dispersive model for material characterization using terahertz time-domain spectroscopy.  _Optics Letters_ 39, 3853–3856 (2014).
  * [5] S, K., M, Y., Rawson, A. & C. K, S.  Recent Advances in Terahertz Time-Domain Spectroscopy and Imaging Techniques for Automation in Agriculture and Food Sector.  _Food Analytical Methods_ 15, 498–526, DOI: 10.1007/s12161-021-02132-y (2022).
  * [6] Oerke, E.-C.  Crop losses to pests.  _The Journal of Agricultural Science_ 144, 31–43, DOI: 10.1017/S0021859605005708 (2006).
  * [7] Bregaglio, S., Donatelli, M., Confalonieri, R., Acutis, M. & Orlandini, S.  Multi metric evaluation of leaf wetness models for large-area application of plant disease models.  _Agricultural and Forest Meteorology_ 151, 1163–1172, DOI: 10.1016/j.agrformet.2011.04.003 (2011).
  * [8] Huber, L. & Gillespie, T. J.  Modeling Leaf Wetness in Relation to Plant Disease Epidemiology.  _Annual Review of Phytopathology_ 30, 553–577, DOI: 10.1146/annurev.py.30.090192.003005 (1992).
  * [9] Goffart, J.-P. _et al._ Potato production in northwestern europe (germany, france, the netherlands, united kingdom, belgium): Characteristics, issues, challenges and opportunities.  _Potato Research_ 65, 503–547, DOI: 10.1007/s11540-021-09535-8 (2022).
  * [10] Li, R., Lu, Y., Peters, J. M. R., Choat, B. & Lee, A. J.  Non-invasive measurement of leaf water content and pressure–volume curves using terahertz radiation.  _Scientific Reports_ 10, 21028, DOI: 10.1038/s41598-020-78154-z (2020).
  * [11] Gente, R. _et al._ Determination of Leaf Water Content from Terahertz Time-Domain Spectroscopic Data.  _Journal of Infrared, Millimeter, and Terahertz Waves_ 34, 316–323, DOI: 10.1007/s10762-013-9972-8 (2013).
  * [12] Gente, R. & Koch, M.  Monitoring leaf water content with THz and sub-THz waves.  _Plant Methods_ 11, 15, DOI: 10.1186/s13007-015-0057-7 (2015).
  * [13] Singh, A. K., Pérez-López, A. V., Simpson, J. & Castro-Camus, E.  Three-dimensional water mapping of succulent Agave victoriae-reginae leaves by terahertz imaging.  _Scientific Reports_ 10, 1404, DOI: 10.1038/s41598-020-58277-z (2020).
  * [14] Rowlandson, T. _et al._ Reconsidering Leaf Wetness Duration Determination for Plant Disease Management.  _Plant Disease_ 99, 310–319, DOI: 10.1094/PDIS-05-14-0529-FE (2015).  Publisher: Scientific Societies.
  * [15] van Mechelen, J. L. M., Frank, A. & Maas, D. J. H. C.  Thickness sensor for drying paints using THz spectroscopy.  _Optics Express_ 29, 7514, DOI: 10.1364/OE.418809 (2021).
  * [16] Park, H. & Son, J.-H.  Machine Learning Techniques for THz Imaging and Time-Domain Spectroscopy.  _Sensors_ 21, 1186, DOI: 10.3390/s21041186 (2021).
  * [17] Wang, Y. _et al._ Terahertz spectroscopic diagnosis of early blast-induced traumatic brain injury in rats.  _Biomedical Optics Express_ 11, 4085, DOI: 10.1364/BOE.395432 (2020).
  * [18] Cao, C., Zhang, Z., Zhao, X. & Zhang, T.  Terahertz spectroscopy and machine learning algorithm for non-destructive evaluation of protein conformation.  _Optical and Quantum Electronics_ 52, 225, DOI: 10.1007/s11082-020-02345-1 (2020).
  * [19] Wang, Q. _et al._ Automatic defect prediction in glass fiber reinforced polymer based on THz-TDS signal analysis with neural networks.  _Infrared Physics & Technology_ 115, 103673, DOI: 10.1016/j.infrared.2021.103673 (2021).
  * [20] Li, R. _et al._ Nondestructive Evaluation of Thermal Barrier Coatings Thickness Using Terahertz Time-Domain Spectroscopy Combined with Hybrid Machine Learning Approaches.  _MDPI_ 12, 1875, DOI: 10.3390/coatings12121875 (2022).
  * [21] Mao, Q. _et al._ Convolutional neural network model based on terahertz imaging for integrated circuit defect detections.  _Optics Express_ 28, 5000, DOI: 10.1364/OE.384146 (2020).
  * [22] Wang, C. _et al._ Convolutional Neural Network-Based Terahertz Spectral Classification of Liquid Contraband for Security Inspection.  _IEEE Sensors Journal_ 21, 18955–18963, DOI: 10.1109/JSEN.2021.3086478 (2021).
  * [23] Koumans, M., Perez-Casanova, A. & Van Mechelen, J. L. M.  Sensing moisture patterns using terahertz spectroscopy.  In _2022 47th International Conference on Infrared, Millimeter and Terahertz Waves (IRMMW-THz)_ , 1–2, DOI: 10.1109/IRMMW-THz50927.2022.9895781 (IEEE, Delft, Netherlands, 2022).
  * [24] Ioffe, S. & Szegedy, C.  Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift (2015).  ArXiv:1502.03167 [cs].
  * [25] Jördens, C., Scheller, M., Breitenstein, B., Selmar, D. & Koch, M.  Evaluation of leaf water status by means of permittivity at terahertz frequencies.  _Journal of Biological Physics_ 35, 255–264, DOI: 10.1007/s10867-009-9161-0 (2009).
  * [26] Geurts, P., Ernst, D. & Wehenkel, L.  Extremely randomized trees.  _Machine Learning_ 63, 3–42, DOI: 10.1007/s10994-006-6226-1 (2006).
  * [27] Kingma, D. P. & Ba, J.  Adam: A Method for Stochastic Optimization (2017).  ArXiv:1412.6980 [cs].

## Acknowledgements

D.v.M. and B.D. are grateful to fruitful discussions with Niels Anten
(Wageningen University, the Netherlands), and for the facilitating role of
research program Synergia - SYstem change for New Ecology-based and Resource
efficient Growth with high tech In Agriculture, financed by NWO, the Dutch
Research Council (project number 17626), industrial and scientific/research
partners, as well as to Rik Vullings (Eindhoven University of Technology) for
proofreading the manuscript.

## Author contributions statement

D.v.M. designed and supervised the project, B.D. conceived the problem
statement, M.K., D.M, and H.M made the experimental setup and conducted the
experiments, M.K. analysed the results. D.v.M, M.K. and B.D. wrote the
manuscript. All authors reviewed the manuscript.

## Additional information

The authors declare no competing interests.


11institutetext: Department of Applied Informatics, University of Macedonia,
Thessaloniki, Greece 11email: {etsoukanara<EMAIL_ADDRESS>
22institutetext: Computer Science & Engineering, University of Ioannina,
Ioannina, Greece 22email<EMAIL_ADDRESS>
# Should I Stay or Should I Go: Predicting Changes in Cluster Membership

Evangelia Tsoukanara 11 Georgia Koloniari 11 Evaggelia Pitoura 22

###### Abstract

Most research on predicting community evolution focuses on changes in the
states of communities. Instead, we focus on individual nodes and define the
novel problem of predicting whether a specific node stays in the same cluster,
moves to another cluster or drops out of the network. We explore variations of
the problem and propose appropriate classification features based on local and
global node measures. Motivated by the prevalence of machine learning
approaches based on embeddings, we also introduce efficiently computed
distance-based features using appropriate node embeddings. In addition, we
consider chains of features to capture the history of the nodes. Our
experimental results depict the complexity of the different formulations of
the problem and the suitability of the selected features and chain lengths.

###### Keywords:

cluster evolution embeddings feature selection classification

##  1 Introduction

For the problem of community evolution, communities are monitored across time
and their properties studied to attain useful conclusions. Such conclusions
can then be exploited so as to predict community changes through time. Most
works model community evolution through a predefined set of events, such as
community growth or shrinkage, community merging or splitting and so on [5, 7,
12, 16]. The problem is then modeled as a classification problem, where given
a community history, the next event in its evolution is predicted.

However, in many applications, it is important not only to predict the
behavior of a community but also of individual community members. Let us
consider customers that are connected via the common products they buy. These
customers can be clustered according to the companies and products they
prefer. Focusing on individuals in such clusters, makes it possible for
companies to identify and reward loyal customers (community members) or take
preemptive measures to change the behavior of the ones that seem less
dedicated.

To this end, we introduce a novel community evolution prediction problem
defined at node level. For each individual node, there are three possible
events: the node may (a) stay in the same cluster, (b) move to a different
cluster, or (c) drop out of the network. We study variations of the problem by
considering combinations and subsets of the possible evolution events, and
evaluate the use of different classification features for solving them.
Firstly, we consider classic features based on popular node measures.
Motivated by the advances of machine learning methods for community detection
that use node embeddings, we also propose features based on distances between
such embeddings. In particular, we deploy the ComE [4] approach that provides
both node and community embeddings. Features of both methods, defined at
cluster and out of cluster or network level, are combined into chains,
modeling the evolution history of the nodes through time [16]. Our
experimental results show that the problems we defined are not trivial, and
that the proposed distance-based embeddings perform almost as well as the
classic ones, while being computed much more efficiently.

The rest of the paper is structured as follows. Section 2 briefly describes
related work. In Section 3, we formulate the novel problem and its variations.
Section 4 presents our classification features and their modeling into chains.
Section 5 includes our experimental results, while Section 6 concludes.

##  2 Related Work

We discuss related research, first, on community evolution and then, on node
and graph embeddings.

Community Evolution. After discovering communities in evolving networks [15],
their properties are studied by mapping corresponding communities through
successive network snapshots [1]. Community evolution is assessed with
measures such as its growth and disappearance rate [20], or its life
expectancy [10].

For predicting community evolution [5, 12, 16], events such as community
growth, shrinkage, merging and splitting are defined. The problem is modeled
as a classification problem in which features based on the structural
properties of communities are exploited, and given the history of a community
the next event in the community’s evolution is predicted. In our work, instead
of communities, we focus on nodes, and introduce a new problem aiming at
predicting changes in the nodes’ cluster membership through time.

Node Embeddings. An embedding is the transformation of a high-dimensional
space to a low-dimension vector. For graphs, the focus is mostly on node
embeddings that preserve network structure. Deepwalk [14] learns node
embeddings that capture second-order proximity (i.e., proximity between shared
neighbors) by simulating short random walks and applying the Skip-gram
algorithm. LINE [17] employs edge-sampling, and to preserve both first-order
(i.e., ties between neighbors) and second-order proximity, it is first trained
separately and then the two resulting embeddings are concatenated. Node2vec
[6] extends Deepwalk by generating biased random walks to explore diverse node
neighborhoods. Finally, GraRep [3] and HOPE [9] derive embeddings that capture
high-order proximity.

Besides, link prediction, node classification and visualization, node and
graph embeddings are also used for community detection [8, 19]. ComE
(Community Embedding) [4] is a framework that jointly solves both community
detection, and learning node and community embeddings. The intuition is that
node embeddings that capture community-aware proximity, can assist community
detection, while community embeddings can in turn improve node embeddings. In
our work, we deploy ComE and explore whether node and community embeddings can
be used to derive predictive features.

##  3 Problem Formulation

A social network is often represented as a graph $G=(V,E)$, where $V$ is the
set of nodes (vertices) and $E$ is the set of edges. A temporal social network
is a network that changes over time and is represented as a sequence of graphs
$\\{G_{1},G_{2},\dots,G_{n}\\}$, where $G_{i}=(V_{i},E_{i})$, represents a
snapshot of graph $G$ at time $i$, and $V_{i}$ and $E_{i}$ are the node and
edge sets at time $i$ respectively. Let
$\mathcal{C}_{i}=\\{C^{1}_{i},C^{2}_{i},\dots,C^{m}_{i}\\}$ be a clustering of
$G_{i}$ consisting of $m$ clusters, such that $C^{j}_{i}\cap
C^{k}_{i}=\emptyset,1\leq j,k\leq m,j\neq k$. For two clusterings
$\mathcal{C}_{i}$ and $\mathcal{C}_{i+1}$ at consecutive timeframes $i$ and
$i+1$, we assume cluster $C^{j}_{i+1}$ is the evolution of cluster
$C^{j}_{i}$. Similarly, for multiple consecutive timeframes,
$C^{j}_{1},C^{j}_{2},\dots,C^{j}_{n}$ denotes the evolution of cluster $C^{j}$
in time period $[1,n]$.

Inspired by the idea of community evolution, we study the problem at node
level by aiming to predict how node memberships in clusters evolve through
time. We discern between different states that a node can have with respect to
its cluster membership in the next timeframe, i.e., a node can stay in the
same cluster, move to another or drop out of the network. Based on the above,
we define our problem as follows.

######  Definition 1 (Stay\Move\Drop ($\mathcal{SMD}$) Problem)

Given a sequence of clusterings $\mathcal{C}_{1},\dots,\mathcal{C}_{i}$,
corresponding to a consecutive set of timeframes for a graph $G$, and node
$v\in C^{j}_{i}$, predict the state of node $v$ regarding its evolution in the
next timeframe $i+1$ as state:

  * •

stay, $\mathcal{S}$: node $v$ stays at the same cluster in $i+1$, that is
$v\in C^{j}_{i+1}$,

  * •

move, $\mathcal{M}$: node $v$ moves to another cluster in $i+1$, that is $v\in
C^{k}_{i+1}$, $k\neq j$, and

  * •

drop, $\mathcal{D}$: node $v$ drops out of the network, that is $v\notin
V_{i+1}$.

Thus, we define a classification problem with 3 classes, labeled stay, move
and drop. Given the history of a node in a given time period, defined as a
sequence of distinct timeframes, we predict its class in the next time frame.

If we are only interested in discerning between loyal cluster members and
members likely to leave, we may simplify our problem to a binary
classification problem. This first alternative problem, is derived by merging
states $move$ and $drop$, in one class $leave$. Thus, the classes are reformed
as follows:
Stay\Leave ($\mathcal{SL}$) Problem: the possible node events are reformed as
state:

  * •

stay, $\mathcal{S}$: node $v$ stays at the same cluster in $i+1$, that is
$v\in C^{j}_{i+1}$, and

  * •

leave, $\mathcal{L}$: node $v$ does not remain in the same cluster, that is
$v\notin C^{j}_{i+1}$.

Finally, we omit the third class of the $\mathcal{SMD}$ problem, and only
provide predictions for nodes that remain in the network in timeframe $i+1$.
For the third variation, we have:
Stay\Move ($\mathcal{SM}$) Problem.

  * •

stay, $\mathcal{S}$: node $v$ stays at the same cluster in $i+1$, that is
$v\in C^{j}_{i+1}$, and

  * •

move, $\mathcal{M}$: node $v$ moves to another cluster in $i+1$, that is $v\in
C^{k}_{i+1}$, $k\neq j$.

(a) $T_{0}$

(b) Syntgen dataset

(c) DBLP dataset

Figure 1: (a) ComE clustering and $\Delta$ distribution for (b) Syntgen and
(c) DBLP.

##  4 Predictive Features & Historical Chains

To solve the classification problems we define, the evaluation and selection
of appropriate predictive features is required. We discern between two basic
types of features, based on structural network measures that are usually
exploited for community evolution predictions, and on distances between
embeddings that we propose. As all three problems we introduce are variations
of the same classification problem, we define and evaluate the same features
for all problems.

Classic Features. Classic features, defined on node level, provide information
about the structural role of a node in the network. In particular, we select:
(a) degree measuring the connections of a node, (b) betweenness, measuring the
number of shortest paths that pass through a node, (c) closeness, that
measures the distance of a node to all other nodes, and (d) eigenvector
centrality measuring the influence of a node in the network, defined on
cluster and network level.

Aggregated at community level, classic features offer insights for predicting
community evolution [12, 16], while for individual nodes, such measures also
contain information about their evolution tendencies. For instance, an
influential central node with high degree is less likely to drop out of the
network compared to a low-degree remote node. Similarly focusing on community
structure, well-connected core nodes within a community are more likely to
stay in their cluster compared to loosely connected border nodes. Thus, we
also differentiate between features defined at cluster level (in) and at
network level (out).

Embeddings-based Features. Since node embeddings are low dimensional vector
representations of nodes that capture structural graph properties, we propose
defining predictive features based on such embeddings. The ComE [4] framework
provides an appropriate solution by solving both community detection and
embeddings learning jointly, focusing on deriving embeddings that capture
community-aware proximity between nodes. ComE uses as input a graph’s edge
list and the number of target clusters $k$ and outputs: (i) for each node its
embedding along with its community membership and (ii) for each community,
which is defined as a multivariate Gaussian distribution, its embedding
parameters, that is, a median vector (i.e., the embedding of the mean of the
community) and a covariance matrix. Fig. 1(a) depicts the node embeddings for
the six clusters detected by ComE for a synthetic dataset generated by the
Syntgen [13] generator for timeframe $T_{0}$, where we notice that ComE
manages to detect reasonable well-separated clusters with similar node
embeddings as projected in 2-d.

We define features on cluster and network level, exploiting the outputs of
ComE. For computing network level features, we exclude the nodes of the
cluster the given node belongs to, to avoid cases of identical network and
cluster features. Let $\phi_{v}$ be the embedding of node $v$, and $\phi_{C}$
the embedding of the median of community $C$. Without loss of generality we
define our features using the euclidean distance ($d$) between pairs of
embeddings. Cosine and $L1$ distances were also evaluated, but euclidean were
used as they performed slightly better. In particular, for a node $v$ such
that $v\in C$, we define the following features:

Cluster level:

  * •

distance from cluster median: $d(\phi_{v},\phi_{C})$

  * •

distance from least similar cluster member: $\max_{\forall u\in
C}d(\phi_{v},\phi_{u})$

  * •

distance from most similar cluster member: $\min_{\forall u\in
C}d(\phi_{v},\phi_{u})$

  * •

avg. distance from all cluster members: $avg_{\forall u\in C}$
$d(\phi_{v},\phi_{u})$

Network level:

  * •

min. distance from other cluster median: $\min_{\forall C^{\prime}\neq
C}d(\phi_{v},\phi_{C^{\prime}})$

  * •

distance from least similar node out of the cluster: $\max_{\forall u\notin
C}d(\phi_{v},\phi_{u})$

  * •

distance from most similar node out of the cluster: $\min_{\forall u\notin
C}d(\phi_{v},\phi_{u})$

  * •

avg. distance from all nodes not in the cluster: $avg_{\forall u\notin C}$
$d(\phi_{v},\phi_{u})$

Our approach is based on the idea that the embeddings of nodes in a cluster
are more similar. Thus, if a node has an embedding that is similar to nodes of
other clusters, it is more likely to change or leave its cluster. In Fig. 1(b)
and Fig. 1(c), we plot the distribution of the difference, $\Delta$, between
$\min_{\forall C^{\prime}\neq C}d(\phi_{v},\phi_{C^{\prime}})$ and
$d(\phi_{v},\phi_{C})$ features, for classes stay and move for a synthetic
dataset generated by Syntgen and a citation DBLP dataset based on [18]. In
both figures, we notice that nodes that move to another cluster tend to have
lower or even negative $\Delta$ compared to the ones that remain in the same
cluster. About 60% of the move nodes have $\Delta$ less than 0.20 for both
datasets, while more than 60% of the stay nodes have $\Delta$ higher than
0.28. Therefore, the values of the various distance measures can provide
insight on the properties and behavior of a node, and thus, are appropriate as
predictive features for our classification problems.

Historical Chains of Features. The evolution of a community is tracked in
successive network snapshots that correspond to successive timeframes. Thus,
all features we describe can be measured for each timeframe. Let us assume a
set of $k$ predictive features for each node $v$, and a time period
$[1,\dots,n]$, where $f^{j}_{i}(v)$ denotes the $j$-th feature of node $v$ at
time $i$. Further, for every pair of consecutive timeframes $i$ and $i+1$ in
the given time period, the state (label), $l_{i+1}(v)$, of node $v$ can be
recorded.

To model the history of the node and exploit it to derive more accurate
predictions, we utilize historical chains of features as defined in [16]. In
particular, we have as final features for $v$:
$\\{f^{1}_{1}(v),\dots,f^{k}_{1}(v)\\}$, $l_{2}(v)$, $\\{f^{1}_{2}(v),$
$\dots,$ $f^{k}_{2}(v)\\}$, $l_{3}(v)$ $\dots$, $\\{f^{1}_{n}(v),$ $\dots,$
$f^{k}_{n}(v)\\}$, while $l_{n+1}(v)$ is the label to be predicted.

Though we have defined here one chain to model the entire node history, the
use of subchains of various lengths can also be deployed. While a longer
history will provide more information regarding the history of a node, it
would limit the number of nodes for which a prediction can be made.

Table 1: Snapshot Structure | DBLP | Email-eu | Syntgen
---|---|---|---
Sn$\\#$ | Nodes | Edges | $|\mathcal{C}|$ | $\mathcal{Q}$ | Nodes | Edges | $|\mathcal{C}|$ | $\mathcal{Q}$ | Nodes | Edges | $|\mathcal{C}|$ | $\mathcal{Q}$
0 | 14731 | 120192 | 17 | 0.681 | 750 | 4740 | 11 | 0.515 | 1583 | 8955 | 6 | 0.525
1 | 16801 | 143404 | 16 | 0.676 | 745 | 5077 | 11 | 0.431 | 1443 | 7936 | 5 | 0.480
2 | 17756 | 156393 | 16 | 0.649 | 742 | 4578 | 10 | 0.528 | 1367 | 7249 | 5 | 0.459
3 | 14765 | 120370 | 17 | 0.659 | 742 | 4886 | 9 | 0.487 | 1567 | 8154 | 5 | 0.476
4 | 10898 | 69005 | 16 | 0.660 | 749 | 5072 | 11 | 0.388 | 1575 | 7993 | 5 | 0.483
5 |  |  |  |  | 739 | 4819 | 10 | 0.521 | 1402 | 6978 | 5 | 0.465
6 |  |  |  |  | 759 | 4846 | 10 | 0.521 | 1526 | 7629 | 5 | 0.472
7 |  |  |  |  | 808 | 5405 | 11 | 0.410 | 1416 | 6973 | 5 | 0.462
8 |  |  |  |  | 772 | 4880 | 13 | 0.391 | 1409 | 6853 | 5 | 0.470
9 |  |  |  |  | 785 | 5169 | 12 | 0.533 | 1594 | 7794 | 6 | 0.490

##  5 Evaluation

We experimentally study all three classification problems while comparing
different sets of predictive features.

Table 2: Performance for the $\mathcal{SM}$ Problem | ComE | Classic
---|---|---
Data | Feat. | Class | P | R | F1 | Acc | P | R | F1 | Acc
DBLP | $in$ | $\mathcal{S}$ | 0.821 | 0.940 | 0.876 | 0.804 | 0.837 | 0.941 | 0.886 | 0.821
$\mathcal{M}$ | 0.715 | 0.425 | 0.533 | 0.746 | 0.485 | 0.587
$out$ | $\mathcal{S}$ | 0.800 | 0.942 | 0.865 | 0.784 | 0.783 | 0.927 | 0.849 | 0.757
$\mathcal{M}$ | 0.679 | 0.339 | 0.452 | 0.577 | 0.279 | 0.376
$all$ | $\mathcal{S}$ | 0.836 | 0.954 | 0.891 | 0.828 | 0.839 | 0.948 | 0.890 | 0.827
$\mathcal{M}$ | 0.787 | 0.476 | 0.592 | 0.769 | 0.490 | 0.599
Email-eu | $in$ | $\mathcal{S}$ | 0.807 | 0.926 | 0.862 | 0.796 | 0.797 | 0.912 | 0.850 | 0.778
$\mathcal{M}$ | 0.757 | 0.507 | 0.606 | 0.712 | 0.483 | 0.574
$out$ | $\mathcal{S}$ | 0.778 | 0.933 | 0.848 | 0.770 | 0.732 | 0.910 | 0.812 | 0.709
$\mathcal{M}$ | 0.733 | 0.409 | 0.524 | 0.569 | 0.261 | 0.357
$all$ | $\mathcal{S}$ | 0.823 | 0.934 | 0.875 | 0.816 | 0.789 | 0.925 | 0.852 | 0.778
$\mathcal{M}$ | 0.792 | 0.554 | 0.652 | 0.731 | 0.452 | 0.558
Syntgen | $in$ | $\mathcal{S}$ | 0.675 | 0.692 | 0.683 | 0.668 | 0.707 | 0.719 | 0.713 | 0.700
$\mathcal{M}$ | 0.661 | 0.643 | 0.652 | 0.693 | 0.680 | 0.687
$out$ | $\mathcal{S}$ | 0.638 | 0.704 | 0.670 | 0.640 | 0.642 | 0.669 | 0.655 | 0.635
$\mathcal{M}$ | 0.643 | 0.572 | 0.606 | 0.628 | 0.599 | 0.613
$all$ | $\mathcal{S}$ | 0.693 | 0.725 | 0.708 | 0.691 | 0.710 | 0.713 | 0.711 | 0.701
$\mathcal{M}$ | 0.690 | 0.656 | 0.672 | 0.691 | 0.687 | 0.689
Table 3: Performance for the $\mathcal{SL}$ Problem | ComE | Classic
---|---|---
Data | Feat. | Class | P | R | F1 | Acc | P | R | F1 | Acc
DBLP | $in$ | $\mathcal{S}$ | 0.626 | 0.676 | 0.650 | 0.922 | 0.715 | 0.687 | 0.700 | 0.937
$\mathcal{L}$ | 0.961 | 0.951 | 0.956 | 0.962 | 0.967 | 0.965
$out$ | $\mathcal{S}$ | 0.632 | 0.664 | 0.647 | 0.922 | 0.680 | 0.648 | 0.664 | 0.929
$\mathcal{L}$ | 0.959 | 0.953 | 0.956 | 0.958 | 0.963 | 0.960
$all$ | $\mathcal{S}$ | 0.680 | 0.710 | 0.695 | 0.933 | 0.728 | 0.690 | 0.709 | 0.939
$\mathcal{L}$ | 0.965 | 0.960 | 0.962 | 0.963 | 0.969 | 0.966
Email-eu | $in$ | $\mathcal{S}$ | 0.796 | 0.923 | 0.855 | 0.820 | 0.799 | 0.892 | 0.843 | 0.809
$\mathcal{L}$ | 0.869 | 0.680 | 0.762 | 0.827 | 0.698 | 0.756
$out$ | $\mathcal{S}$ | 0.771 | 0.915 | 0.837 | 0.795 | 0.725 | 0.882 | 0.796 | 0.740
$\mathcal{L}$ | 0.847 | 0.633 | 0.724 | 0.775 | 0.548 | 0.642
$all$ | $\mathcal{S}$ | 0.809 | 0.934 | 0.867 | 0.836 | 0.793 | 0.892 | 0.840 | 0.804
$\mathcal{L}$ | 0.889 | 0.703 | 0.784 | 0.826 | 0.686 | 0.749
Syntgen | $in$ | $\mathcal{S}$ | 0.656 | 0.654 | 0.655 | 0.708 | 0.696 | 0.690 | 0.693 | 0.741
$\mathcal{L}$ | 0.747 | 0.748 | 0.747 | 0.774 | 0.778 | 0.776
$out$ | $\mathcal{S}$ | 0.634 | 0.660 | 0.646 | 0.694 | 0.631 | 0.625 | 0.628 | 0.686
$\mathcal{L}$ | 0.742 | 0.720 | 0.731 | 0.726 | 0.731 | 0.729
$all$ | $\mathcal{S}$ | 0.670 | 0.682 | 0.676 | 0.723 | 0.693 | 0.684 | 0.688 | 0.738
$\mathcal{L}$ | 0.763 | 0.753 | 0.758 | 0.770 | 0.778 | 0.774

###  5.1 Datasets

We use three datasets for our evaluation, two real and one synthetic.

DBLP: DBLP is a citation network111https://www.aminer.org/citation that
includes additional information about the publications, such as year of
publication and fields of study they belong to [18]. Similarly to [4], we
filter papers with primary of study as NLP, Databases, Networking, Data Mining
and Computer Vision. We construct five snapshots, for years 2015 to 2019, by
maintaining publications of the given year and adding cited papers that belong
to the selected fields regardless of their publication year. To attain a
denser network, we sample nodes with degree $\geq$ 20 and build the induced
undirected subgraph.

Email-eu: The Email-eu222https://snap.stanford.edu/data/email-Eu-core-
temporal.html [11] dataset consists of incoming and outgoing e-mails between
members of a large European institution in a period of 803 days. We consider
the network undirected and split the data into 10 balanced snapshots.

Syntgen: To further investigate the impact of network properties on our
problems, we use the Syntgen generator [13] that creates temporal undirected
networks simulating real networks using explicit specifications, like degree
distributions and cluster sizes, as well as implicitly controlling the
perseverance of nodes popularity over time. The intra-cluster to total degree
ratio determines cluster density. A high ratio leads to dense well-separated
communities, while low values exhibit no clustering. We set the default ratio
to 0.7.

We apply ComE [4] to detect communities at each snapshot. To determine an
appropriate number of clusters $k$ as input for ComE, for DBLP and Email-eu,
we first apply the Louvain [2] community detection method that selects the $k$
that maximizes modularity. As to Syntgen, we use as $k$ the number of clusters
obtained from the generator. Clusters are mapped across different snapshots
based on the majority of their common nodes. Table 1 presents the number of
nodes and edges as well as the number of clusters ($|\mathcal{C}|$) and
modularity ($Q$) for each network snapshot for all datasets.

With regards to historical chains, DBLP with 5 snapshots can form 2-length up
to 4-length chains, while Email-eu and Syntgen with 10 snapshots can form from
2-length up to 9-length chains. Trying to balance between more information
that longer chains provide and the ability to provide predictions for more
nodes, we use as default chain length 5 for both Email-eu and Syntgen, and 2
for DBLP.

(a) $\mathcal{SM}$ Problem

(b) $\mathcal{SL}$ Problem

(c) $\mathcal{SMD}$ Problem

Figure 2: Macro f1 score per class.

###  5.2 Experimental Results

We report results using the Random Forest classifier, which performed better
compared to other classifiers we tried (e.g. Logistic Regression, Naive
Bayes). We apply stratified 5-fold cross-validation to preserve the same class
distribution in both train and test sets. Tables 2 to 4 illustrate precision,
recall and f1-score per class, and accuracy for each problem when using ComE
and Classic features at cluster ($in$) and network level ($out$), and their
combination ($all$) for each dataset. We denote with bold the best f1-score
for each class per dataset and problem. As structural features are typically
used in community evolution prediction, we consider Classic features as a
baseline method to compare the proposed ComE-based features.

General observations. A first observation derived by Table 2 to Table 4 is
that while the three classification problems are not trivial, both Classic and
ComE based features show promising initial results, achieving high accuracy
but showcasing that some classes are more difficult to predict than others. We
also notice that $all$ features perform better in most cases, and use them as
our default features for the rest of our study.

Table 4: Performance for the $\mathcal{SMD}$ Problem | ComE | Classic
---|---|---
Data | Feat. | Class | P | R | F1 | Acc | P | R | F1 | Acc
DBLP | $in$ | $\mathcal{S}$ | 0.589 | 0.818 | 0.685 | 0.905 | 0.666 | 0.778 | 0.718 | 0.917
$\mathcal{M}$ | 0.538 | 0.358 | 0.430 | 0.587 | 0.315 | 0.410
$\mathcal{D}$ | 0.973 | 0.941 | 0.957 | 0.962 | 0.962 | 0.962
$out$ | $\mathcal{S}$ | 0.587 | 0.821 | 0.685 | 0.904 | 0.627 | 0.744 | 0.680 | 0.911
$\mathcal{M}$ | 0.491 | 0.250 | 0.331 | 0.463 | 0.161 | 0.239
$\mathcal{D}$ | 0.972 | 0.944 | 0.958 | 0.961 | 0.966 | 0.963
$all$ | $\mathcal{S}$ | 0.621 | 0.853 | 0.718 | 0.914 | 0.682 | 0.778 | 0.727 | 0.921
$\mathcal{M}$ | 0.601 | 0.384 | 0.468 | 0.616 | 0.310 | 0.411
$\mathcal{D}$ | 0.976 | 0.945 | 0.960 | 0.963 | 0.967 | 0.965
Email-eu | $in$ | $\mathcal{S}$ | 0.795 | 0.932 | 0.858 | 0.816 | 0.783 | 0.909 | 0.841 | 0.790
$\mathcal{M}$ | 0.748 | 0.517 | 0.611 | 0.653 | 0.468 | 0.545
$\mathcal{D}$ | 1 | 0.881 | 0.936 | 0.996 | 0.881 | 0.935
$out$ | $\mathcal{S}$ | 0.762 | 0.925 | 0.836 | 0.785 | 0.718 | 0.897 | 0.798 | 0.732
$\mathcal{M}$ | 0.686 | 0.413 | 0.515 | 0.524 | 0.272 | 0.357
$\mathcal{D}$ | 1 | 0.881 | 0.936 | 0.993 | 0.881 | 0.933
$all$ | $\mathcal{S}$ | 0.797 | 0.939 | 0.862 | 0.821 | 0.775 | 0.905 | 0.835 | 0.784
$\mathcal{M}$ | 0.765 | 0.520 | 0.619 | 0.647 | 0.451 | 0.530
$\mathcal{D}$ | 1 | 0.881 | 0.937 | 0.990 | 0.881 | 0.932
Syntgen | $in$ | $\mathcal{S}$ | 0.645 | 0.680 | 0.662 | 0.693 | 0.683 | 0.707 | 0.695 | 0.724
$\mathcal{M}$ | 0.634 | 0.646 | 0.640 | 0.668 | 0.696 | 0.681
$\mathcal{D}$ | 1 | 0.828 | 0.906 | 1 | 0.828 | 0.905
$out$ | $\mathcal{S}$ | 0.620 | 0.692 | 0.654 | 0.680 | 0.620 | 0.658 | 0.638 | 0.670
$\mathcal{M}$ | 0.628 | 0.599 | 0.613 | 0.603 | 0.610 | 0.606
$\mathcal{D}$ | 1 | 0.828 | 0.906 | 1 | 0.828 | 0.906
$all$ | $\mathcal{S}$ | 0.661 | 0.712 | 0.685 | 0.711 | 0.685 | 0.706 | 0.695 | 0.725
$\mathcal{M}$ | 0.659 | 0.657 | 0.658 | 0.668 | 0.698 | 0.682
$\mathcal{D}$ | 1 | 0.828 | 0.906 | 1 | 0.828 | 0.906

$\boldsymbol{\mathcal{SM}}$ vs. $\boldsymbol{\mathcal{SL}}$ vs.
$\boldsymbol{\mathcal{SMD}}$. To compare the three problems, we illustrate the
f1-score for each class for both types of features (ComE, denoted as CM, and
Classic denoted as CL) for the three datasets in Fig. 2. For the
$\mathcal{SM}$ problem, Fig. 2(a) shows that class stay performs better than
move for DBLP and Email-eu. DBLP achieves the best performance with 0.891 for
stay and 0.592 for move respectively for the ComE features and similar results
for the Classic ones (Table 2). This is due to the imbalance in the real
datasets between the two classes, with the majority of the nodes in class
stay. In contrast, in the Syntgen dataset, classes are well-balanced and we
notice similar performance for both. For the $\mathcal{SL}$ problem (Fig.
2(b)), we notice significant difference mainly on DBLP. In this case, class
stay is underrepresented due to the network construction. Looking at Table 3,
class stay achieves f1 0.695 and leave 0.962, for DBLP with feature type
$all$. Finally, as we can see in Table 4 and Fig. 2(c) for the $\mathcal{SMD}$
problem, class move is heavily underrepresented in both DBLP and Email-eu
resulting in a rather low f1-score. Summing up, the $\mathcal{SL}$ problem
seems to have the best overall performance, while $\mathcal{SM}$ appears to
have the worst. Apparently, class move is overall the most difficult to
predict. Besides class imbalance that makes the problem more difficult, this
occurs especially when the characteristics across communities are not
significantly different, which is a similarity indicator between communities.

(a) $\mathcal{SM}$ Problem

(b) $\mathcal{SL}$ Problem

(c) $\mathcal{SMD}$ Problem

Figure 3: Macro f1 score per feature category.

Selecting appropriate features. Next, we focus on the comparison between the
different types of features $in$ (I), $out$ (O) and $all$ (A) for ComE (CM)
and Classic (CL) features and how they perform at each problem. As we can see
in Fig. 3, $in$ and $all$ outperform $out$ features. For ComE features $all$
is the best choice for all problems, while for Classic features $in$ performs
sometimes better. This seems to depend on the dataset and not the problem we
study, as we see that for Email-eu Classic $in$ features perform better for
all problems (with the exception of class stay on the $\mathcal{SM}$ problem).
For the $\mathcal{SM}$ problem and ComE features, we notice a significant
difference between $out$ and $all$ features (Fig. 3(a)). Macro f1 is 0.613 for
$out$ and 0.744 for $all$ for the DBLP dataset and 0.584 for $out$ and 0.704
for $all$ for the Email-eu dataset. Overall, we do not observe significant
differences between ComE and Classic features, deducing that the ComE features
offer satisfying performance while being more efficiently computed. For
instance, for the DBLP dataset, the computation time is 5493s for ComE
features, while Classic features are slower by an order, requiring 34618s.

Influence of the length of historical chains. Depicted in Fig. 4(a), 4(b), and
4(c), we explore the effect of chains of features with varying length for the
DBLP (D), Email-eu (E) and Syntgen (S). Both Email-eu and Syntgen show that f1
generally improves as chain length grows for all problems. The Syntgen dataset
follows the same pattern for the $\mathcal{SL}$ and $\mathcal{SMD}$ problems,
increasing and reaching its peak at chain length 8. For the Email-eu, we
observe a temporary drop at length 7, while the highest score is reached at
length 9, with 0.804 for the $\mathcal{SM}$ (Fig. 4(a)), 0.883 for the
$\mathcal{SL}$ (Fig. 4(b)) and 0.850 for the $\mathcal{SMD}$ (Fig. 4(c))
problem respectively. As we have mentioned, while longer chains provide more
information and more accurate predictions, they are not available for a large
number of nodes. In particular, Email-eu consists of 6034 instances at
2-length chain and only 750 instances at its longest chain. Similarly, Syntgen
consists of 11879 instances at 2-length and 1583 instances at 9-length chain.
With regards to DBLP, the best score, 0.828, is achieved at the $\mathcal{SL}$
problem with chain length 2 with 0.828, but the history is too limited to
derive safe comparative conclusions.

Influence of intra-cluster to total degree ratio. In the last experiment, we
focus on the Syntgen generator producing different datasets with varying
intra-cluster to total degree ratio, which determines the density within the
constructed clusters compared to the overall network. In Fig. 4(d), we notice
a sharp increase on macro f1 for ratio 0.6 up to 0.8 for all problems. Lower
ratio indicates poor clustering and thus is not suitable for our context.
Beyond 0.8, behavior diverges. In such tightly connected communities, most
nodes have similar roles in their community, making it difficult for a
classifier to determine their behavior. As a conclusion, cases with ratio
close to 0.5 that exhibit no locality, or close to 0.9 indicating almost
disconnected communities, fail to achieve good results.

(a) $\mathcal{SM}$ Problem

(b) $\mathcal{SL}$ Problem

(c) $\mathcal{SMD}$ Problem

(d) Increasing ratio

Figure 4: Macro f1 score (a), (b), (c) per chain length and (d) per degree
ratio.

##  6 Conclusions

In this paper, we defined a novel problem, related to community evolution,
that focuses on nodes and aims at predicting whether they will stay in their
cluster, move to another or drop out of the network. We modeled the problem as
a classification problem and with three variations. We determined appropriate
features, based on both local and global node measures, and formed chains of
features to take advantage of node history. We also proposed exploiting node
and community embeddings derived by the ComE [4] framework to define distance
based features. Our experimental results showed that the novel problem is not
trivial, and the distance-based features performed similarly to the Classic
ones, while requiring far less computation time. Next, we will consider
alternative community learning approaches to derive node embeddings and define
appropriate features.

## References

  * [1] Aynaud, T., Fleury, E., Guillaume, J.L., Wang, Q.: Communities in evolving networks: Definitions, detection, and analysis techniques. In: Dynamics On and Of Complex Networks, Volume 2: Applications to Time-Varying Dynamical Systems, pp. 159–200. Springer New York (2013)
  * [2] Blondel, V.D., Guillaume, J.L., Lambiotte, R., Lefebvre, E.: Fast unfolding of communities in large networks. Journal of Statistical Mechanics: Theory and Experiment 2008(10), P10008 (2008)
  * [3] Cao, S., Lu, W., Xu, Q.: Grarep: Learning graph representations with global structural information. In: Proc. of the 24th ACM CIKM. p. 891–900 (2015)
  * [4] Cavallari, S., Zheng, V.W., Cai, H., Chang, K.C.C., Cambria, E.: Learning community embedding with community detection and node embedding on graphs. In: Proc. of the 2017 ACM CIKM. p. 377–386 (2017)
  * [5] Gliwa, B., Bródka, P., Zygmunt, A., Saganowski, S., Kazienko, P., Koźlak, J.: Different approaches to community evolution prediction in blogosphere. In: Proc. of the 2013 IEEE/ACM ASONAM. p. 1291–1298 (2013)
  * [6] Grover, A., Leskovec, J.: Node2vec: Scalable feature learning for networks. In: Proc. of the 22nd ACM SIGKDD. p. 855–864 (2016)
  * [7] Ilhan, N., Öğüdücü, Ş.G.: Predicting community evolution based on time series modeling. In: Proc. of the 2015 IEEE/ACM ASONAM. p. 1509–1516 (2015)
  * [8] Kozdoba, M., Mannor, S.: Community detection via measure space embedding. In: Advances in Neural Inf. Proc. Sys. 28: NIPS 2015. pp. 2890–2898 (2015)
  * [9] Ou, M., Cui, P., Pei, J., Zhang, Z., Zhu, W.: Asymmetric transitivity preserving graph embedding. In: Proc. of the 22nd ACM SIGKDD. p. 1105–1114 (2016)
  * [10] Palla, G., Barabási, A.L., Vicsek, T.: Quantifying social group evolution. Nature 446, 664–667 (2007)
  * [11] Paranjape, A., Benson, A.R., Leskovec, J.: Motifs in temporal networks. In: Proc. of the Tenth ACM WSDM. p. 601–610 (2017)
  * [12] Pavlopoulou, M.E.G., Tzortzis, G., Vogiatzis, D., Paliouras, G.: Predicting the evolution of communities in social networks using structural and temporal features. In: 12th Int. Work. on Semantic and Social Media Adaptation and Personalization. pp. 40–45 (2017)
  * [13] Pereira, L.R., Lopes, R.J., Louçã, J.: Syntgen: a system to generate temporal networks with user-specified topology. Journal of Complex Networks 4(0), 1–26 (2019)
  * [14] Perozzi, B., Al-Rfou, R., Skiena, S.: Deepwalk: Online learning of social representations. In: Proc. of the 20th ACM SIGKDD. p. 701–710 (2014)
  * [15] Rossetti, G., Cazabet, R.: Community discovery in dynamic networks: A survey. ACM Comput. Surv. 51(2) (2018)
  * [16] Saganowski, S.: Predicting community evolution in social networks. In: Proc. of the 2015 IEEE/ACM ASONAM. p. 924–925 (2015)
  * [17] Tang, J., Qu, M., Wang, M., Zhang, M., Yan, J., Mei, Q.: Line: Large-scale information network embedding. In: Proc. of the 24th WWW. p. 1067–1077 (2015)
  * [18] Tang, J., Zhang, J., Yao, L., Li, J., Zhang, L., Su, Z.: Arnetminer: Extraction and mining of academic social networks. In: Proc. of the 14th ACM SIGKDD. p. 990–998 (2008)
  * [19] Tian, F., Gao, B., Cui, Q., Chen, E., Liu, T.Y.: Learning deep representations for graph clustering. In: Proc. of the 28th AAAI. p. 1293–1299 (2014)
  * [20] Toyoda, M., Kitsuregawa, M.: Extracting evolution of web communities from a series of web archives. In: Proc. of the 14th ACM Conf. on Hypertext and Hypermedia. p. 28–37 (2003)


# PREGNANCY LOSS AND UNETHICAL ALGORITHMS: ETHICAL ISSUES IN TARGETED
ADVERTISING

Fatemeh Golpayegani1 1Fatemeh is with School of Computer Science, University
College Dublin, Dublin, Ireland<EMAIL_ADDRESS>
##  I EXTENDED ABSTRACT

Pregnancy loss which happens for one in four pregnancies can be a traumatic
experience that impacts women’s personal and social life in several ways.
Studies have shown that women who have experienced pregnancy loss are at a
great risk of depression and anxiety and they are more likely to have major
social difficulties. This can be followed by social life distortion and
isolation leaving social media and disclosure to anonymous online support
groups a way to grief as well as seeking some help. However, so-called
intelligent targeted advertising will follow them on every on-line media.
Targeted advertisement on baby products or pregnancy clothes can target
pregnant women as soon as they do a few searches on the web, or when the first
pregnancy tracking app is installed. Although targeted advertising algorithms
are intelligent enough to identify pregnant women in very early stages through
their search history, the algorithms fail to adapt their mechanism when a
pregnancy loss occurs. Despite all the search history on pregnancy loss
symptom and reasons, and activities such as joining support groups with
distinctive names, the algorithms continue to target these women . It will be
a false claim that recognising such situation is not possible despite all the
advancement in the field of social media behaviour tracking, emotion analysis,
and history of the on-line activities. Extensive effort and research has been
conducted on ethics in advertisement to regulate the content and reach of
advertising to protect the public over the past years. This includes ethical
issues in controversial advertising, matters of deception, representation,
targeting of vulnerable population. However, there is no particular research
on how these issues will be addressed when pregnant women are targeted and how
algorithms can be designed to behave differently when they figure out that a
loss has happened.

##  II ETHICAL ISSUES OF TARGETED ADVERTISING ON PREGNANT WOMEN

It is common amongst pregnant women to use on-line resources, pregnancy-
tracking Apps [1] and social media to access pregnancy-related information [2,
3]. Such resources remain as primary resources that women use even after a
pregnancy loss or miscarriage [4]. Pregnancy is not a similar experience for
all women, which can depend on many factors such as the family support, women
economic and social support, or the burden of their jobs. Both physical
changes and the whole unknown that a new child will bring to a woman’s life
sometimes makes pregnancy an overwhelming experience. The targeted advertising
can make this experience even worse by featuring glowing in shape pregnant
women, or pregnant mothers with their kids in clean and organised houses. Such
ads will start as soon as the first pregnancy tracking app is installed or
after a few searches on the internet, and will follow the pregnant women
recommending them Yoga classes, anti-stretchmarks body oils and towards the
end of the pregnancy showing off all the beautiful nurseries and baby
products. This might seem intelligent as the advertisement algorithm will keep
the product recommendation updated as the pregnancy progresses, however this
will be a torturous experience for a woman who has lost a pregnancy but keeps
receiving all of these updated targeted ads. Many of the general ethical
issues such as gender stereotypes [5], and female role portrayals [6, 7], of
targeted advertising on women apply to pregnant women as well. Although the
literature studies many of such ethical issues raised by advertisement
targeted at women, they lack the analysis of such advertisements targeted at
pregnant women. This includes both general products and specific maternity or
baby products targeted at pregnant women. The very important question to be
asked is why these ads keep following women experiencing a loss, despite the
fact that it is easy to identify individuals experiencing pregnancy loss or
miscarriage [8], either through their social media activity or their search
history. This question would not be raised in the first place if the
algorithms deciding who to target and what ad to show were designed ethically.
Algorithms are designed to learn about us and structure our lives in several
ways such as targeted advertisement, and determining our search results which
might follow a decision that we will make [9]. Such algorithms are
continuously learning about us through the available Big Data. For example,
they can know people’s ethnicity, gender and religion, if somebody is getting
a loan, has got fired [10], the places that has already visited or going to
visit, or even more personal information such as if getting married or
pregnant. Although algorithms can shape our lives in many beneficial ways,
their ethical implications must be discussed [11], especially for more
delicate cases such as advertisements targeted at pregnant women. Advertising
algorithms are getting more accurate having access to more and more data and
computational resources to process those data [12, 13]. They aim to target the
best advertisement to the best audience to maximize businesses’ profit.
However, they raise multiple ethical issues that has to be considered. As the
literature suggests there are several types of ethical concerns that can be
raised by algorithms:

  * •

Traceability, how a result of an algorithm can be traced to its source.

  * •

Transformative effects, how the original data is transformed to derive to a
particular result.

  * •

Unfair outcomes, an algorithm can produce biased results.

  * •

Misguided evidence, as the produced evidence will depend on the quality of the
input data.

  * •

Inscrutable evidence, as often it is not possible to understand the relation
between the data point being used and the drawn conclusion when machine
learning algorithms and big data are used.

  * •

Inconclusive evidence, as they are probable but uncertain knowledge that are
derived from machine learning or inferential statistics.

Additionally, the following concerns must be particularly addressed in this
context:

  * •

Ads’ timeliness, how relevant is the outcome of the algorithm during the time,
for example a pregnant woman might experience a very early miscarriage, should
the ads go on for months without taking into account any further data or user
activity?

  * •

Ads’ reach, who is seeing the ads? For example, Should the ads featuring
active glowing pregnant women be targeted at the ones experiencing
complications during their pregnancy? This also concerns about the reach of
the inappropriate ads to single mothers, or ones with weaker social or
economical situations.

##  III CONCLUSION AND FUTURE WORK:

In this paper, the ethical issues and the importance of ethical algorithm
design for target ads were briefly discussed. The targeted ads can become
helpful tools if they can recognise cases of loss, and recommending relevant
mind or body care for women experiencing a loss. This requires the algorithms
to always adapt and learn from the changes and use the feedback provided by
the user either directly or indirectly through search histories or online
activities.

## References

  * [1] D. Rodger, A. Skuse, M. Wilmore, S. Humphreys, J. Dalton, M. Flabouris, and V. Clifton, “Pregnant women’s use of information and communications technologies to access pregnancy-related health information in south australia,” _Australian journal of primary health_ , vol. 19, no. 4, pp. 308–312, 2013.
  * [2] C. Zhu, R. Zeng, W. Zhang, R. Evans, and R. He, “Pregnancy-related information seeking and sharing in the social media era among expectant mothers: qualitative study,” _Journal of medical Internet research_ , vol. 21, no. 12, p. e13694, 2019.
  * [3] F. Robinson and C. Jones, “Women’s engagement with mobile device applications in pregnancy and childbirth.” _The practising midwife_ , vol. 17, no. 1, pp. 23–25, 2014.
  * [4] N. Andalibi, G. Marcu, T. Moesgen, R. Mullin, and A. Forte, “Not alone: Designing for self-disclosure and social support exchange after pregnancy loss,” in _Extended Abstracts of the 2018 CHI Conference on Human Factors in Computing Systems_ , 2018, pp. 1–6.
  * [5] E. Plakoyiannaki, K. Mathioudaki, P. Dimitratos, and Y. Zotos, “Images of women in online advertisements of global products: does sexism exist?” _Journal of business ethics_ , vol. 83, no. 1, pp. 101–112, 2008.
  * [6] L. Tuncay Zayer and C. A. Coleman, “Advertising professionals’ perceptions of the impact of gender portrayals on men and women: a question of ethics?” _Journal of Advertising_ , vol. 44, no. 3, pp. 1–12, 2015.
  * [7] S. L. Grau and Y. C. Zotos, “Gender stereotypes in advertising: a review of current research,” _International Journal of Advertising_ , vol. 35, no. 5, pp. 761–770, 2016.
  * [8] P. Pang, M. Temple-Smith, C. Bellhouse, V.-H. Trieu, L. Kiropoulos, H. Williams, A. Coomarasamy, J. Brewin, A. Bowles, and J. Bilardi, “Online health seeking behaviours: what information is sought by women experiencing miscarriage,” _Stud Health Technol Inform_ , vol. 252, pp. 118–25, 2018\.
  * [9] K. Martin, “Ethical implications and accountability of algorithms,” _Journal of Business Ethics_ , vol. 160, no. 4, pp. 835–850, 2019.
  * [10] C. O’Neil, “Weapons of math destruction: How big data increases inequality and,” 2016.
  * [11] B. D. Mittelstadt, P. Allo, M. Taddeo, S. Wachter, and L. Floridi, “The ethics of algorithms: Mapping the debate,” _Big Data & Society_, vol. 3, no. 2, p. 2053951716679679, 2016.
  * [12] K. A. Thompson, “Commercial clicks: Advertising algorithms as commercial speech,” _Vand. J. Ent. & Tech. L._, vol. 21, p. 1019, 2018.
  * [13] J. Naor and D. Wajc, “Near-optimum online ad allocation for targeted advertising,” _ACM Transactions on Economics and Computation (TEAC)_ , vol. 6, no. 3-4, pp. 1–20, 2018.