{admonition}
:class: tip
This chapter was written and developed by {{ barros }} <img src="../img/authors/ricardo.jpg" alt="Ricardo Barros" width="50" height="50">

Support Vector Machine (SVM)¶

This tutorial features the Support Vector Machine (SVM) technique for classifying riverbed morphology types as a function of hydraulic and sediment transport data.

{admonition}
+ Pre-process a dataset containing hydrodynamic measurements and five river bed morphology classes according to {cite:t}`montgomery_channel-reach_1997`.
+ Use Variance Inflation Index (VIF) to evaluate Multicollinearity among features.
+ Train and test a Support Vector Machine (SVM) model to classify river bed morphology using a Radial Basis Function (RBF) kernel.

{admonition}
:class: attention

Read and understand data handling with {ref}`numpy`, {ref}`pandas`, {ref}`matplotlib`, and {ref}`functions <chpt-functions>`. We assume that you are familiar with basic statistic parameters, such as Gaussian probability distributions, their mean, and their standard deviation values.

Enable interactive reading and executing code blocks with and find morphology-predictor-svm.ipynb Alternatively, install {ref}Python <install-python> and {ref}JupyterLab <jupyter> locally and download this Jupyter notebook.

Theory¶

In the field of restoration science, the identification of river morphology provides essential insights, for instance, to guide target terraforming aiming at re-instating a near-census natural state of a fluvial landscape. In addition, morphological pattern may serve as a predictor for estimating fluvial {term}sediment transport <Sediment transport> and vice versa. Thus, there is a bidirectional relationship between fluvial hydraulics, sediment transport, and morphological pattern, which represent physical habitat for aquatic species and may even affect energy production {cite:p}recking2016quantifying.

Yet, classifying and predicting fluvial morphodynamics is challenging because of the complexity of river ecosystems and every river being a unique environment. However,there are repeating morphological units of rivers with similar characteristics. A basic set of repetitive morphological features was introduced by {cite:t}montgomery_channel-reach_1997, and we will focus in this tutorial on the following five morphological units (click on the items to read more):

• {term}Step pool
• {term}Plane bed
• {term}Riffle pool
• {term}Braiding
• {term}Sand bed

These morphological units can be found in the above-listed order along a river, starting at the source and ending with its estuary. Thus, step pool units can typically be found in upstream (mountain) river section, while sand beds are predominantly existent in shallower river sections, close to estuaries or confluences. Note that there are many other morphological units besides, such as slackwater, swales, or bars {cite:p}wyrick_geospatial_2014. Fig. {numref}fig-bed-morphology illustrates some morphological units of near-census natural rivers.

{figure}
  :alt: riverbed morphologies
  :name: fig-bed-morphology


  a) A colluvial headwater stream (Furtschaglbach, Austria), b) a cascade stream (Torrent des Favrands, France), c) a bedrock stream (Anse St-Jean, Québec, Canada), d) a step-pool stream (Dessoubre, France), e) a plane-bed stream (Dranse, Switzerland), f) a riffle-pool stream (Le diable, Québec, Canada), g) a braided stream (Jenbach, Germany). Source: {cite:t}`schwindt_hydro-morphological_2017`

To guide restoration actions, experts often visually classify morphology units on-site, though some of them may also be distinguishable on aerial imagery. On-site (in-situ) expert assessments may also serve as ground truth for machine learning models. In this exercise, we will use hydrodynamic parameters with expert-based ground truth from the database https://bedloadweb.com. Thus, at the end of this tutorial, you may be able to reduce a necessary amount of fieldwork, just by using hydraulic parameters that can sometimes also be retrieved from gauging stations. To this end, we will train and test an SVM model that classifies morphological units as a function of data that can be collected in-situ or derived from aerial imagery.

Pre-processing¶

Load Data¶

We downloaded a dataset from https://bedloadweb.com (download the bedload_dataset sample):

• Every sample is derived from the cross-sections of 146 different rivers around the globe.
• The dataset is embraces hydraulic and sediment transport (e.g., grain sizes) though in this tutorial, we will only use the following parameters to train the SVM model:
  • W: Channel width (m)
  • S: Channel slope (m)
  • Q: Discharge (m$^3$ s$^{-1}$)
  • U: Bulk flow velocity (m s$^{-1}$)
  • H: Water depth (m)

In a machine learning context, these parameters for predicting morphological units are also called features. The first step is to load (and view) the data as a pandas.DataFrame:

Clean the Dataset¶

To avoid the SVM model being affected by statics made based on nonsense entries (e.g., Not-a-Number NaN) values or extreme outliers, for instance, resulting from typos, we will make apply some cleaning methods.

Remove NaN Values¶

The first cleaning step is to remove NaNs that may lead the SVM model later to making bad conclusions.

Search Outliers¶

Outliers may occur in extreme environment or under extreme conditions (e.g., intense precipitation), but some of them also stem from equipment failures that will affect the performance of the SVM model. Thus, recognizing and removing such device-related outliers is an essential step.

Visualizing the plotted features in a scatter plot can be helpful to spot outliers candidates. The following code block creates scatter plots of every feature, marks the expected value, and the distance of 2 standard deviations (assuming a Gaussian-like distribution) from it.

Remove Outliers¶

Outliers can be removed with two different methods:

1. Manual (i.e., expert judgment): identify outliers of the features and remove them by setting limits of their interval.
2. Automated: assume a probability distribution and set limits based on low probability of observation occurrence.

{admonition}
:class: warning

The automated method is more prone to deleting real, extreme observations. So pay attention in practice!

The following code blocks implements these two options (manual-expert and automated) methods. The manual outlier removal is hard coded because the limits should be given by experts.

The automated approach (below else statement) assumes the data having a Gaussian distribution, where points with a so-called z-score greater than 4 (four) are removed. Data points with a z-score less or equal to four have a 99.9%-chance of being inside an interval of four standard deviations away from the expected value. Thus, a sample that has a 0.01%-chance of occurrence is considered automatically an outlier and removed from the dataset. Thus, in this case, the z-score represents the standard deviation, and we use the zscore function of the scipy library.

Verify Class Proportions¶

A machine learning model requires the ground truth data to comply with some quality standard. Those involve (among others):

• Every combination of models needs a different minimum number of samples to have statistical significance. It is not a trivial task to identify the minimum number of samples that yield a trustworthy model. However, to represent a class properly (whatever that means), the so-called frequency distribution of every feature must be similar to the population frequency distribution.
• Having more than 30 samples when dealing with normally distributed variables is a good-practice rule of thumb.
• Ideally, the classes used for a classification problem should be represented by a similar quantity of samples. An unbalanced dataset may produce a model that gives preference to classifying correctly the class with a larger number of samples. This is a consequence of using accuracy as a performance metric. Accuracy is the ratio of correct and the total number of predictions ($accuracy = n_{correct\_prediction}/ n_{of\_samples}$). By preferencing the total number of correct predictions over a balanced one, accuracy masks the importance of getting the under-represented classes correctly.
• Solutions for balancing unbalanced samples embrace:
  • Collecting or generating new data (e.g., Monte Carlo methods). If available, precision is preferable for measuring model peformance (instead of accuracy).
  • Weighting class features, which means to assign weights to increase the importance of under-represented classes correctly. In classification problems, testing different combinations of weight yield an ROC curve, which shows the effect of varying weights on model performance.

In this tutorial, the samples can be considered balanced, except for the sand bed class (see output below). Weighting and ROC curves are not considered here. If sensitivity and specification of the sand bed class were acceptable when classifying test data, the model could also be used for sand bed classification without the need for using balancing techniques.

Derive New Predictors¶

Sometimes it is possible to compute new meaningful features from the existing ones. For instance, known predictors for morphology are:

• the product of (energy) slope $S$ (-), bulk flow velocity $U$ (m/s), and water depth $H$ (m), which basically represents {term}shear stress <Dimensionless bed shear stress>
• the ratio of discharge $Q$ (m$^3$/s) and width $W$ (m) (i.e., Q/W)

{warning}
Deriving new features from existing ones through simple operations may introduce spurious correlation, which results in unreliable models.

Verify Multicolinearity¶

Correlated features (i.e., statistically dependent measurement parameters) may result in unreliable models. For instance, consider a scientist who wants to use machine learning (ML) models to predict the number of tourists that get sun-burned in Brazil during their vacation. She wants to use solar radiation intensity and the number of icecream scoops sold as features. The scientist may find that both features are relevant to predict the number of sun-burned tourist, but the reason for this finding is that the number of icecream scoops sold, and sun intensity are correlated. Thus, the scientist has no way to know what is the real importance of each feature in predicting sun-burned tourists because the features grow and decrease simultaneously. By the way, icecream consumption and sun burns have no causality in common, but that is another subject.

To this end, it is important to investigate the correlation among features and consider eliminating the correlated features to produce a more reliable model. One way to visually spot colinearity is plotting the combination of two variables in a scatter plot. For this purpose, the below figure aids in spotting the following aspects of the bedload_dataset:

• There is a strong linear correlation between the $H$ and $Q/W$ features.
• There is no linear correlation between $W$ and $U$.
• The flow velocity $U$ distribuiton is approximately normal (see main diagonal).

In addition, a heatmap is another efficient method to visualize linear correlation among features. Also computing linear correlation enables quantitative analysis.

Variance Inflation Factor (VIF)¶

A straightforward method to evaluate multicollinerity is computing the Variance Inflation Factor (VIF) for every feature, which quantifies the correlation among features.

The first step to compute the VIF is treating a feature as a dependent variable and fit a linear regression model using the other features as independent variable. Second, store the corresponding mean square errors $MSE$ yielded by the linear model.

Finally, VIF can be computed for each feature by the following equation:

• $VIF_{feature(n)}$ is the variance inflation factor of the $n_{th}$ feature, and
• $MSE^2_{feature(n)}$ is the mean square error resulting from the linear regression of the $n_{th}$ feature.

Note that the smaller the MSE (better linear regression model), the greater is the inflation. In other words, large VIF values indicate high correlation of the feature under consideration with at least one of the other features. Still, it is not possible to say based on VIF only, which other feature is causing VIF to inflate.

As rule of thumb, a VIF greater than 10 (ten) indicates an unsuitable set of features to produce a reliable ML model. Ideally, the VIF should remain smaller than 5 (five), which is considered moderate correlation {cite:p}franke2010multicollinearity. However, these reference values are not a generally and always applicable rule. VIF equal to or greater than 10 can still yield a good model.

Further analysis to study the impact of a high VIF on correlated features can be carried out by removing each feature at a time and investigate how the feature removal affects the variation of the standard error and p-value of the model.

In this tutorial, 2 features have a VIF $> $10 (notably, $H$ and $Q/W$), as the following code block shows.

Remove Features to Reduce VIF¶

With the goal of reducing the feature VIFs to less than 5, we use a trial-and-error approach in this tutorial. Modifying the below code block shows that removing $H$ and $Q$ is enough to reduce the VIFs of the remaining features to less than 5. In particular, these two removed variables also go into above-introduced features $Q/W$ and $SUH$. Thus, although $H$ and $Q$ are removed, they still "contribute information" to the final model.

Combining correlated features to procduce a new one, and later removing them, is a common approach to reduce VIFs. Since we yielded good VIF values by removing some of the features, no further analysis is necessary.

{admonition}
Alternatively, a [Principal Component Analysis (PCA)](https://www.youtube.com/watch?v=FgakZw6K1QQ) is an effective method for reducing dimensionality (i.e., the number of dataset features going into a model) to increase the reliability of a model, and at the same time, minimize information loss.

Finally, the features that are used to train and test the SVM model are shown by the output of the code block below.

Split Test and Training Datasets¶

The last preparation step for building the SVM model is to derive a training and a testing dataset from the entire bedload_dataset. Choosing a share for data to compose a training and a testing dataset is commonly a heuristic choice. Here, we use 33% of the data for testing a final hypothesis later. Note that inside the split datasets, the index of the samples (rows of the tabular data) have their corresponding class inside the target datasets.

Build the SVM Model¶

Working Principle¶

A Support Vector Machine (SVM) is a learning algorithm that searches for a so-called hyperplane that optimally separates feature classes. To find an optimal separation, the algorithm tries to maximize the distance (or margin) from the hyperplane to the nearest data points. The data points near to the limits of the margin are also denominated support vectors because they work as a reference to draw the hyperplane (see Fig. {numref}svm-margin-morphology below).

{figure}
:alt: datascience
:name: svm-margin-morphology


Illustration of a hyperplane and margins of an SVM. Source: Ricardo Barros

For instance, assume that two of the above classes are linearly separable. We can define the regions that contain each class by means of the following hyperplane:

$ \vec{w} \vec{X} + b > 1 $ --> Class Red, and $ \vec{w} \vec{X} + b < 1 $ --> Class Blue

where

• $\vec{w}$ is the vector of weights,
• $\vec{X}$ is the matrix of features,
• $b$ is the independent variable (also called bias), and
• the value of $1$ is chosen as a representative simplification, since $b$ can basically be any value.

Finally (after some steps that are not discussed here), we can compute an optimal hyperplane by minimizing the following expression:

$min (0.5 \vec{w^t} \vec{w})$

where $\vec{w^t}$ is the transposed weight vector.

In the real world, however, it is common to find problems where classes are not linearly separable. In such situations, so-called kernels can be applied. A kernel is a function that transforms $\vec{X}$ to a matrix with higher dimensions by combining the available features to create new features (i.e., new dimensions).

Here, we apply the Radius Base Function (RBF) kernel. The RBF is a popular kernel because it transforms the $\vec{X}$ into an infinite dimension matrix without the risk of overfitting the model.

{admonition}
:class: hint
+ For more details on the _SVM_  calculation procedure, visit [analyticsvidhya.com](https://www.analyticsvidhya.com/blog/2021/10/support-vector-machinessvm-a-complete-guide-for-beginners/).
+ For a complete lecture on _SVM_ (provided by the [MIT](https://mit.edu)), visit their [YouTube channel](https://www.youtube.com/watch?v=eHsErlPJWUU).

k-fold Cross Validation¶

SVMs composed by an RBF kernel can be validated by tuning two parameters:

1. C: determines how much are the margins are permissive. It differentiates between a hard margin (i.e., a margin that enables little or no penetration of data points) and a soft margin (i.e., a margin that enables penetration of data points inside the margin). A higher C corresponds to a harder the margin.
2. gamma: determines the region of similarity among samples. It is the only variable of the RBF kernel. The higher gamma is, the broader is the influence of a labeled sample on classifying a new sample.

Hard margins can also be interpreted as the sensitivity of the hyperplane when fitting it with the data. The harder the margin, the more the hyperplane twists to not allow any data point to penetrate the margin.

In this tutorial, we define C based on a try-and-error approach and a gamma stemming from a five-fold cross-validation. The following code block builds the SVM using the sklearn library, and prints the validation and training error curves resulting from varying gamma.

{admonition}
:class: hint
For a _C_ and _gamma_ search code example, have a look into the [scikit-learn docs](https://scikit-learn.org/stable/auto_examples/svm/plot_rbf_parameters.html#:~:text=Intuitively%2C%20the%20gamma%20parameter%20defines,the%20model%20as%20support%20vectors.).

Select Optimum Parameters¶

The decision on optimum parameters for the SVM depends on your opinion. For the purpose of selecting optimum parameters in this tutorial, the below code block quantifies the above-shown plot (in particular, the blue curve) in tabular form and saves the gamma that yields the smallest error (gamma_opt).

Train Optimum Hypothesis SVM Model¶

After defining a gamma yielding the smallest validation error (gamma_opt), we can train an optimal hypothesis (h_opt) model on the entire training dataset. In addition, we can evaluate the accuracy of the model by using it to classify the testing dataset. The below code block features the training procedure, fitting of the SVM (i.e., an optimum hypothesis), and prints the score of the optimum hypothesis.

Performance Evaluation (Confusion Matrix)¶

Finally, the following code block generates the so-called confusion matrix to evaluate the optimum hypothesis model performance when classifying the samples with distinct morphologies. The confusion matrix uses the total number of samples and visualizes their true and predicted morphological units.

The confusion matrix of an exact model only has diagonal entries larger than zero, while all others are zero. Thus, the above-shown confusion matrix shows that our SVM model performs very well in predicting step-pool and riffle-pool morphologies, but unreliably predicts sand bed morphologies because of not many observations. However, the model also made wrong predictions of step-pool and riffle pool features. For instance, according to the confusion matrix, the model predicted 6 times step-pool, where in reality (true label), plane bed morphological units prevailed, and it wrongly predicted one braiding observation to be a step-pool observation. The sand bed predictions are unreliable because there are only 25 true-label sand bed observations (opposed to 66 step-pool observations). Thus, the sand bed observations are unbalanced and a solution to removing the unreliability could be to sample more sand bed morphological units (as above explained).