Drivers of harmful algal blooms in coastal areas of Eastern Mediterranean: a machine learning methodological approach

1.   Introduction

Machine Learning (ML) algorithms are considered as powerful and reliable tools for prediction, applied so far in many scientific disciplines ^[1]. ML is a collection of algorithms of different types that improve prediction by gaining knowledge from data ^[2]. The main advantage of ML techniques is their ability to model multidimensional interactions between nonlinear, noisy, correlated or partly missing data. The field of marine ecology, characterized by complex interactions among a large number of variables, is a specific field counting several applications of ML algorithms ^[3].

Phytoplankton communities usually follow a rather typical annual pattern on biomass variations and species succession ^[4]. However, under a combination of forcing factors, it is possible that a regime shift in the system may be established. These factors can be the hydrodynamic conditions in the area ^[5], the buffering capacity of the system ^[6], salinity, temperature, nutrients, organic compounds ^[7], the presence of zooplankton and bacterial populations as well as species interactions ^[8]. However, work on those aspects ^[9], has shown that regime shifts are closely connected with the presence of nutrients; high sensitivity was also found in temperature fluctuations. These regime shifts in the case of phytoplankton are expressed as algal blooms. In many cases these blooms are due to exponential growth and dominance of toxic species of phytoplankton ^[10]. As the shift from the pristine conditions to an undesirable bloom does not follow a linear pattern, is of utmost importance to develop a probabilistic approach so that possible harmful effects on shellfish and fish culture as well as on human health can be indicated.

Algal blooms seem to be a problem of great concern among scientists due to the considerable economic, health and ecological effects on the neighboring coastal zones ^[11,12]. Τhe various causes of blooms, often linked with eutrophication in coastal ecosystems, is a matter under consideration ^[13,14]. Many methods have therefore been applied aiming at associating different parameters expressing land use ^[15], human activities, nutrient loadings, climate changes or other environmental processes with excessive algal growth ^[13,16,17]. These include the use of ML algorithms. In most cases, ML was applied to model and predict primary production (in terms of chlorophyll-α or biomass) using physical and chemical variables, as temperature, salinity, inorganic nitrogen and phosphorus concentrations ^[18,19,20]. In addition, there are studies that have modelled different dimensions of primary production such as diversity ^[21,22] or algal growth time ^[23]. A recent attempt has been made by Yu et al. (2021) ^[24] to predict phytoplankton blooms, using Machine Learning algorithms. Their argument is that "when planktonic algae proliferate over certain limits, harmful algal blooms (HABs) will occur". This is not necessarily the case ^[25]; increased concentrations of toxic algae can occur without a bloom and vice versa: not all blooms end to be toxic. Although algal blooms are entirely natural phenomena and have appeared throughout historical times ^[10], the frequency of occurrence has increased at an alarming level, over the last few decades. Land based sources of eutrophication, even when they are not the main drivers for triggering mechanisms for bloom formation, it is beyond any doubt that input of nutrients into the marine environment supports bloom formation. However, the principal question remains: will toxic microalgae grow? The answer can be given only if any predictions based on Machine Learning techniques, are targeted towards the exponential growth of toxic species, already recorded in the areas under study.

HABs are being considered as a specific form of coastal eutrophication, causing serious impact including death of organisms, toxins production, hypoxia and often human poisoning due to the consumption of infected fish or shellfish ^[14,26,27]. ML algorithms have been used to assess the formation of HABs in both fresh and marine waters ^[28]. Abundance as well as presence/absence of specific harmful algal species have been studied aiming to reveal possible associations with several biotic and abiotic variables. HAB forming species recently studied with selected ML techniques include Karenia brevis ^[29], Planktothrix rubescens ^[30] and Dinophysis acuminata ^[31], while Bourel et al. ^[32] have used ensembles i.e., combinations of classifiers to predict the presence-absence of 8 harmful marine phytoplankton species. Considering the various ML algorithms, special attention has been given to Neural Networks (NNs) and Random Forests (RFs), two commonly used algorithms with great performance even when dealing with complex ecological processes as HAB formation ^[33,34,35]. Some applications of NNs include the assessment of algal population dynamics ^[36,37,38], the classification of different algal types ^[39,40] and detection of HABs from satellite images ^[41,42]. On the other hand, the increased effectiveness of RFs has been denoted on the assessment of the distribution of phytoplankton cell size ^[43] and on the efficient prediction of chlorophyll-a concentration ^[44,45], of the abundance of some toxic genotypes ^[46] and of the presence of some toxic algal species ^[47].

Variables used as drivers in the application of ML algorithms for the prediction of the abundance or presence/absence of HAB forming microalgae, include both biotic and abiotic variables. Abiotic variables used so far include site, season, calendar day, distance from the coast, meteorological variables (air temperature, wind speed and direction, cloud cover), physical variables of seawater (temperature, salinity, conductivity, turbidity, oxygen saturation, pH, Secchi depth, photosynthetic radiation, remote sensing reflectance) and nutrients in both dissolved and particulate forms ^{[29,32,38,46]}. Considering biotic variables, chl-a is often used ^[24,29,47], as well as the abundances of microcrustaceans, ciliates, tintinids, microheterotrophs, cladocerans or copepodes ^[29,36]. It seems that the selection of cause variables in the existing literature mainly aims to improving the prediction of HABs, although not necessarily exploring a cause-and-effect relationship from the ecological point of view.

In the present study, a screening of several ML algorithms and ensembles is performed aiming to detect the presence/absence of 18 harmful or potentially harmful microalgae at genus level, using existing data collected from six coastal areas in the Aegean Sea, Greece, representing different productivity levels. The overall performance of each algorithm is assessed considering all the 18 studied microalgae. The cause variables include physical variables (temperature and salinity) as well as nutrients, already established in previous studies ^[16,20] as drivers of algal growth. The objectives of the study are: (a) to identify the most efficient algorithms or ensembles based on the overall prediction of presence/absence of the 18 microalgae, (b) to quantify the effect of cause variables on microalgal growth and (c) to attempt possible grouping of the studied HAB forming microalgae, considering their response to the abiotic drivers. Since all cause variables are easily measured on a routine basis, the results of this study can form the basis of an operational system for the prediction of HABs in coastal waters, eliminating possible adverse effects on ecosystem and human health.

2.   Materials and methods

2.1.   Study area and sources of data

The data set used in the present work was set up by compilation of six sets collected from 42 sampling sites in three different areas of the Aegean Sea, Greece. Five out of six sets originate from coastal waters whereas one set (Rhodes offshore) is characteristic of pelagic waters (Figure 1). The sampling areas are: (a) Kalloni Gulf (Island of Lesbos): eight sampling stations K1-K8 ^[48] (b) Gulf of Gera (Island of Lesbos): Sampling sites G1-G8 ^[49] (c) the coastal area of the city of Mytilini: sampling sites M1 and M2 ^[50] (d) Saronikos Gulf (near Athens): sampling sites S1-S9 ^[51] (e) offshore waters NW of the city of Rhodes R1-R5 ^[52] and (f) coastal waters of the city of Rhodes: ten sampling sites RH1-RH10 ^[51].

Figure 1.  Map of the six sampling areas.

DownLoad: Full-Size Img PowerPoint

2.2.   Variables used

The variables included in the dataset are: (a) temperature (T) and salinity (S) ^[20], also expressing seasonality and (b) nutrient concentrations, namely dissolved inorganic nitrogen (DIN), phosphates (${\mathrm{P}\mathrm{O}}_{4}^{3-}$) and silicates (${\mathrm{S}\mathrm{i}\mathrm{O}}_{2}$), often recognized as drivers of algal growth and increased primary productivity in coastal waters ^[53]. These variables are used as inputs (cause variables) in the ML algorithms. The target variable trying to estimate is the presence or absence of 18 genera of harmful or potentially harmful algae estimated by the list of toxic microalgae of IOC-UNESCO ^[54], shown in Table 1. The number of occurrences of each genus in the 889 samples, the sites where found and summary statistics of the abiotic conditions under which each genus appeared, are shown in Table 1. Redfield' s ratio, i.e., nitrogen to phosphorus ratio (N:P), a characteristic proxy of the nutrient limiting primary productivity, is also shown in Table 1.

Table 1.  The studied HAB forming microalgae, their occurrence and summary statistics (mean and range) of the corresponding abiotic variables related to each particular genus.

Genus	Occurrences	Sampling Areas	T (℃)	S (psu)	DIN (μM)	PO43- (μΜ)	SiO2 (μΜ)	N:P
Alexandrium	25	K, R1	14.76 (9.8-23.9)	37.81 (34.0-40.1)	5.78 (0.5-39.8)	0.12 (0.0-1.5)	15.97 (3.4-86.6)	196.64 (7.8-1217.5)
Amphidinium	71	K, S, R1, R2	17.04 (9.8-23.9)	38.26 (34.0-40.2)	4.78 (0.5-39.8)	0.13 (0.0-1.5)	14.34 (2.4-86.6)	232.41 (1.7-6020.0)
Cryptomonas	168	K	17.10 (9.5-27.7)	38.57 (34.0-41.1)	4.13 (0.5-45.2)	0.09 (0.0-1.6)	14.36 (1.7-94.0)	180.16 (6.81-6020.0)
Dictyocha	48	K	13.93 (9.4-23.3)	38.29 (34.8-40.7)	4.46 (0.6-33.2)	0.10 (0.0-1.3)	16.22 (4.3-70.9)	102.69 (7.4-692.2)
Dinophysis	228	K, G, M, S, R1, R2	18.27 (9.7-28.2)	38.29 (28.3-40.6)	2.70 (0.1-38.0)	0.16 (0.0-1.8)	8.85 (0.4-80.7)	43.76 (2.8-1712.0)
Goniodoma	21	G, M, S, R1, R2	19.20 (13.4-27.4)	38.35 (37.3-39.1)	2.84 (0.4-11.7)	0.45 (0.0-6.0)	9.20 (0.8-15.1)	25.41 (1.5-116.7)
Gonyaulax	37	K, G, M, S, R1, R2	18.41 (14.37-22.11)	38.82 (36.4-40.6)	1.66 (0.3-4.1)	0.19 (0.0-0.85)	7.19 (2.1-20.9)	30.05 (2.7-379.0)
Gymnodinium	279	K, M, S, R1, R2	19.58 (9.6-27.7)	38.83 (28.0-40.2)	1.76 (0.1-27.8)	0.07 (0.0-4.1)	9.66 (1.6-53.9)	68.13 (0.2-1660.0)
Gyrodinium	177	K, G, M, S, R1, R2	18.23 (9.7-27.7)	38.51 (34.7-40.3)	2.78 (0.4-45.2)	0.16 (0.0-1.8)	8.63 (1.5-94.0)	70.01 (1.7-1660.0)
Karenia	394	S, R1, R2	20.00 (13.1-27.6)	38.55 (28.0-39.7)	2.01 (0.1-38.0)	0.15 (0.0-6.0)	8.02 (0.4-53.9)	33.04 (1.4-1643.0)
Karlodinium	133	K	17.20 (9.9-27.7)	38.68 (34.8-40.7)	3.51 (0.5-33.2)	0.07 (0.0-1.3)	12.64 (1.7-70.9)	180.74 (6.8-6020.0)
Lingulodinium	23	K, R1	14.81 (9.9-25.6)	37.67 (34.8-40.0)	6.00 (1.4-33.2)	0.10 (0.0-1.3)	11.81 (2.4-68.7)	282.78 (9.6-1660.0)
Navicula	165	G, M, S, R1, R2	19.27 (10.0-26.5)	39.00 (37.1-40.3)	2.22 (0.3-23.8)	0.14 (0.0-4.1)	9.22 (1.1-53.9)	44.78 (0.2-1643.0)
Peridinium	301	K, G, M, S, R1, R2	18.92 (9.6-27.7)	38.62 (28.0-40.3)	2.60 (0.1-39.8)	0.11 (0.0-2.3)	9.87 (1.2-86.6)	71.58 (1.7-1712.0)
Phaeocystis	71	S, R1	16.4 (13.1-26.7)	38.27 (37.6-39.1)	2.98 (0.5-14.1)	0.36 (0.0-6.0)	5.44 (0.4-17.2)	16.24 (1.4-72.5)
Polykrikos	38	K, M, R1	14.88 (10.0-25.6)	38.7 (34.9-40.4)	3.77 (0.6-22.8)	0.06 (0.0-0.9)	9.56 (1.1-70.9)	346.81 (11.9-6020.0)
Prorocentrum	777	K, G, M, S, R1, R2	19.10 (9.4-28.2)	38.61 (28.0-41.1)	2.26 (0.1-38.0)	0.13 (0.0-6.0)	8.43 (0.4-83.0)	51.69 (1.4-1660.0)
Pseudo-nitzschia	491	K, G, M, S, R1, R2	17.39 (9.4-27.7)	38.5 (28.3-41.1)	2.89 (0.4-45.2)	0.14 (0.0-6.0)	8.92 (0.4-94.0)	78.51 (1.4-1712.0)

---

 | Show Table

DownLoad: CSV

2.3.   Machine learning algorithms

A variety of ML algorithms was used to classify the presence or absence of the 18 genera of HAB forming microalgae (Table 2). The selected classifiers cover all basic supervised ML categories i.e., rules, trees, lazy, functions and bayes, whereas classic ensemble methods, already used in marine ecology ^[3], were also applied. In order to optimize the performance of the basic classifiers, different values of the crucial (hyper) parameters of each algorithm were tested during training. Thus, the training included at least four different versions of each basic algorithm (Table 2, last column) eliminating the relation of the parameters' selection to the final algorithm performance. Random Forest (RF) was considered as an ensemble, since it uses a set of decision trees to provide classification. For some ensembles (e.g., Voting), which use techniques trained by the combination of the results of basic classifiers, the number of combined classifiers was limited to three in order to avoid ties. The selected three classifiers were the three best performing basic classifiers, as previously defined.

Table 2.  ML algorithms used in the current study: category, abbreviation, short description and hyper parameter values.

Category	Abbreviation	Algorithm Description	Hyper parameter values
Rules	Jrip	Repeated incremental pruning to produce error reduction (RIPPER)	Batch size = 50,100 Min total weight = 2, 5
	Part	PART decision list	Batch size = 50,100 Min number of instances per rule = 2, 5
Trees	J48	C4 pruned decision tree	Batch size = 50,100 Min number of instances per leaf = 2, 5 Number of folds = 3, 5
	Rep	Decision tree using reduced error pruning with backfiting	Batch size = 50,100 Min number of instances per leaf = 2, 5
Lazy	IBk	The k-nearest neighbors using Euclidean distance	Number of neighbors = 1, 5, 10, 20
	KStar	Nearest neighbor with entropic distance	Global blend = 5, 10, 20, 50
Functions	Log	Multinomial logistic regression
	MLP	Multilayer Perceptron using backpropagation	Number of neurons = 2, 5, 7, 10
Bayes	NB	Naïve Bayes using estimator classes	Batch size = 50,100
	BN	Bayes network	Batch size = 50,100
Ensembles	RF	Forest of random decision trees	Batch size = 50,100 Number of iterations = 20, 50,100
	Bagging	Bagging classifiers to reduce variance	Classifiers: The best basic one Number of iterations = 10
	Boosting	Boosting classifiers using Adaboost M1 method	Classifier: The best basic one
	CVR	Classification using regression methods	Classifier: M5 model tree
	RC	Randomizabling classifiers (Random Committee)	Classifier: Random tree
	Stacking	Combining classifiers using stacking method	Classifiers: The 3 best basic ones Number of folds = 10
	Voting	Combining classifiers using votes	Classifiers: The 3 best basic ones Combination rule: Average of Probabilities

---

 | Show Table

DownLoad: CSV

2.4.   Algorithm evaluation and selection of abiotic parameters

The efficiency of algorithms was evaluated using the 10-fold cross validation (10-fold CV) procedure ^[55]. The RWeka interface written in R ^[56], of Weka ML techniques ^[57], was used to run and test the algorithms. The number of correctly classified instances of either the presence or absence of each genus as labels of the target class over the total number of water samples (i.e., accuracy), was used to determine the performance effectiveness of each algorithm. Additionally, some other measures of predictive performance were used in order to better evaluate the algorithms' predictions ^[58]. More precisely, sensitivity (or recall) expressing algorithm completeness, is the fraction of the correct genus presences over the total number of predicted presences (true presences plus false absences) in the total samples. A similar measure for the prediction of absences was used as a second measure of completeness, since the conditions related to the absence of a genus in a sample are also considered crucial ^[59]. This measure is specificity which is the same fraction calculated on absences (number of correct predicted absences over the sum of true absences plus false presences). Moreover, precision was used, expressing the power of algorithm's correctness as it measures in how many instances that the algorithm predicted as "genus present", the genus was actually present. The classic measure of kappa statistic that represents the degree of accuracy and reliability in classification problems, was also included ^[60]. Kappa statistic ranges from -1 (total disagreement) through 0 (random classification) to 1 (perfect agreement). Finally, due to the existence of imbalanced data, i.e., low number of appearances for eight genera (i.e., Alexandrium, Amphidinium, Dictyocha, Goniodoma, Gonyaulax, Lingulodinium, Phaeocystis and Polykrikos) in the dataset (see Table 1), one more specialized evaluation measure was used; the discriminant power is a measure that combines sensitivity and specificity and evaluates how well an algorithm distinguishes the presences and absences of a rare genus in case of imbalanced data ^[59]. The thresholds for this measure are: "poor" for values less than 1, "limited" for values between 1 and 2, "fair" for values between 2 and 3 and "good" for values higher than 3 ^[61].

The performance of an algorithm can be crucially affected by the variables included in the input dataset ^[24]. In order to select the optimal subset of the input abiotic variables offering the higher predictive performance to each trained algorithm, an exhaustive search was implemented. All possible subsets of the five available abiotic variables containing 5 (all input variables), 4 out of 5 (5 subsets), 3 out of 5 (10 subsets) and 2 out of 5 (10 subsets) variables were used to find the best combination.

2.5.   Estimation of the importance of the cause variables

Some ML algorithms estimate directly the importance of the input variables (e.g., most rules or trees). Especially for tree algorithms, applying the split-criterion in each node, the relative importance of each input variable can be easily determined. The most common method to measure a node's effectiveness in trees is the Gini Impurity which measures the probability of misclassification for a new instance (i.e., sample) by its specific tree node. A related measure for RFs is the Mean Decrease of Impurity (MDI) which weighs the impact of each input variable to the final prediction (presence of a genus) by measuring the effectiveness of each variable at reducing the uncertainty during tree induction ^[62]. In this study the MDI measure was used to assess the importance of each abiotic variable to the presence or absence of the potentially harmful algae genera. Furthermore, this measure was also used for grouping genera affected by similar abiotic conditions. Hierarchical clustering was applied using Euclidean distance and Ward's minimum variance method.

3.   Results

3.1.   Performance of basic classifiers and ensembles

The performance of the ML algorithms measured by the accuracy of correctly classified presences/absences in the 889 seawater samples is presented in Table 3 while Table 4 contains the best performing algorithm along with the optimal values of hyper parameters and the optimal subset of the abiotic variables. The predictive performance for all genera is above 70%, being over 90% for 10 out of 18 genera. The higher percentage of correctly classified instances was for Lingulodinium (98.1%) and the lowest for Peridinium (72.2%). The performance of the various algorithms did not vary significantly within each genus, but mainly among genera. Indeed, the range of predictive performance among different algorithms for the same genus ranged from 32.2% for Gymnodinium to 1.4% for Goniodoma (9.1% the mean difference for all genera), while the corresponding range among genera for the same algorithm ranged from 44.2% for Naïve Bayes to 25.8% for RF technique (30.6% mean difference for all algorithms).

Table 3.  Predictive performance (best performing single and ensemble in bold) of different ML techniques in terms of accuracy in HA genera presence/absence prediction as evaluated by 10-CV.

	Rules		Trees		Lazy		Functions		Bayes		Ensembles
HA genus	Jrip	Part	J48	Rep	IBK	KStar	Logistic	MLP	NB	BN	RF	Bagging	Boosting	CVR	RC	Stacking	Voting
Alexandrium	0.970	0.973	0.972	0.972	0.976	0.968	0.971	0.972	0.957	0.965	0.973	0.972	0.968	0.972	0.972	0.972	0.972
Amphidinium	0.918	0.919	0.916	0.917	0.916	0.917	0.919	0.916	0.907	0.910	0.925	0.921	0.917	0.919	0.922	0.920	0.917
Cryptomonas	0.925	0.885	0.899	0.873	0.749	0.907	0.858	0.894	0.838	0.910	0.936	0.936	0.942	0.882	0.904	0.811	0.940
Dictyocha	0.954	0.943	0.947	0.943	0.955	0.953	0.943	0.943	0.936	0.948	0.965	0.946	0.946	0.940	0.943	0.946	0.937
Dinophysis	0.759	0.771	0.756	0.748	0.741	0.727	0.727	0.744	0.731	0.780	0.769	0.760	0.780	0.767	0.736	0.744	0.770
Goniodoma	0.976	0.976	0.976	0.976	0.979	0.970	0.975	0.975	0.965	0.976	0.976	0.978	0.979	0.977	0.970	0.976	0.976
Gonyaulax	0.957	0.955	0.958	0.958	0.955	0.954	0.956	0.958	0.926	0.926	0.961	0.960	0.961	0.958	0.953	0.958	0.958
Gymnodinium	0.772	0.780	0.793	0.772	0.692	0.803	0.680	0.742	0.524	0.753	0.846	0.811	0.773	0.812	0.802	0.686	0.827
Gyrodinium	0.793	0.802	0.802	0.800	0.729	0.801	0.800	0.798	0.775	0.799	0.827	0.812	0.790	0.799	0.790	0.801	0.808
Karenia	0.813	0.762	0.757	0.765	0.779	0.791	0.588	0.729	0.650	0.666	0.848	0.791	0.794	0.772	0.799	0.683	0.805
Karlodinium	0.901	0.879	0.880	0.889	0.748	0.898	0.857	0.884	0.852	0.916	0.917	0.915	0.916	0.899	0.911	0.850	0.916
Lingulodinium	0.980	0.981	0.981	0.978	0.974	0.972	0.973	0.974	0.966	0.974	0.980	0.976	0.975	0.974	0.981	0.980	0.981
Navicula	0.812	0.818	0.814	0.805	0.803	0.792	0.814	0.832	0.791	0.814	0.822	0.814	0.800	0.838	0.799	0.814	0.814
Peridinium	0.714	0.706	0.687	0.685	0.586	0.679	0.658	0.674	0.654	0.650	0.722	0.719	0.717	0.708	0.690	0.661	0.717
Phaeocystis	0.947	0.944	0.945	0.952	0.951	0.934	0.918	0.961	0.916	0.951	0.958	0.961	0.956	0.947	0.954	0.920	0.961
Polykrikos	0.957	0.962	0.961	0.961	0.964	0.955	0.958	0.958	0.945	0.958	0.962	0.961	0.958	0.962	0.958	0.960	0.963
Prorocentrum	0.876	0.875	0.874	0.875	0.876	0.880	0.874	0.874	0.862	0.874	0.877	0.877	0.872	0.875	0.854	0.874	0.874
Pseudo-nitzschia	0.757	0.756	0.761	0.754	0.771	0.765	0.665	0.699	0.565	0.737	0.791	0.768	0.771	0.762	0.750	0.653	0.786
Mean	0.877	0.872	0.871	0.868	0.841	0.870	0.841	0.863	0.820	0.862	0.892	0.882	0.879	0.876	0.872	0.845	0.885

---

 | Show Table

DownLoad: CSV

Table 4.  Best perfomed algorithme along with the optimal parameter combination for each genera.

HA genus	Best Algorithm	Hyper Paremeters values	Best Parameter Combination
Alexandrium	IBk	Number of neighbors=5	T+S+PO43-
Amphidinium	RF	Batch size=50 Number of iterations=50	T+S+DIN+ PO43-+SiO2
Cryptomonas	Boosting	Classifier: Jrip (Batch size=100 Min total weight=5)	T+S+DIN+ PO43-+SiO2
Dictyocha	RF	Batch size = 100 Number of iterations = 50	T + S + PO4 + SiO2
Dinophysis	BN/ Boosting	Batch size = 100/ Classifier: IBk (Number of neighbors = 10)	T + S + DIN + PO43- + SiO2
Goniodoma	IBk	Number of neighbors = 10	T + S + DIN + PO43- + SiO2
Gonyaulax	RF	Batch size = 100 Number of iterations = 100	T + S + DIN + PO43- + SiO2
Gymnodinium	RF	Batch size = 50 Number of iterations = 50	T + S + DIN + PO43- + SiO2
Gyrodinium	RF	Batch size = 50 Number of iterations = 100	T + S + DIN + PO43- +SiO2
Karenia	RF	Batch size = 50 Number of iterations = 100	T + S + PO43- + SiO2
Karlodinium	RF	Batch size = 50 Number of iterations = 100	T + S + DIN+ PO43- +SiO2
Lingulodinium	Part/ RC/ Voting	Batch size = 50, Min number of instances per rule =2/ Classifier: Random tree/ Classifiers: Jrip (Batch size = 100, Min total weight = 5), Part (Batch size = 100, Min number of instances per rule = 2), J48 (Batch size = 100, Min number of instances per leaf = 5)	T + S + DIN + PO43- +SiO2
Navicula	CVR	Classifier: M5 model tree	T + S + DIN + PO43- + SiO2
Peridinium	RF	Batch size = 50 Number of iterations = 100	T + S + DIN + PO43- + SiO2
Phaeocystis	MLP/Bagging/Voting	Number of neurons = 7/ Classifier:MLP/ Classifiers: MLP, Rep (Batch size = 50, Min number of instances per leaf = 2), BN (Batch size = 100)	T + S + DIN + PO43- + SiO2
Polykrikos	IBk	Number of neighbors = 10	T + S + PO43-
Prorocentrum	KStar	Global blend = 50	T + S + PO43-
Pseudo-nitzschia	RF	Batch size = 50 Number of iterations = 100	T + S + DIN + PO43- + SiO2

---

 | Show Table

DownLoad: CSV

Among the basic algorithms, Jrip showed the best performance with mean predictive percentage for all genera equal to 87.7%. Part rule algorithm, both trees, KStar, MLP and BN also showed high performance exceeding 86%. Although IBK and MLP are the best performing algorithms regarding 6, and 3 genera respectively, they achieve lower overall performance compared to other basic algorithms. Finally, the classical statistical procedure of Logistic function and NB give predictions of rather low quality.

As expected, the ensembles generally exceed the performance of single classifiers, with four ensembles (RF, Bagging, Boosting, and Voting) achieving higher performances than the best basic, when considering means for all genera. The most efficient ensemble is RF with 89.2% mean predictive performance (1.5% higher than Jrip), providing the highest prediction for 9 genera (Amphidinium, Dictyocha, Gonyaulax, Gymnodinium, Gyrodinium, Karenia, Karlodinium, Peridinium and Pseudo-nitzschia). The second-best performing ensemble is voting with almost 0.7% lower mean predictive percentage (88.5%) compared to RF, offering the most successful prediction for Phaeocystis. Both Bagging and Boosting also performed sufficiently well, the former achieving an overall performance above 88.2% and the latter providing the best prediction for 4 genera. CVR, RC and Stacking show the same or lower performance compared with the basic algorithms and thus are considered as inefficient. For most genera (13 out of 18), the performance of the algorithms was optimal when all variables were included in the input data set. For Dictyocha and Karenia the best performance was achieved when removing DIN, whereas for the two Lazy algorithms (IBk and KStar) the prediction for Alexandrium, Polykrikos and Prorocentrum was optimal when T, S and ${\mathrm{P}\mathrm{O}}_{4}^{3-}$ were only included in the input variables. For RF in particular, the use of all five abiotic variables optimized the prediction of 7 genera (Ampidinium, Gonyaulax, Gymnodinium, Gyrodinium, Karlodinium, Peridinium and Pseudo-nitzschia). Therefore in the rest of the analysis aiming to propose a single algorithm for all genera and to assess the relative importance of the input variables for prediction, RF the best performing ensemble was applied on the initial input data set of all five abiotic variables.

3.2.   Assessing the performance of RF

Since RF seems to be the overall best performing algorithm, more details based on RF classification results using the initial dataset (optimal parameter combination) are shown in Table 5. The mean precision is 0.723, meaning that in 72.3% of the cases in which RF predicted the presence of a genus, the genus is actually present. For some genera this measure is rather high, being 91.5% for Karlodinium and 88.1% for Prorocentrum, showing an upward trend for genera with high number of occurrences in the whole database. On the other hand, the rate of correct classifications of presences that is sensitivity is rather moderate since only 43.1% of the total genus presences is classified correctly by the algorithm. The predictions are even worse for genera present in less than 10% of the samples (rare species), as Alexandrium, Goniodoma, Gonyaulax or Polykrikos. The specificity of the RF classification (correctly predicted absences) is generally high with a mean of 92.1%. This finding implies that absences are more correctly classified than presences, especially for rare genera. On the other hand, the mean value of Kappa statistic is equal to 0.469 showing a rather moderate predictive performance of the RF algorithm. The value of Kappa statistic is lower for the rare genera and varies from 0.287 (fair agreement) for Gonyaulax to 0.723 (substantial agreement) for Cryptomonas. However, the discriminant power, a special measure for the predictive performance for rare genera, is greater than 2 for 6 rare genera, showing that the prediction by RF is "fair" in these cases. Considering each specific genus, the prediction is "good" being above 3 for Lingulodinium and "limited" being 1.632 for Amphidinium.

Table 5.  Predictive details of RF (the best performing algorithm) for each genus. The discriminant power refers only to genera with imbalanced data (rare genera).

Algae genus	N# presences	Precision	Sensitivity	Specificity	Kappa statistic	Discriminant power
Alexandrium	25	0.556	0.200	0.995	0.312	2.154
Amphidinium	71	0.567	0.239	0.984	0.386	1.632
Cryptomonas	168	0.878	0.774	0.974	0.723
Dictyocha	48	0.615	0.333	0.988	0.503	2.049
Dinophysis	228	0.647	0.407	0.912	0.500
Goniodoma	21	0.667	0.095	0.999	0.324	2.565
Gonyaulax	37	0.667	0.108	0.998	0.287	2.261
Gymnodinium	279	0.814	0.659	0.931	0.655
Gyrodinium	177	0.682	0.353	0.970	0.509
Karenia	394	0.758	0.787	0.800	0.604
Karlodinium	133	0.915	0.489	0.992	0.595
Lingulodinium	23	0.833	0.217	0.999	0.291	3.100
Navicula	165	0.550	0.200	0.963	0.480
Peridinium	301	0.642	0.405	0.884	0.518
Phaeocystis	71	0.867	0.549	0.993	0.633	2.840
Polykrikos	38	0.667	0.158	0.996	0.349	2.119
Prorocentrum	777	0.881	0.985	0.429	0.375
Pseudo-nitzschia	490	0.811	0.808	0.769	0.668
Mean values		0.723	0.431	0.921	0.469	2.340

---

 | Show Table

DownLoad: CSV

3.3.   Weighing the effect of abiotic variables

According to our findings, RF seems to be the most efficient algorithm among single classifiers and ensembles. Since RF is also showing a general immediacy in weighing the effects of input variables based on the composition of the constructed trees, it is selected for providing further insights on the role of the abiotic drivers for the presence/absence of potentially harmful microalgae.

The relative importance of each abiotic variable in the trees constructed for each genus by the RF algorithm, is assessed with the MDI measure (Figure 2). Temperature (T) has the most powerful effect with an MDI median of 0.325, whereas its importance varies among genera (minimum of 0.12 for Pseudo-nitzschia and maximum of 0.40 for Prorocentrum). On the other hand, salinity (S) seems to have the lowest effect, with an MDI median equal to 0.169 (minimum of 0.11 for Amfidinium and maximum of 0.19 for Lingulodinium). ${\mathrm{P}\mathrm{O}}_{4}^{3-}$is the most important variable among nutrients, with a median of 0.26, against 0.21 and 0.22 for DIN and ${\mathrm{S}\mathrm{i}\mathrm{O}}_{2}$, respectively. Moreover, the importance of nutrients did not vary considerably among genera (standard deviation equal to 0.39, 0.39 and 0.34 for DIN, ${\mathrm{P}\mathrm{O}}_{4}^{3-}$and ${\mathrm{S}\mathrm{i}\mathrm{O}}_{2}$, respectively), implying the similar effect of each nutrient on all species considered.

Figure 2.  Box-and-whisker plots of MDI, measuring the relative importance of the various abiotic variables on the construction of RF trees.

DownLoad: Full-Size Img PowerPoint

3.4.   Clustering of genera according to the similarity of the effects of abiotic variables

Possible similarity in the importance of the effects of abiotic variables among the various genera in terms of MDI is assessed with hierarchical clustering (Figure 3). It seems that two groups are formed when considering the effect of the abiotic variables on trees' construction. The first cluster is formed off 8 genera, being Pseudo-nitzschia, Dictyocha, Peridinium, Phaeoystis Gymnodinium, Cryptomonas, and two genera belonging to the Brachidiniales (BR) order (Karenia and Karlodinium). The second cluster includes 10 genera, i.e., Dinophysis, Navicula, Prorocentrum, genera of the Gonyaulacales (GO) order (Alexandrium, Goniodoma, Gonyaulax and Lingulodinium) and of the Gymnodiniales (GY) order (Amphidinium, Gyrodinium and Polykrikos). The members of the first cluster are generally less affected by the abiotic variables compared to the members of the second cluster. Nutrients seem to have the main role for the occurrence of the members of the first cluster, compared to the physical variables (temperature and salinity). On the other hand, the members of the second cluster are mostly affected by temperature and nutrients, whereas salinity has the lowest effect.

Figure 3.  (a) Hierarchical clustering tree of the genera based on the relative importance of the abiotic variables based on MDI (Mean Decrease of Impurity) (b) The relative importance of the abiotic variables for each cluster. Brachidiniales (BR), Gonyaulacales (GO) and Gymnodiniales (GY) orders.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

In this study, we assessed the overall performance of various ML techniques (both single and ensemble) for the prediction of the presence/absence of 18 harmful or potentially harmful marine microalgae at genus level. We further optimized the performance of the algorithms using different values of initial parameters using 10-fold CV. Jrip was the best performing basic classifier implying that single rules involving the five input variables can rather effectively model the response variable. The effectiveness of Jrip single rules has been also identified in previous studies ^[63,64]. The best ensemble method was RF achieving an overall mean accuracy of 89.2% and outperforming the rest algorithms by at least 1.5% in accuracy. This is in accordance with many previous studies that revealed the exceptional predictive power of RFs ^[65,33], including marine ecological studies ^[66,44]. A possible explanation for the successful performance of RFs is that they retain the benefits of the embedded decision trees, while they also combine tree results exploiting the great effectiveness of the majority of voting schema ^[67]. Thus, the recurring discrimination of samples based on the abiotic variables during the tree induction process crucially supports the detection of genera presences. The bootstrapping methods used for the construction of different trees protect RFs from overfitting issues while the majority of voting schema combines the predictions by taking into advantage the overlapped conditions that drive genus's occurrence in each tree ^[68]. Results on accuracy were remarkably high (89.2%), although the detailed analysis of the RF results revealed some predicting weaknesses. The sensitivity of the model was found rather moderate, implying that some real presences were not correctly identified by the algorithm. On the other hand, the specificity was above 90% indicating that absences were successfully determined. The precision of the predictions was found relatively high (72.3%), implying that the predicted presences are in a high percentage real presences. High specificity combined with moderate sensitivity and low kappa statistic is a weakness of RFs, usually arising when occurrences are rare ^[69,70]. Data bootstrapping is biased towards the major class (absences) leading RF technique to over-predict the major and under-predict the minor class ^[71]. This imbalance was observed in the present study since 8 genera (i.e., Alexandrium, Amphidinium, Dictyocha, Goniodoma, Gonyaulax, Lingulodinium, Phaeocystis and Polykrikos) were only present in less than 8% of the total samples. Discriminant power for these genera showed at least a "fair" prediction for RFs with the exception of Amphidinium for which prediction was characterized as "limited". This finding combined with the higher mean sensitivity for the other more common genera (being 58.7%), which improves as the balance between presences and absences increases (80.8% and 78.7% for the genera with the most presences Pseudo-nitzschia and Karenia, respectively) render RF results as satisfying. Therefore, special attention should be paid to the prediction of occurrences of rare genera using RFs or other ML techniques. This can be improved by using specific performance measures or by appropriate manipulations on the dataset, as oversampling or repeated random sampling ^[72]. The exhaustive search for the optimal subset of abiotic variables to be used for input in the ML algorithms, showed that the best performance was achieved when all five variables were used for most genera. As it was found in a similar study ^[24], a low number of input variables improves prediction, however it seems that in our case the number of five variables is already low and the removal of variables results to information loss. Aiming to exploit the advantages of RFs, we weighted the relative importance of the input variables (i.e., T, S, DIN, ${\mathrm{P}\mathrm{O}}_{4}^{3-}$ and ${\mathrm{S}\mathrm{i}\mathrm{O}}_{2}$) using the initial dataset, in terms of MDI, for the presence/absence of the 18 harmful microalgal genera. Temperature was the most important driver of occurrence, although its power varied among genera. The key-role of temperature for the appearance and abundance of phytoplankton species has been recognized by various studies in Mediterranean Sea ^[20,73,48] and its effect is considered as mainly indirect, associated with seasonal changes and stratification which are identified as drivers of primary production in both ocean and coastal waters ^[74]. Due to the recent climate changes and the tendency of seawater temperature to increase, the effect of temperature is currently under thorough investigation ^[75]. It has been reported since a long time ^[76] that higher temperatures favor the growth of flagellates, whereas diatoms are well adapted at lower temperatures. According to Wyatt ^[77] "diatoms have seasonal standing-crop maxima in spring and autumn in middle latitudes while lower numbers occur during the intervening months". On the contrary blooms of dinoflagellates are more characteristic during summer time, especially in temperate and subtropical seas. The latter favors dinoflagellate dominance over diatoms during summer period. It is obvious that algal species succession is highly depended on temperature changes and this fact is clearly reflected in the results of the present work.

Although temperature is an important driver of HABs, salinity had the less significant role (lower mean MDI value) for the studied microalgae. There are two possible explanations for this lack of significance. Salinity acts on marine organisms through the osmotic pressure exercised on their cellular fluids. In the present work the range of salinity values was fairly limited (around 38 psu) and therefore it cannot be concluded whether this variable is not of importance to phytoplankton or the effects were not shown due to the narrow range of the salinity values used in this work.

Nutrients were also found to affect the presence of the studied microalgae, phosphates having a principal effect. The mean N:P values referring to specific genera are much higher than 16 with only one exception: the microalga Phaeocystis (with a corresponding N:P value 16:24). This value indicates that the conditions were phosphorus limited (Table 1). Nutrient concentrations in the surface waters in the Eastern Mediterranean are low due to mixing processes with nutrient poorer basin water and biological activity ^[78]. In addition to this shortcoming that limits phytoplankton abundance, the Eastern Mediterranean Sea seems to be phosphorus limited ^[79]. This is a rather peculiar characteristic, since the oceans, including the Atlantic, are nitrogen limited and the Atlantic Ocean supplies the Mediterranean with water masses through the Strait of Gibraltar. The interpretation of the results in the present work is also complicated because there is an additional problem regarding phosphates. There are some algal species that given a supply of phosphates, they can store it within their cells in the form of polyphosphate (volutin) granules ^[7]. This "luxury phosphorus" as it is known in the literature, can be consumed by the algae and support algal growth in the absence of an external supply source. This way any relationship between phosphorus concentration dissolved in the marine environment and algal standing stock is weakened due to the presence of stored phosphorus in some species.

According to the MDI values, expressing the effect of abiotic drivers for the appearance of the studied harmful microalgae, a possible grouping of those organisms was attempted. Two clusters were formed, the first of eight genera (Pseudo-nitzschia, Dictyocha, Peridinium, Phaeoystis Gymnodinium, Cryptomonas, Karenia and Karlodinium) and the second of ten genera (Dinophysis, Navicula, Prorocentrum, Alexandrium, Goniodoma, Gonyaulax, Lingulodinium, Amphidinium, Gyrodinium and Polykrikos). The first group is rather diverse, including members from several phyla, as Bacillariophyta (Pseudo-nitzschia), Ochrophyta (Dictyocha), Miozoa (Peridinium, Gymnodinium, Karenia, Karlodinium), Haptophyta (Phaeocystis) and Cryptophyta (Cryptomonas). The second group mainly includes members of the Miozoa phylum, except of Navicula belonging to Bacillariophyta. Therefore, common functional characteristics of genera belonging to the same phyla seem to play a role for their appearance ^[80], however this role is not so clear-cut. Members of both clusters are influenced by nutrients, and mainly by phosphorus, which is possibly related to the phosphorus limitation often observed in Eastern Mediterranean waters. For some of the studied genera, the role of phosphorus as a driver for presence of some genera can be found in the existing literature, as Karlodinium in Li et al. ^[81] and Cryptomonas in Gasol et al. (1993) ^[82]. Considering the physical drivers, salinity plays a minor role for both groups, although its role is important for specific microalgae as Pseudo-nitzchia ^[83,84] and Karenia ^[85]. Temperature is the main driver for the presence of the members of the second group, mainly involving Miozoa. Changes of temperature express seasonality, so related processes as stratification or factors as light availability ^[86] may be also important for the proliferation of the members of the second group.

Results from this study showed that it is possible to anticipate harmful algal blooms and design management practices to mitigate their effects on marine life and humans. This could be done by monitoring abiotic drivers that were identified in this work and issue appropriate warnings that may suggest some sort of action on behalf of appropriate management authorities. This study and the proposed methodology may form the basis for an effort to improve predictability of these occasionally devastating events.

5.   Conclusions

In the present study various ML techniques were assessed for their efficiency to predict the appearance of 18 potentially harmful marine microalgae using data from 6 coastal sites in Eastern Mediterranean. RF algorithm using five abiotic drivers was the most efficient algorithm for prediction when considering overall the 18 studied microalgae. Moreover RF results were useful to enlighten the role of the various abiotic drivers in an ecological context. Although the overall performance of RF was satisfying in terms of predictability, a more exhaustive training of the algorithm with large number of samples is always desirable. Since the five abiotic variable are easily measured in a routine basis, the proposed methodology may form the basis of an operational system to be used for the prediction of HABs and therefore eliminate effects on marine life and human health.

Acknowledgments

This research is funded in the context of the project "Predicting harmful algae blooms using machine learning techniques" (MIS 5047187) under the call for proposals "Supporting researchers with an emphasis on new researchers" (EDULLL 103). The project is co-financed by Greece and the European Union (European Social Fund- ESF) by the Operational Programme Human Resources Development, Education and Lifelong Learning 2014–2020.

The authors wish to thank two anonymous reviewers for their insightful comments on a previous version of the manuscript.

Conflict of interest

There is no conflict of interests.