[1]Sergey Kucheryavskiy

[1]Department of Chemistry and Bioscience, Aalborg University, Niels Bohrs vej
8, Esbjerg, 6700, Denmark

2]CSort, LLC., Germana Titova st. 7, Barnaul, 656023, Russia

# Collinear datasets augmentation using Procrustes validation sets

<EMAIL_ADDRESS>Sergei Zhilin<EMAIL_ADDRESS>*  [

###### Abstract

In this paper, we propose a new method for the augmentation of numeric and
mixed datasets. The method generates additional data points by utilizing
cross-validation resampling and latent variable modeling. It is particularly
efficient for datasets with moderate to high degrees of collinearity, as it
directly utilizes this property for generation. The method is simple, fast,
and has very few parameters, which, as shown in the paper, do not require
specific tuning. It has been tested on several real datasets; here, we report
detailed results for two cases, prediction of protein in minced meat based on
near infrared spectra (fully numeric data with high degree of collinearity)
and discrimination of patients referred for coronary angiography (mixed data,
with both numeric and categorical variables, and moderate collinearity). In
both cases, artificial neural networks were employed for developing the
regression and the discrimination models. The results show a clear improvement
in the performance of the models; thus for the prediction of meat protein,
fitting the model to the augmented data resulted in a reduction in the root
mean squared error computed for the independent test set by 1.5 to 3 times.

###### keywords:

data augmentation, artificial neural networks, Procrustes cross-validation,
latent variables, collinearity

##  1 Introduction

Modern machine learning methods that rely on high complexity models, such as
artificial neural networks (ANN), require a large amount of data to train and
optimize the models. Insufficient training data often lead to overfitting
problems, as the number of model hyperparameters to tune is much larger than
the number of degrees of freedom in the dataset.

Another common issue in this case is the lack of reproducibility because the
ANN training procedure is not deterministic, given the random selection of
initial model parameters and the stochastic nature of their optimization.
Consequently, it never leads to a model with the same parameters and
performance, as different training trials can result in different models. This
variability becomes large if the training set is too small.

This problem is particularly urgent in the case of fitting the experimental
data, as it is often expensive and time-consuming to run many experimental
trials, making it simply impossible to collect thousands of measurements
needed for proper training and optimization. There can also be other
obstacles, such as paperwork related to permissions in medical research.

One way to overcome the problem of insufficient training data is to
artificially augment it by either simulating new data points or making small
modifications to existing ones. This technique is often referred to as “data
augmentation”. Data augmentation has proved to be particularly efficient in
image analysis and classification, with a large body of research reporting
both versatile augmentation methods [1] [2], [3] and methods that are
particularly effective for specific cases [4] [5]. Augmentation methods for
time series data are also relatively well developed [6].

However, there is a lack of efficient methods that can provide decent data
augmentation for numeric datasets with a moderate to high degree of
collinearity. Such datasets are widespread in experimental research, including
various types of spectroscopic data, results of genome sequencing (e.g., 16S
RNA), and many others. Many tabulated datasets also exhibit internal
structures where variables are mutually correlated. Currently available
methods for augmentation of such data mostly rely on adding various forms of
noise [7] to the existing measurements, which is not always sufficient. There
are also promising methods that utilize variational autoencoders by random
sampling from their latent variable space [8], or methods based on generative
adversarial networks [4]. The downsides are that both approaches require
building and tuning a specific neural network model for the data augmentation
and hence need a thorough and resource demanding optimization process and a
relatively large initial training set.

In this paper, we propose a simple, fast, versatile, yet efficient method for
augmenting numeric and mixed collinear datasets. The method is based on an
approach that was initially developed for other purposes, specifically for
generating validation sets, and hence is known as Procrustes cross-validation
[9] [10]. However, as demonstrated in this paper, it effectively addresses the
data augmentation problem, resulting in models with significantly improved
prediction or classification performance.

Our method directly leverages collinearity in the generation procedure. It
fits the training data with a set of latent variables and then employs cross-
validation resampling to measure variations in the orientation of the
variables. This variation is then introduced to the training set as sampling
error, resulting in a new set of data points.

Two fitting models can be employed — singular value decomposition (SVD) and
partial least squares (PLS) decomposition. The choice of the fitting model
allows the user to prioritize a part of covariance structure, which will be
used for generation of the new data.

Both fitting models have two parameters — the number of latent variables and
the number of segments used for cross-validation resampling. The experiments
show though that the parameters do not require specific tuning. Any number of
latent variables large enough to capture the systematic variation of the
training set values serve equally well. As well as any number of segments
starting from three.

The proposed method is versatile and can be applied to both fully numeric data
as well as to tabulated data where one or several variables are qualitative.
This opens another perspective, namely data mocking, which can be useful,
e.g., for testing of high loaded software systems, although we do not consider
this aspect here.

The paper describes the theoretical foundations of the method and illustrates
its practical application and performance based on two datasets of different
nature. It provides comprehensive details on how the method can be effectively
applied to diverse datasets in real-world scenarios.

We have implemented the method in several programming languages, including
Python, R, MATLAB, and JavaScript, and all implementations are freely
available in the GitHub repository (https://github.com/svkucheryavski/pcv).
Additionally, we provide an online version where one can generate new data
points directly in a browser (https://mda.tools/pcv).

##  2 Methods

This chapter presents the general idea behind the Procrustes cross-validation
approach, provides the necessary mathematical background and describes two
implementations of the approach in detail.

Let $\\{\mathbf{X},\mathbf{Y}\\}$ be two matrices of size $I\times J$ and
$I\times M$ correspondingly, representing the original training set to be
augmented. Columns of matrix $\mathbf{X}$ are predictors — variables measured
or observed directly, for example, light absorbance at different wavelengths
in the case of spectra, abundance of operational taxonomic units (OTU) in the
case of 16S RNA sequencing, results of clinical tests for patients etc. If the
original dataset has qualitative predictors, they should be replaced by
corresponding dummy variables, which are combined with quantitative predictors
to form $\mathbf{X}$.

Columns of matrix $\mathbf{Y}$ are responses — variables whose values should
be predicted. There are situations when matrix $\mathbf{Y}$ is obsolete, for
example, in the case of developing of one class classifier, only matrix with
predictors is needed.

To augment the data a specific model should be developed. It aims at capturing
the variance-covariance structure of $\mathbf{X}$, or its part, directly
related to the variance of values in $\mathbf{Y}$. This model will be referred
to as the _PV-model_ and denoted as $\mathcal{M}$.

The selection of the methods for creating the PV-model depends on the
objectives. Thus, if data augmentation is needed for developing a one class
classifier, when the matrix with responses, $\mathbf{Y}$, is obsolete, we
suggest using singular value decomposition (SVD), which is efficient in
capturing the whole variance-covariance structure of ${\mathbf{X}}$.

If the augmented data are used for developing of regression models, partial
least squares (PLS) decomposition, which accounts for the variance-covariance
structure of $\mathbf{X}^{\textrm{T}}\mathbf{Y}$ will be the favorable choice.

In the case of augmenting data for discrimination models, there can be two
solutions. One solution is to consider each class as an independent population
and use SVD decomposition to capture the variance-covariance structure of
$\mathbf{X}$ within each class. The second solution will be to create a matrix
with dummy predictors, $\mathbf{Y}$, where each column corresponds to a
particular class, and employ PLS decomposition. More details about using SVD
and PLS decompositions are provided in the next two sections.

A general definition of the approach can be formulated as follows.

Let us create a _global_ PV-model $\mathcal{M}$ by fitting all data points
from the training set $\\{\mathbf{X},\mathbf{Y}\\}$. If this model is then
applied to a new dataset, $\\{\mathbf{X}^{*},\mathbf{Y}^{*}\\}$, it will
result in a set of outcomes $\mathbf{R}^{*}$:

${\mathcal{M}(\mathbf{X}^{*},\mathbf{Y}^{*})\rightarrow\mathbf{R}^{*}}$ (1)

The nature of the outcomes depends on the method used for developing of
$\mathcal{M}$, and can include, for example, predicted response values in the
case of PLS decomposition and scores and distances in the case of SVD
decomposition.

Let us consider cross-validation resampling with $K$ segments. For any segment
$k=1,...,K$, the dataset $\\{\mathbf{X},\mathbf{Y}\\}$ is split into a local
validation set $\\{\mathbf{X}_{k},\mathbf{Y}_{k}\\}$ with $I_{k}=I/K$ rows and
a local training set, $\\{\tilde{\mathbf{X}}_{k},\tilde{\mathbf{Y}}_{k}\\}$
with $I-I_{k}$ rows. The local training set is used to create a local PV-model
$\mathcal{M}_{k}$, which is then applied to the local validation set to obtain
the outcomes $\mathbf{R}_{k}$:

${\mathcal{M}_{k}(\mathbf{X}_{k},\mathbf{Y}_{k})\rightarrow\mathbf{R}_{k}}$
(2)

Now let us consider another matrix with predictors, $\mathbf{X}_{pv}$, which
we will call a _PV-set_. It has the same size as $\mathbf{X}$ and the rows of
the two matrices are related.

This matrix should be created to obtain the outcomes, computed by applying the
global model $\mathcal{M}$ to the same $I_{k}$ rows of the PV-set,
$\mathbf{X}_{pv_{k}}$, to be as close as possible to the outcomes computed by
applying the local model $\mathcal{M}_{k}$ to the local validation set:

${\mathcal{M}(\mathbf{X}_{pv_{k}},\mathbf{Y}_{k})\rightarrow\mathbf{R}_{pv_{k}}}$
(3)

${\mathbf{R}_{pv_{k}}\approx\mathbf{R}_{k}}$ (4)

for any $k=1,...,K$.

This criterion is called a _Procrustean rule_. The main idea of the proposed
approach is that the PV-set, $\mathbf{X}_{pv}$, created using the general
procedure described above, on the one hand will have a variance-covariance
structure (full or partial), similar to the variance-covariance structure of
the predictors from the original training set, $\mathbf{X}$. On the other
hand, it will have unique values, as if it represents another sample taken
from the same population as $\mathbf{X}$.

Using cross-validation with random splits enables the generation of a very
large number of the PV-sets, which will hold the properties described above,
but will be different from each other. Combining the generated PV-sets with
the original training data will result in the augmented training set we are
aiming at.

The next two sections describe the details of PV-set generation using SVD and
PLS decompositions. More information about the approach can be found in [10].

###  2.1 Generation of PV-sets based on Singular Value Decomposition

In the case of truncated SVD with $A$ latent variables, the PV-model
$\mathcal{M}$ consists of:

  * •

number of latent variables (singular vectors), $A$

  * •

$1\times J$ row-vector for mean centering, $\mathbf{m}$

  * •

$1\times J$ row-vector for standardization (if required), $\mathbf{s}$

  * •

$J\times A$ matrix with right singular vectors, $\mathbf{V}$

  * •

$1\times A$ vector with corresponding singular values, $\mathbf{\sigma}$.

The SVD decomposition of $\mathbf{X}$ is defined as:

${\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\textrm{T}}+\mathbf{E}}$ (5)

Here, $\mathbf{U}$ is an $I\times A$ matrix with left singular vectors and
$\mathbf{\Sigma}$ is an $A\times A$ diagonal matrix with singular values
$\mathbf{\sigma}=\\{\sigma_{1},...,\sigma_{A}\\}$. Element $\mathbf{E}$ is an
$I\times J$ matrix with residuals, which can be used to estimate the lack of
fit both for individual data points and for the whole dataset. For the sake of
simplicity we assume that columns of $\mathbf{X}$ are already mean centered
and standardized.

The right singular vectors $\mathbf{V}$ act as the orthonormal basis of the
$A$-dimensional latent variable space located inside the original
$J$-dimensional variable space ($A\leq J$). The product of the left singular
vectors and the singular values, $\mathbf{T}=\mathbf{U}\mathbf{\Sigma}$, forms
_scores_ — coordinates of data points from $\mathbf{X}$ being projected to the
latent variable space.

The expression can be also be written as:

${\mathbf{X}=\mathbf{TV}^{\textrm{T}}+\mathbf{E}=\hat{\mathbf{X}}+\mathbf{E}}$
(6)

where $\hat{\mathbf{X}}$ is a part of $\mathbf{X}$ that is captured (or
explained) by the SVD model. In the case of full decomposition, when
$A=rank(\mathbf{X})$, $\mathbf{X}=\hat{\mathbf{X}}$.

Any new data point, $\mathbf{x}=\\{x_{1},...,x_{J}\\}$ can be projected to the
model $\mathcal{M}$ by mean centering and (if required) standardizing its
values using $\mathbf{m}$ and $\mathbf{s}$, and then projecting the resulting
values to the latent variable space defined by the columns of $\mathbf{V}$:

${\mathbf{t}=\mathbf{x}\mathbf{V}}$ (7)

The score vector $\mathbf{t}$ can be normalized to obtain the corresponding
left singular vector, $\mathbf{u}$:

${\mathbf{u}=\mathbf{t}\mathbf{\Sigma}^{-1}=\\{t_{1}/\sigma_{1},...,t_{A}/\sigma_{A}\\}}$
(8)

The explained part of $\mathbf{x}$ can be computed as:

${\hat{\mathbf{x}}=\mathbf{t}\mathbf{V}^{T}=\mathbf{u}\mathbf{\Sigma}\mathbf{V}^{T}}$
(9)

The residual part is:

${\mathbf{e}=\mathbf{x}-\hat{\mathbf{x}}=\mathbf{x}-\mathbf{u}\mathbf{\Sigma}\mathbf{V}=\mathbf{x}(\mathbf{I}-\mathbf{VV}^{\textrm{T}})}$
(10)

One can define a relationship between the data point $\mathbf{x}$ and the PV-
model $\mathcal{M}$ using two distances. A squared Euclidean distance, $q$,
between $\mathbf{x}$ and $\hat{\mathbf{x}}$ in the original variable space:

${q=\sum_{j=1}^{J}(x_{j}-\hat{x}_{j})^{2}=\sum_{j=1}^{J}e_{j}^{2}}$ (11)

And a squared Mahalanobis distance, $h$, between the projection of the point
in the latent variable space and the origin:

${h=\sum_{a=1}^{A}u_{a}^{2}}$ (12)

The first distance, $q$, is a measure of lack of fit for the point, while the
second distance, $h$, is a measure of extremeness of the point (as the
majority of the data points will be located around the origin due to mean
centering).

If we apply cross-validation resampling to the SVD decomposition, it will
result in a set of $K$ local PV-models, $\mathcal{M}_{k}$, created using the
local training sets, $\tilde{\mathbf{X}}_{k}$ ($k=1,...K$).

The main difference between the global PV-model, $\mathcal{M}$, and each local
PV-model, $\mathcal{M}_{k}$, is the orientation of the right singular vectors,
columns of $\mathbf{V}$ and $\mathbf{V}_{k}$, in the variable space. If we
take two singular vectors, the $a$-th vector from the global model,
$\mathbf{v}_{a}$, and the $a$-th vector from the $k$-th local model,
$\mathbf{v}_{ak}$, then the angle between the vectors in the original variable
space can be considered as a measure of sampling error, emulated by the cross-
validation resampling.

If we introduce this error to the local validation set $\mathbf{X}_{k}$, we
will create a new set of measurements, $\mathbf{X}_{pv_{k}}$. In the case of
full rank decomposition this can be done as:

${\mathbf{X}_{pv_{k}}=\mathbf{X}_{k}\mathbf{V}_{k}\mathbf{V}^{\textrm{T}}}$
(13)

because the dot product $\mathbf{V}_{k}\mathbf{V}^{\textrm{T}}$ is a
rotational matrix between the two orthonormal spaces, defined by
$\mathbf{V}_{k}$ and $\mathbf{V}$.

It can be shown easily that in this case the following Procrustean rules will
be hold:

${\mathbf{q}_{pv_{k}}=\mathbf{q}_{k}}$ (14)

${\mathbf{h}_{pv_{k}}=\mathbf{h}_{k}}$ (15)

where $\mathbf{q}_{pv_{k}}$ and $\mathbf{h}_{pv_{k}}$ are vectors with the two
distances, $q$ and $h$, computed for each data point from the
$\mathbf{X}_{pv_{k}}$ being projected to the global model $\mathcal{M}$, while
$\mathbf{q_{k}}$ and $\mathbf{h_{k}}$ are vectors with the two distances
computed for each data point from the $\mathbf{X}_{k}$ being projected to the
local model $\mathcal{M}_{k}$.

Hence we meet the general Procrustes rule requirement defined by Equation 4
with the outcomes $\mathbf{R}$ consisting of the two distances. Moreover, this
rule holds true for the distances computed using any number of latent
variables $a=1,...,A$.

The generation of the PV-set in the case when $A<rank(\mathbf{X})$ is similar
but requires an additional step to hold the rule defined by Equation 14,
further details can be found in [10].

By repeating this procedure for all segments, $k=1,...,K$, we form the
complete PV-set, $\mathbf{X}_{pv}$, which will have the same number of data
points as the original training set $\mathbf{X}$ and will hold the same
variance-covariance structure; however it will have its own sampling error
estimated by the cross-validation resampling.

Repeating the generation using random cross-validation splits provides a large
number of the unique PV-sets that can be used to augment the original training
set. Therefore, the SVD based data augmentation procedure has two parameters:

  * •

number of segments, $K$

  * •

number of latent variables, $A$

As it will be shown in the results section, neither parameter significantly
influences the quality of the augmented data. In the case of the number of
latent variables, $A$, any number large enough for capturing the systematic
variation in $\mathbf{X}$, will work well. The SVD based PV-set generation
does not suffer from overfitting, and optimization of $A$ is not needed.

###  2.2 Generation of PV-sets based on PLS decomposition

Partial least squares is a decomposition used for solving multiple linear
regression problem in case when predictors (and responses if they are
multiple) are collinear. The decomposition can be expressed in the following
set of equations:

${\mathbf{T}=\mathbf{X}\mathbf{W}}$ (16)

${\mathbf{X}=\mathbf{T}\mathbf{P}^{\textrm{T}}+\mathbf{E}_{x}=\hat{\mathbf{X}}+\mathbf{E}_{x}}$
(17)

${\mathbf{Y}=\mathbf{T}\mathbf{C}^{\textrm{T}}+\mathbf{E}_{y}=\hat{\mathbf{Y}}+\mathbf{E}_{y}}$
(18)

Here, matrix $\mathbf{W}$ is a matrix of PLS-weights, which are computed based
on covariance matrix $\mathbf{X}^{\textrm{T}}\mathbf{Y}$ and the consequent
deflation procedure using the SIMPLS algorithm [11]. Every column of
$\mathbf{W}$ is a unit vector, defining the orientation of latent variables in
the predictors space, similar to the columns of $\mathbf{V}$ in SVD
decomposition. Therefore, the matrix of PLS-scores, $\mathbf{T}$, contains
coordinates of the original data points being projected to the latent
variables defined by the columns of $\mathbf{W}$. The number of latent
variables, $A$, is a decomposition parameter.

In the case of one response, matrix $\mathbf{Y}$ can be replaced by a column
vector $\mathbf{y}$ containing $I$ response values. Then the Equation 18 can
be rewritten in the following compact form:

${\hat{\mathbf{y}}=\mathbf{T}\mathbf{c}^{\textrm{T}}=c_{1}\mathbf{t}_{1}+c_{2}\mathbf{t}_{2}+...+c_{A}\mathbf{t}_{A}}$
(19)

Matrix $\mathbf{C}$ in this case is represented by a column vector
$\mathbf{c}=\\{c_{1},...,c_{A}\\}$ and vectors
$\mathbf{t}_{1}$,…,$\mathbf{t}_{A}$ are columns of the score matrix
$\mathbf{T}$.

In contrast to SVD, the latent variables in PLS are oriented along directions
in the column space of $\mathbf{X}$, which give the largest covariance between
the scores (columns of matrix $\mathbf{T}$) and the values of $\mathbf{y}$.
Hence, by using PLS, we can prioritize the variance-covariance structure of
$\mathbf{X}$ directly related to $\mathbf{y}$.

The difference between the original and estimated response values can be
defined as:

${\mathbf{e}=\mathbf{y}-\hat{\mathbf{y}}=\mathbf{y}-\mathbf{T}\mathbf{c}^{\textrm{T}}}$
(20)

If we apply cross-validation resampling, then for each segment $k$
($k=1,...,K$), a local PLS-model, $\mathcal{M}_{k}$ is computed. The model
consists of, among others, the matrix with weights, $\mathbf{W}_{k}$, and the
vector with y-loadings, $\mathbf{c}_{k}$, obtained using the local training
set $\\{\tilde{\mathbf{X}}_{k},\tilde{\mathbf{y}}_{k}\\}$.

Applying the local model to the local validation set
$\\{\mathbf{X}_{k},\mathbf{y}_{k}\\}$ results in the following:

${\mathbf{T}_{k}=\mathbf{X}_{k}\mathbf{W}_{k}}$ (21)

${\hat{\mathbf{y}}_{k}=\mathbf{T}_{k}\mathbf{c}_{k}}$ (22)

${\mathbf{e}_{k}=\mathbf{y}_{k}-\hat{\mathbf{y}}_{k}=\mathbf{y}_{k}-\mathbf{T}_{k}\mathbf{c}_{k}^{\textrm{T}}}$
(23)

Now, let us assume that we have a PV-set, $\\{\mathbf{X}_{pv},\mathbf{y}\\}$.
If we take $I_{k}$ rows from the PV-set,
$\\{\mathbf{X}_{pv_{k}},\mathbf{y}_{k}\\}$, which correspond to the rows of
the local validation set, $\\{\mathbf{X}_{k},\mathbf{y}_{k}\\}$, and apply the
global model to the rows we obtain:

${\mathbf{T}_{pv_{k}}=\mathbf{X}_{pv_{k}}\mathbf{W}}$ (24)

${\hat{\mathbf{y}}_{pv_{k}}=\mathbf{T}_{pv_{k}}\mathbf{c}}$ (25)

${\mathbf{e}_{pv_{k}}=\mathbf{y}_{k}-\hat{\mathbf{y}}_{pv_{k}}=\mathbf{y}_{k}-\mathbf{T}_{pv_{k}}\mathbf{c}^{\textrm{T}}}$
(26)

We can define the Procrustean rule for the regression as:

${\mathbf{e}_{pv_{k}}=\mathbf{e}_{k}}$ (27)

Taking into account the equations above, this equation can be simplified to:

${\hat{\mathbf{y}}_{pv_{k}}=\hat{\mathbf{y}}_{k}}$ (28)

or in the case of PLS:

${\mathbf{T}_{pv_{k}}\mathbf{c}=\mathbf{T}_{k}\mathbf{c}_{k}}$ (29)

which gives us the expression for computing $\mathbf{T}_{pv_{k}}$:

$\mathbf{T}_{pv_{k}}=\mathbf{T}_{k}\mathbf{D}_{k}=\mathbf{X}_{k}\mathbf{W}_{k}\mathbf{D}_{k}$

where $\mathbf{D}_{k}$ is a diagonal matrix with scalars
$\\{c_{k1}/c_{1},...,c_{kA}/c_{A}\\}$ along the main diagonal. Finally, taking
into account Equation 17, $\mathbf{X}_{pv_{k}}$ can be computed as:

$\hat{\mathbf{X}}_{pv_{k}}=\mathbf{T}_{pv_{k}}\mathbf{P}^{\textrm{T}}=\mathbf{T}_{k}\mathbf{D}_{k}\mathbf{P}^{\textrm{T}}=\mathbf{X}_{k}\mathbf{W}_{k}\mathbf{D}_{k}\mathbf{P}^{\textrm{T}}$

In contrast to SVD, when points corresponding to $\mathbf{X}_{pv}$, are
computed mostly by rotations of the points from $\mathbf{X}$, in this case a
more advanced procedure, which includes both scaling and rotations, is
employed.

In the case of full decomposition, when $A=rank(\mathbf{X})$,
$\mathbf{X}_{pv_{k}}=\hat{\mathbf{X}}_{pv_{k}}$. When the decomposition is not
full, an additional step is needed in order to compute the matrix with
residuals, $\mathbf{E}_{x}$. This step is similar to the one used for SVD
decomposition and can be found in [10]. The procedure can also be generalized
to the case of multiple response variables, which is also described in [10].

By repeating this procedure for all $k=1,...,K$, the whole PV-set,
$\mathbf{X}_{pv}$, is computed. Similar to SVD, cross-validation with random
splits creates large number of PV-sets that can be used to augment the
original matrix with predictors, $\mathbf{X}$. Note, that the response values
for the augmented part will be identical to the training response
($\mathbf{y}_{pv}=\mathbf{y}$).

Therefore, PV-set generation in the case of PLS has the same two parameters as
the SVD-based procedure:

  * •

number of segments, $K$

  * •

number of latent variables, $A$

It must also be noted that the number of latent variables, $A$, in case of PLS
must be limited (in contrast to the SVD based algorithm). If too many latent
variables are selected, the variation in scalar values, $c_{k}/c$, can be very
large, which leads to noisy PV-sets. This happens because higher latent
variables do not have any covariance with the response variable, $y$; hence,
$c_{k}$ values vary chaotically. It is recommended to keep all latent
variables with absolute values of $c_{k}/c$ within the $[0,2]$ range.

As it is shown in the experimental part, if this limit is introduced, the
number of latent variables does not have a substantial influence on the PV-
sets quality and does not require specific optimization. Same can be concluded
about the number of segments, $K$.

##  3 Results

###  3.1 Datasets

Two datasets were selected to demonstrate the capabilities of PV-sets to be
used for data augmentation. Since both SVD and PLS decompositions model the
variance-covariance structure of $\mathbf{X}$, the proposed augmentation
approach will work best for datasets with moderate to high degree of
collinearity, which was taken into account for the selection of the datasets.

It must be noted, that optimization of the ANN architecture was outside the
scope of this work. In both examples a very simple architecture with several
layers was used.

All calculations were carried out using Python 3.10.12 supplemented with
packages _PyTorch_[12] 2.0.1, _NumPy_[13] 1.26.0 and _prcv_ 1.1.0. Statistical
analysis of the results was performed in R 4.3.1.

Python scripts used to obtain the results, presented in this chapter, are
freely available so everyone can reproduce the reported results.

####  3.1.1 Tecator

The _Tecator_ dataset consists of spectroscopic measurements of finely minced
meat samples with different moisture, protein and fat content. The
measurements were taken by Tecator Infratec Food and Feed Analyzer working by
the Near Infrared Transmission (NIT) principle. The dataset was downloaded
from the StatLib public dataset archive (http://lib.stat.cmu.edu/datasets/).

The dataset contains the spectra and the response values for 215 meat samples
in total (subsets labeled as _C_ , _M_ and _T_ in the archive), of which 170
samples are used as the training set and 45 are used as the independent test
set to assess the performance of the final optimized models. Fat content (as %
w/w) was selected as the response variable for this research.

The matrix of predictors $\mathbf{X}$ contains absorbance values measured at
100 wavelengths in the range of 850–1050 nm.

Figure 1: Tecator spectra colored by fat content in the samples (colorbar
legend maps the colors of the spectra to the fat content in %w/w).

According to the Beer-Lambert law, there is a linear relationship between the
concentration of chemical components and the absorbance of the electromagnetic
radiation. Therefore, the prediction of fat content can in theory be performed
by fitting the dataset with a multiple linear regression model. However, the
shape of the real spectra suffers from various side effects, including noise,
light scattering from an uneven surface of the samples, and many others.

Figure 1 shows the spectra from the training set in the form of line plots
colored according to the corresponding fat content using a color gradient from
blue (low content) to red (high content) colors. As one can see, there is no
clear association between the shape of the spectra and the response value.

Handling such data necessitates the careful selection and optimization of
preprocessing methods to eliminate undesirable effects and reveal the
information of interest. However, the use of nonlinear models, such as
artificial neural networks, has the capability to automatically unveil the
needed information [14], bypassing the need for extensive preprocessing steps.

Nonlinearity is not the sole factor in play. Preprocessing can be considered
as a feature extraction procedure. It is known that in networks with many
hidden layers, the initial layers act as feature extractors. The greater the
number of layers, the more intricate features can be derived from raw data.
However, it is worth noting that such models demand a substantial number of
measurements or observations to achieve satisfactory performance.

The _Tecator_ dataset will be employed for a thorough investigation of how the
use of PV-set based data augmentation can attack the problem and how the two
main parameters of the PV-set generation procedure (number of latent
variables, number of segments) as well as the number of generated PV-sets
influence the performance of the ANN regression models.

Heart disease

The _Heart_ dataset came from a study of 303 patients referred for coronary
angiography at the Cleveland Clinic. The patients were examined by a set of
tests, including physical examination, electrocardiogram at rest, serum
cholesterol determination and fasting blood sugar determination. The dataset
consists of the examination results combined with historical records. More
details about the data can be found elsewhere [15] [16]. The dataset is
publicly available from the UC Irvine Machine Learning Repository [17]. The
original data include records from several hospitals, in this research only
data from the Cleveland Clinic are used.

Eleven records with missing values were removed from the original data,
resulting in 292 rows. The dataset consists of 14 numeric and categorical
variables, the overview is given in Table 1. The _Class_ variable is used as a
response, and the rest are attributed to predictors.

Table 1: Variables of the _Heart_ dataset, their type and value range. Attribute | Type | Values
---|---|---
Age | numeric, years | 29–77
Sex | categorical | male (200), female (92)
Chest pain type | categorical | abnang (49), angina (23), asympt (138), notang (82)
Resting blood pressure | numeric | 94–200
Cholesterol | numeric | 126–564
Sugar | categorical |  $<120$ (43), $\geq 120$ (249)
Resting ECG | categorical | hyper (145), normal (147)
Max heart rate | numeric | 71–202
Exercise induced angina | categorical | true (95), false (197)
Oldpeak | numeric | 0–6.2
Slope | categorical | down (20), flat (134), up (138)
Number of vessels colored | discrete | 0–3
Thal | categorical | fix (17), normal (161), rev (114)
Class | categorical | healthy (159), sick (133)
|  |

This dataset is used to demonstrate that PV-set augmentation can also be
applied to mixed datasets with categorical variables, and to show how SVD and
PLS versions work for solving binary classification tasks.

###  3.2 ANN regression of Tecator data

As with most of the spectroscopic measurements, the _Tecator_ data are highly
collinear, and SVD decomposition of the matrix with spectra resulted in the
first eigenvalue of 25.6, while starting from the 6th all eigenvalues are
significantly below 0.001. PLS decomposition of the training set resulted in
similar outcomes, therefore, one can assume that using 5-6 latent variables is
enough to capture the majority of the predictors’ variation.

At the same time, preliminary investigation has shown that $c_{k}/c$ values
are within the desired limit of $[0,2]$ for all latent variables up to $A=50$.
Thus, using any number of latent variables between $5$ and $50$ will produce a
reasonable PV-set.

To predict fat content a simple ANN model was employed. The model consisted of
six fully connected (linear) layers of the following sizes: $(100,150)$,
$(150,200)$, $(200,150)$, $(150,100)$, $(100,50)$, and $(50,1)$. The first
value in the parenthesis denotes the number of inputs and the second value
denotes the number of outputs in each layer. The first five layers were
supplemented with the ReLU activation function, while the last layer was used
as the output of the model. All other characteristics of the model are shown
in Table 2.

The model contains approximately 95000 tunable parameters in total, which is
much larger than the number of samples in the training set (170). However, as
already noted, first several hidden layers of ANN usually serve as feature
extractors, so not all parameters will be a part of the regression model
itself.

Table 2: Characteristics of the ANN model for Tecator experiments. Characteristic | Value
---|---
Layer 1 | (100, 150) + ReLU
Layer 2 | (150, 200) + ReLU
Layer 3 | (200, 150) + ReLU
Layer 4 | (150, 100) + ReLU
Layer 5 | (100, 50) + ReLU
Layer 6 | (50, 1)
Optimizer | Adam
Learning rate | 0.0001
Number of epochs | 300
Batch size | 10
|

The columns of the predictors were standardized using mean and standard
deviation values computed for the original training set. The response values
were mean centered only to obtain errors directly comparable in the original
units (%w/w).

First, the ANN training procedure was run 30 times using the original training
set without augmentation. Repeated experiments are necessary because the ANN
training procedure is not deterministic, as explained in the introduction, and
the performance varies from model to model.

At each run the resulting model was evaluated using the independent test set.
The median root mean squared error for the test set (RMSEP) was $4.17$ %w/w
and the coefficient of determination, $R^{2}$, was $0.91$. The best values for
RMSEP and $R^{2}$ were $3.69$ and $0.93$ respectively.

After that, a full factorial experiment was set up to determine how the
augmentation of the training data with PV-sets influences the performance of
the ANN model and which parameters for PV-set generation as well as the number
of generated sets have the largest effect on the performance. The PV-sets were
generated using PLS decomposition. The generation parameters are shown in
Table 3.

All possible combinations of the values from the table were tested (60
combinations in total). For every combination the training/test procedure was
repeated 5 times with full reinitialization, resulting in 300 outcomes.

Table 3: Levels of PV-set generation parameters and number of sets used in the factorial experiments. Parameter | Levels
---|---
Number of PV-sets | 1, 5, 10, 20, 50
Number of LV-s | 5, 10, 20
Number of segments | 2, 4, 10, 20
|
Figure 2: Box and whiskers plots showing variation of RMSEP values of the ANN
model depending on the parameters of PV-sets generation and the number of the
sets (left — number of PV-sets, middle — number of latent variables in PV-
model, right — number of segments in cross-validation resampling). Black
points and numbers show median values.

Figure 2 graphically illustrates the outcomes of the experiment in the form of
boxplots. Every plot shows the variation in the RMSEP for a given value of one
of the three tested parameters. The black point inside each box supplemented
with text shows the corresponding median value.

It is clear that the number of PV-sets used for data augmentation (first plot
in the figure) has the largest effect on the RMSEP. Thus, augmenting the
original training set with one PV-set reduces the median test set error by
approximately 18% (from $4.17$ to $3.41$ %w/w). Adding five PV-sets reduces
the median RMSEP down to $1.58$ %w/w (62% reduction). Using 10 PV-sets further
reduces the RMSEP down to $1.36$ %w/w; however, the effect is smaller.

Changing the number of latent variables does not give a significant difference
in the performance of the models, while the number of segments in cross-
validation resampling shows a small effect. The best model was obtained by
using data augmented with 20 PV-sets generated with 10 latent variables and 10
segments. It could predict fat content in the test set with RMSEP = 0.94 %w/w
($R^{2}=0.995$) which is more than 3 times smaller compared to the best model
obtained without augmentation.

Figure 3: Box and whiskers plots showing variation of RMSEP values of the ANN
model depending on the parameters of PV-sets generation and the number of the
generated sets, obtained using larger learning rate (0.001).

Statistical analysis of the outcomes carried out by N-way analysis of variance
(ANOVA) and Tukey’s honest significance difference (HSD) test confirmed the
significant difference between average RMSE values obtained using different
numbers of PV-sets (ANOVA p-value $\ll 0.01$). However, pairwise comparison
shows no significant difference between 10 vs. 20, 10 vs. 50, and 20 vs 50 PV-
sets (p-value $>0.20$ in for all three pairs).

The number of latent variables has no significant effect (ANOVA p-value
$\approx 0.71$). The number of segments shows only a significant difference
between 2 segments and other choices, favoring the use of 4 segments or above
(p-values for all pairs with 2 segments $\ll 0.01$, for 4 segments vs. 10
segments p-value $\approx 0.26$)

Optimization of the ANN learning parameters (by varying the learning rate and
the batch size) toward reducing RMSEP values for the models trained without
augmentation, made this gap smaller. Figure 3 shows the results obtained using
a learning rate = 0.001. Training the model without augmentation results with
median RMSEP = $2.84$ %w/w. Training the model on the augmented data with 20
PV-sets resulted in a median RMSEP = $1.80$ %w/w. This is still approximately
37% smaller; however, the overall performance of the model trained with this
learning rate is worse.

Figure 4: Box and whiskers plots showing variation of RMSEP values of the ANN
models trained using original training data and augmented data depending on
the learning rate.

This effect is understandable as changing the learning rate and batch size to
improve the performance of the model trained on the original data makes the
model less sensitive to overfitting and local minima problems but also less
flexible, which, in turn, makes the use of augmented data less efficient.

Figure 4 further illustrates this effect. The plot shows the variation in
RMSEP values for ANN models trained using different learning rates. Each
learning rate is represented by a pair of box and whiskers series. The left
(blue) series in each pair illustrates the results obtained using 30 models
trained with the original data. The right (red) series in each pair shows the
results for 30 models trained using the augmented data (20 PV-sets computed
using 10 LVs and 4 segments). New PV-sets were computed at each run to
eliminate possible random effects.

As one can see, at small learning rates ($10^{-5}$—$10^{-3}$) the models
trained with augmented data clearly outperform the models trained using the
original training set. Not only in mean or median RMSEP values but also model-
to-model variation of RMSEP is much smaller, despite the additional
uncertainty introduced by the augmentation (each PV-set is generated using
random splits; hence, the augmented datasets are always different from each
other). Starting from LR = $10^{-2}$ and larger data augmentation does not
show any benefits, but the model trained using the original data also clearly
loses the performance.

This experiment was also repeated using other ANN architectures, including
networks with convolution layers. The gap between the performance varied
depending on the architecture and the learning parameters; however, in all
experiments, a minimum of 20% remained, clearly indicating the benefits of the
PV-set augmentation in this particular case.

###  3.3 ANN classification of Heart data

All categorical variables from the _Heart_ dataset with $L$ levels were
converted to $L-1$ dummy values $[0,1]$, so the matrix with predictors,
$\mathbf{X}$ contained 17 columns in total. The columns of the matrix were
mean centered and standardized. SVD decomposition of the matrix indicates
moderate collinearity with eigenvalues ranging from 3.43 to 0.77 for the first
10 latent variables.

ANN model for classification of the patient conditions used in this research
included six linear layers and one sigmoid layer. The first five layers were
supplemented with the ReLU activation function, while the last layer was used
as the output of the model. The main characteristics of the model are shown in
Table 4. The model has approximately 14700 tunable parameters, while the
original dataset has only 292 objects.

Table 4: Characteristics of the ANN model for Heart dataset classification. Characteristic | Value
---|---
Layer 1 | (17, 34) + ReLU
Layer 2 | (34, 68) + ReLU
Layer 3 | (68, 68) + ReLU
Layer 4 | (68, 68) + ReLU
Layer 5 | (68, 34) + ReLU
Layer 6 | (34, 1)
Layer 7 | Sigmoid
Optimizer | Adam
Learning rate | 0.000001
Number of epochs | 300
Batch size | 10
|

Since this dataset does not contain a dedicated test set, the following
procedure was employed. At each run, the data objects were split into two
groups with healthy persons in one group and sick persons in the other. Then,
75% of records were selected randomly from each group and merged together to
form a training set. The remaining 25% of records formed the test set for the
run. Classification accuracy, which was computed as a ratio of all correctly
classified records to the total number of records in the test set, was used as
the performance indicator.

The experimental design was similar to the _Tecator_ experiments. First, the
training/test procedure was repeated 30 times using the data without
augmentation. Then the ANN models were trained using the augmented data and
tested using the randomly selected test set. The PV-sets for augmentation were
computed using an algorithm based on PLS decomposition applied to the randomly
selected training set. As in the previous chapter, all possible combinations
of PV-set generation parameters from Table 2 were tested with 5 replications
(300 runs in total).

Figure 5: Box and whiskers plots showing variation of test set accuracy of the
ANN models depending on PV-set generation parameters (PV-sets computed using
algorithm based on PLS decomposition, ANN learning rate = $10^{-6}$).

Figure 5 presents the outcomes of the experiment in the form of boxplots,
which show a variation of the test set accuracy computed at each run. One can
clearly notice that the accuracy of the models trained on the data without
augmentation is very low (median accuracy is 0.50), while augmenting the
training data with 20 PV-sets raises the median accuracy to 0.84. The best
model obtained in these experiments is also trained on the augmented data and
has an accuracy of 0.91.

N-way ANOVA and Tukey’s tests for the outcomes have shown that two parameters
— number of PV-sets as well as number of segments in cross-validation
resampling have a significant influence on the accuracy (p-value for both
factors was $\ll 0.001$ in ANOVA test). However, a significant difference was
observed only for the number of segments equal to two, similar to the
_Tecator_ results. Using four or more segments for the generation of PV-sets
shows statistically similar performance to the models trained on the augmented
data.

Figure 6: Box and whiskers plots showing variation of test set accuracy of the
ANN models depending on parameters of PV-set generation (PV-sets computed
using algorithm based on PLS decomposition, ANN learning rate = $10^{-3}$).

The large gap in the performance of the models trained using the original and
augmented data is explained by the initial selection of a low learning rate
($10^{-6}$), which was justified by the results of experiments with _Tecator_
data. Optimization of the learning rate to obtain the best accuracy for the
data without augmentation eliminates this gap as shown in Figure 6. However,
as in the _Tecator_ experiments, the overall performance of the models in this
case is slightly lower, and the median accuracy was approximately 0.79.

Since the _Heart_ data consist of two classes and the classes are balanced
(the number of observations in the classes is comparable), one can apply an
SVD-based algorithm for PV-set generation. In this case, at each run the PV-
sets are generated separately for the sick and healthy groups, and then are
merged together to form the augmented training set. The results obtained for
the learning rate of $10^{-6}$ are shown in Figure 7.

Figure 7: Box and whiskers plots showing variation of test set accuracy of the
ANN models depending on PV-set generation (PV-sets computed using algorithm
based on SVD decomposition, ANN learning rate = $10^{-6}$).

The results are very similar to those obtained using the PLS-based algorithm.
However, in this case, the overall performance of the models trained on the
augmented data was a slightly higher, with a median accuracy of 0.84 for the
models trained with on data augmented with 10, 20 and 50 PV-sets. The best
model had an accuracy of 0.95.

Reducing the learning rate down to 0.001 had the same effect as for the data
augmented using the PLS-based algorithm — the gap between the models trained
with and without augmented data is eliminated however the overall performance
also gets lower with, median a accuracy of approximately 0.79.

##  4 Discussion

The experimental results confirm the benefits of PV-set augmentation, however
optimization of ANN learning parameters is needed to make the benefits
significant. At the same time, optimization of the PV-set generation algorithm
is not necessary for most of the cases. Based on our experiments, we advise
using cross-validation resampling with 5 or 10 splits and a number of latent
variables large enough to capture the majority of variation in $\mathbf{X}$.
In some specific cases one can use tools for quality control of generated PV-
sets described in [10].

It must also be noted that the use of PV-sets for data augmentation is not
always beneficial. Thus, according to our experiments, which are not reported
here, in methods that are robust to overfitting, such as, for example, random
forest (RF), increasing the training set artificially does not have a
significant effect on the model performance. In the case of eXtreme Gradient
Boosting changing training parameters, which regulate the overfitting, such as
a learning rate, maximum depth and minimum sum of instance weight, can have an
effect, but most of the time the effect observed in our experiments was
marginal.

##  5 Conclusions

This paper proposes a new method for data augmentation. The method is
beneficial specifically for datasets with moderate to high degree of
collinearity as it directly utilizes this feature in the generation algorithm.

Two proposed implementations of the method (SVD and PLS based) cover most of
the common data analysis tasks, such as regression, discrimination and one-
class classification (authentication). Both implementations are very fast —
the generation of a PV-set for $\mathbf{X}$ of $200\times 500$ with 20 latent
variables and 10 segments splits requires several seconds (less than a second
on a powerful PC), much less than the training of an ANN model with several
layers.

The method can work with datasets of small size (from tens observations) and
can be used for both numeric and mixed datasets, where one or several
variables are categorical.

## References

  *   * [1] Ratner, A. J., Ehrenberg, H. R., Hussain, Z., Dunnmon, J. & Ré, C.  Learning to compose domain-specific transformations for data augmentation (2017).  1709.01643.
  * [2] Goodfellow, I. J. _et al._ Generative adversarial networks (2014).  1406.2661.
  * [3] Dao, T. _et al._ A kernel theory of modern data augmentation (2019).  1803.06084.
  * [4] Pérez, E. & Ventura, S.  Progressive growing of generative adversarial networks for improving data augmentation and skin cancer diagnosis.  _Artificial Intelligence in Medicine_ 141, 102556 (2023).  URL https://www.sciencedirect.com/science/article/pii/S0933365723000702.
  * [5] Perez, F., Vasconcelos, C., Avila, S. & Valle, E.  Stoyanov, D. _et al._ (eds) _Data augmentation for skin lesion analysis_.  (eds Stoyanov, D. _et al._) _OR 2.0 Context-Aware Operating Theaters, Computer Assisted Robotic Endoscopy, Clinical Image-Based Procedures, and Skin Image Analysis_ , 303–311 (Springer International Publishing, Cham, 2018).
  * [6] Iglesias, G., Talavera, E., González-Prieto, Á., Mozo, A. & Gómez-Canaval, S.  Data augmentation techniques in time series domain: a survey and taxonomy.  _Neural Computing and Applications_ 35, 10123–10145 (2023).  URL https://doi.org/10.1007/s00521-023-08459-3.
  * [7] Sáiz-Abajo, M., Mevik, B.-H., Segtnan, V. & Næs, T.  Ensemble methods and data augmentation by noise addition applied to the analysis of spectroscopic data.  _Analytica Chimica Acta_ 533, 147–159 (2005).  URL https://www.sciencedirect.com/science/article/pii/S000326700401428X.
  * [8] Chadebec, C. & Allassonnière, S.  Data augmentation with variational autoencoders and manifold sampling (2021).  2103.13751.
  * [9] Kucheryavskiy, S., Zhilin, S., Rodionova, O. & Pomerantsev, A.  Procrustes cross-validation—a bridge between cross-validation and independent validation sets.  _Analytical Chemistry_ 92, 11842–11850 (2020).
  * [10] Kucheryavskiy, S., Rodionova, O. & Pomerantsev, A.  Procrustes cross-validation of multivariate regression models.  _Analytica Chimica Acta_ 1255, 341096 (2023).  URL https://www.sciencedirect.com/science/article/pii/S0003267023003173.
  * [11] de Jong, S.  Simpls: An alternative approach to partial least squares regression.  _Chemometrics and Intelligent Laboratory Systems_ 18, 251–263 (1993).
  * [12] Paszke, A. _et al._ in Pytorch: An imperative style, high-performance deep learning library 8024–8035 (Curran Associates, Inc., 2019).  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * [13] Harris, C. R. _et al._ Array programming with NumPy.  _Nature_ 585, 357–362 (2020).  URL https://doi.org/10.1038/s41586-020-2649-2.
  * [14] Borggaard, C. & Thodberg, H. H.  Optimal minimal neural interpretation of spectra.  _Analytical Chemistry_ 64, 545–551 (1992).
  * [15] Detrano, R. _et al._ International application of a new probability algorithm for the diagnosis of coronary artery disease.  _The American Journal of Cardiology_ 64, 304–310 (1989).  URL https://www.sciencedirect.com/science/article/pii/0002914989905249.
  * [16] Detrano, R. _et al._ Bayesian probability analysis: a prospective demonstration of its clinical utility in diagnosing coronary disease.  _Circulation_ 69, 541–547 (1984).
  * [17] Janosi, A., Steinbrunn, W., Pfisterer, M. & Detrano, R.  Heart Disease Dataset (1988).