Pre-processing¶

Load Data¶

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

• The samples stem from cross-sections of 146 different rivers around the globe.

• The dataset embraces hydraulic and sediment transport parameters (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/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 statistics based on nonsense entries (e.g., Not-a-Number NaN values) or extreme outliers, for instance, resulting from typos, we will apply some cleaning methods.

Remove NaN Values¶

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

Search Outliers¶

Outliers may occur in extreme environments 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 features in scatter plots can be helpful to spot outlier candidates. The following code block creates a scatter plot of every feature, and marks its mean value (horizontal line) and the range of ±2 standard deviations around it (vertical line, assuming a Gaussian-like distribution).

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 a low probability of occurrence.

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

The automated approach (below the else statement) assumes that the data follow a Gaussian distribution, and removes data points with a so-called z-score of more than 4 (four) in absolute terms. The z-score expresses by how many standard deviations a value deviates from the mean. For a Gaussian distribution, approximately 99.994 % of all values lie within four standard deviations from the mean. Thus, a sample with an occurrence probability of less than 0.01 % is automatically considered an outlier and removed from the dataset. To this end, 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 standards. Those involve (among others):

• Every model 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 number of samples. An unbalanced dataset may produce a model that gives preference to correctly classifying 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 predictions and the total number of predictions ($accuracy = n_{correct\_predictions}/ n_{samples}$). By favoring the total number of correct predictions over a balanced one, accuracy masks the importance of correctly predicting the under-represented classes. Therefore, class-sensitive metrics, such as precision and recall, are preferable for measuring the performance of models trained on unbalanced data.

• Solutions for balancing unbalanced samples embrace:

  • Collecting new data or generating synthetic data (e.g., by oversampling the under-represented classes or with Monte Carlo methods).

  • Weighting classes, which means to assign higher weights to under-represented classes to increase their importance in the training process. The effect of weights or decision thresholds on the trade-off between true and false positive rates can be evaluated with an ROC curve.

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 the sensitivity and specificity 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 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 is proportional to the specific stream power $\omega = \rho g U H S$ (W/m$^2$), where $\rho$ is the density of water and $g$ is gravitational acceleration; similarly, the product of $H$ and $S$ is proportional to the bed shear stress $\tau = \rho g H S$

• the ratio of discharge $Q$ (m$^3$/s) and width $W$ (m), which is the unit discharge $q = Q/W$ (m$^2$/s)

Verify Multicollinearity¶

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 sunburned in Brazil during their vacation. She wants to use solar radiation intensity and the number of ice cream scoops sold as features. The scientist may find that both features are relevant to predict the number of sunburned tourists, but the reason for this finding is that the number of ice cream scoops sold and solar radiation intensity are correlated. Thus, the scientist has no way to know the real importance of each feature in predicting sunburned tourists because the features grow and decrease simultaneously. By the way, ice cream consumption does not cause sunburns (nor vice versa); both are driven by a common cause (sunshine), but that is another subject.

To this end, it is important to investigate the correlation among features and consider eliminating correlated features to produce a more reliable model. One way to visually spot collinearity 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, which is not surprising because $Q/W = U \cdot H$.

• There is only a weak linear correlation between $W$ and $U$ (correlation coefficient of approximately 0.37, see the heatmap below).

• The flow velocity $U$ distribution is roughly bell-shaped (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 multicollinearity is computing the Variance Inflation Factor (VIF) for every feature, which quantifies how well a feature can be explained by the other features.

The first step to compute the VIF is treating a feature as a dependent variable and fitting a linear regression model using the other features as independent variables. Second, store the corresponding coefficient of determination $R^2$ yielded by the linear model.

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

• $VIF_{feature(n)}$ is the variance inflation factor of the $n$-th feature, and

• $R^2_{feature(n)}$ is the coefficient of determination resulting from the linear regression of the $n$-th feature on all other features.

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

As a 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 Franke, 2010. However, these reference values are not a generally and always applicable rule. A 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 one feature at a time and investigating 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 the 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 produce 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.

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 of the data to compose a training and a testing dataset is commonly a heuristic choice. Here, we use 30 % of the data for testing a final hypothesis later. Note that inside the split datasets, the indices of the samples (rows of the tabular data) correspond to the indices of their classes in 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 located on (or inside) the limits of the margin are also denominated support vectors because they work as a reference to draw the hyperplane (see Figure 2 below).

For instance, assume that two of the above classes are linearly separable. The separating hyperplane is defined by $\vec{w} \cdot \vec{x} + b = 0$, and the regions that contain each class can be defined as:

$\vec{w} \cdot \vec{x} + b \geq +1$ --> Class Red, and $\vec{w} \cdot \vec{x} + b \leq -1$ --> Class Blue

where

• $\vec{w}$ is the vector of weights (normal to the hyperplane),

• $\vec{x}$ is the feature vector of a sample,

• $b$ is the intercept (also called bias), and

• the value of 1 is a scaling convention, since $\vec{w}$ and $b$ can be multiplied by any positive factor without changing the hyperplane.

With this convention, the margin width is $2 / \lVert \vec{w} \rVert$. Thus, maximizing the margin (after some steps that are not discussed here) corresponds to minimizing the following expression:

$\min \left(0.5 \, \vec{w}^T \vec{w}\right)$

where $\vec{w}^T$ is the transposed weight vector, subject to the constraint that all training samples lie on the correct side of the margin.

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 implicitly maps the features $\vec{x}$ into a space with higher dimensions (i.e., new features that are combinations of the available features), where the classes may become linearly separable. The so-called kernel trick makes this possible without explicitly computing the new features.

Here, we apply the Radial Basis Function (RBF) kernel. The RBF is a popular kernel because it implicitly maps $\vec{x}$ into an infinite-dimensional feature space, which makes it very flexible. However, this flexibility also means that an RBF-SVM can overfit the training data when its parameters are not carefully tuned (see the next section).

k-Fold Cross Validation¶

SVMs with an RBF kernel can be tuned with two parameters:

1. C: determines how permissive the margins are. 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 margin.

2. gamma: determines the region of similarity among samples. It is the only parameter of the RBF kernel. The higher gamma is, the narrower is the influence of a labeled sample on classifying a new sample (i.e., only very close samples are considered similar).

Hard margins can also be interpreted as the sensitivity of the hyperplane when fitting it to the data. The harder the margin, the more the hyperplane twists to not allow any data point to penetrate the margin. Thus, both a high C and a high gamma increase the risk of overfitting.

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

Select Optimum Parameters¶

The decision on optimum parameters for the SVM depends on your judgment. 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 error of the optimum hypothesis on the testing dataset.

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. The above-shown confusion matrix shows that our SVM model performs well in predicting step-pool (62 out of 66 true-label observations) and riffle-pool (62 out of 69) morphologies. Despite the smaller number of observations, the model also predicts most sand bed observations correctly (22 out of 25). However, these 25 true-label sand bed observations (opposed to 66 step-pool or 86 plane bed observations) are few, which makes the sand bed statistics less robust. Thus, the sand bed observations are unbalanced, and a solution for increasing the reliability could be to sample more sand bed morphological units (as explained above).

The weakest performance concerns plane bed morphologies: only 54 out of 86 true-label plane bed observations were predicted correctly. In particular, the model predicted 18 times riffle-pool and 6 times step-pool, where in reality (true label), plane bed morphological units prevailed. In addition, the model wrongly predicted one braiding observation to be a step-pool observation. These confusions are plausible from a physical point of view because plane beds represent a transitional morphology between step-pool and riffle-pool channels Montgomery & Buffington, 1997.