Application of Multidimensional Statistical Methods to the Hydrochemical Study with R Software

2.1 Environment and study site

The study is only conducted in the Oued Inaouen watershed, which has a total measured area of 3,396 km² and is located upstream of the Idriss I dam (Figure 1) [15]. The Sebou basin upstream of the Idriss I dam makes up 8.3% of the total area. Its boundaries are the Upper Sebou in the south, the Upper Ouergha in the north, and the Middle Mouloya watershed in the east. On the right bank of this basin, it primarily drains the marly formations of the Pre-Rifain relief. On the left bank, it primarily drains the carbonate formations of the Middle Atlas Causse. The Oued Inaouen receives the Oued Flows of Amlil, Larbaa, and Lahdar on its right bank. These tributaries gather the runoff from the pre-rific hills. On its left bank, the Oued Inaouen receives tributaries from the Middle Atlas limestones Matmata, Zerg, Bouzemlane, and Bouhlou, which are often karstic in this area and are partly fed by the Tazekka central massif. Fifty-two wastewater discharges from the city of Taza are spread. There could be up to 5,000,000 m³ of wastewater released into waterways each year. Given the population growth in the city of Taza, this flow is increasing year over year [3].

1.png

Figure 1. Geographical location of the study area

2.2 Samples

The sampling campaigns were carried out mainly on fixed dates throughout the study period from February 2014 to December 2015. The frequency of measurement is seasonal (spring, summer, autumn and winter) to obtain a reasonably representative picture of water quality and its seasonal, annual and multi-year evolution [16].

Our database contains 100 surface water samples (observations) collected in the province of Taza; water samples were collected, transported, and stored according to the National Drinking Water Office's (ONEP) policy and procedures [17]. A portion of the analysis was done there. Another portion was done at the Regional University Centre of Interface (CURI) attached to the Sidi Mohamed Ben Abdellah University (USMBA) of Fez at the time of sampling to avoid the change and temperature, conductivity, pH using a multi-parameter analyzer Type CONSORT - Model C535, and dissolved oxygen by the titration method of Winkler. The methods used at CURI Laboratory: volumetry for bicarbonates, chlorides, calcium and magnesium; molecular absorption spectrophotometry for sulfates, nitrates, nitrites, ammonium ions and orthophosphates; and flame spectrophotometry for sodium and potassium. The four heavy metals are analyzed by atomic absorption spectrum (AAS).

2.3 Physicochemical variables studied

The variables retained for this statistical study are the eight major ions: Sodium (Na⁺), Bicarbonate (HCO₃^-), Magnesium (Mg²⁺), Chlorides (Cl^-), Potassium (K+), Calcium (Ca²⁺), Sulfates (SO₄^2-) and Total Alkalinity (as CaCO₃). The means and standard deviations of these studied parameters are reported in the Table 1.

Table 1. The statistical summary of data

2.4 Data standardization

Variables in the principal component analysis are frequently normalized [18, 19]. This is especially advised when the variables are measured in various units (such as kilograms, kilometres, centimetres, etc.); otherwise, the PCA result will be significantly impacted.

Making the variables similar is the goal. The variables are often normalized to have a mean of zero, a standard deviation of one, and a standard deviation of one [20, 21].

Technically speaking, the method entails converting the data by subtracting a reference value (the variable's mean) from each value and dividing the result by the standard deviation. The data left over after this transformation is referred to as centred-reduced data. Normalized PCA is the name of the PCA used to analyze these modified data [22].

Before PCA and clustering studies, data normalization was widely utilized in gene expression data analysis.

The information can be altered as follows after normalizing the variables [23]: $\frac{x_i-\operatorname{mean}(x)}{s d(x)}$.

Where mean(x) is the average of the values of x, and sd(x) is the standard deviation.

2.5 Statistical processing of data

As a multidimensional descriptive statistical technique, Principal Component Analysis (PCA) enables the analysis and presentation of a data set containing individuals characterized by quantitative variables [24, 25]. A statistical approach that allows the investigation of multivariate data is this one (data with several variables). It is possible to think of each variable as a different dimension [26]. It could be very challenging to show your dataset in a multidimensional hyperspace if it has more than three variables. The crucial data in a multivariate data table are extracted and visualized using principal component analysis [27]. This data is combined via PCA into a small number of new variables called principal components. The original variables are linearly combined to create these new variables [28]. There are fewer or the same number of principal components as there were original variables. The total variance or inertia of a data set is the information it holds [29]. The goal of PCA is to find the principal axes or principal components [30], otherwise known as the directions along which the data vary most. In other words, PCA preserves as much information as possible while condensing the dimensions of multivariate data to two or three principal components that can be graphically represented [31]. Giving a graphical representation of the data and the relationships between the variables is challenging when there are many observed variables [32].

The most important parameters that describe the quality of surface waters and illustrate their variability can be found using this method, which has been widely used to interpret hydrochemical data of hydro systems. Principal component analysis was the foundation for the statistical analysis.

R software (version 3.6.3) was used to generate the intermediate correlation matrix and project the variables into the space of the Dim1 and Dim2 axes.

2.6 PCA procedure

PCA replaces a family of variables with new variables called principal components (PC) [33-35]. The latter is of maximum variance and uncorrelated two by two. They are linear combinations of the original variables. Let us consider a set of data collected during a normal operation of the system under study. A matrix can represent these data [36].

$\mathbf{X}=[\mathbf{x}(1), \cdots, \mathbf{x}(N)]^T \in \mathbb{R}^{N \times m}$        (1)

where, $N$ is the number of observations and $m$ is the number of measured variables. Each row of the data matrix $\mathbf{X}$ represents an observation in the form of a vector of measurements collected at a time $k$ generally centred [37].

$\mathbf{x}(k)=\left[x_1(k), \cdots, x_m(k)\right]^T \in \mathbb{R}^m$          (2)

where, $x_j(k)$ with $j=\{1, \cdots, \mathrm{m}\}$ represents the measurement of the variable $j$ at time $k$. By definition, the covariance matrix is given in the research [37]:

$\Sigma=\mathbb{E}\left\{\mathbf{x} \mathbf{x}^T\right\}=\frac{1}{N} \mathbf{X}^T \mathbf{X} \in \mathbb{R}^{m \times m}$              (3)

According to the principle of PCA, it is assumed that a vector of components $\hat{\mathbf{t}} \in \mathbb{R}^{\ell}$ is associated with each observation vector whose representation it optimizes in the sense of minimizing the estimation error of $x$ or the maximization of the variance of $\hat{\mathbf{t}}$. At each time $k$, the vectors $\hat{\mathbf{t}}$ and $\mathbf{x}$ are linked by a linear transformation of type $\hat{\mathbf{t}}(k)=$ $\hat{P}^T \mathbf{x}(k)$ such that the transformation matrix $\hat{P} \in \mathbb{R}^{m \times \ell}$ verifies the orthogonality condition $\hat{P}^T \hat{P}=\mathbf{I}_{\ell} \in \mathbb{R}^{\ell \times \ell}$.

The columns of the matrix $\hat{P}$ are the vectors of an orthonormal basis of a subspace $\mathbb{R}^{\ell}$ of reduced representation of the original data. The linear transformation results in the projection of the original data expressed in the space of dimension $m$ to an orthogonal subspace of dimension $\ell$. The components $t_j(k)$ with $j=\{1, \cdots, \ell\}$ of the vector $\hat{\mathbf{t}}(k)$ are the projections of the elements of the data vector $\boldsymbol{x}(k)$ in the subspace $\mathbb{R}^{\ell}$.

The representation optimization based on the projection matrix $\hat{P}$ is obtained by minimizing the squared estimation error $o \mathbf{x}$ f. Let us note by $\hat{P}$ the optimal representation matrix, which can be given in study [37]:

$\hat{P}=\arg \min _{\grave{P}}\left\{J_e(\grave{P})\right\}$                (4)

where, $J_e$ is the criterion for the estimation error by PCA, which should be minimized [38]. Under the Orthogonality constraint of the projection matrix $\hat{P}$ We can write [39]:

$\begin{aligned} J_e(\hat{P}) & =\mathbb{E}\left\{\|\mathbf{x}-\hat{\mathbf{x}}\|^2\right\}=\mathbb{E}\left\{\left\|\mathbf{x}-\hat{P} \hat{P}^T \mathbf{x}\right\|^2\right\} \\ & =\mathbb{E}\left\{(\mathbf{x}-\hat{P} \hat{\mathbf{t}})^T(\mathbf{x}-\hat{P} \hat{\mathbf{t}})\right\}=\mathbb{E}\left\{\mathbf{x}^T \mathbf{x}-\hat{\mathbf{t}}^T \hat{\mathbf{t}}\right\} \\ & =\mathbb{E}\left\{\operatorname{tr}\left(\mathbf{x x}^T\right)-\hat{\mathbf{t}}^T \hat{\mathbf{t}}\right\}=\operatorname{tr}\{\Sigma\}-\mathbb{E}\left\{\hat{\mathbf{t}}^T \hat{\mathbf{t}}\right\} \\ & =\operatorname{tr}\{\Sigma\}-J_v(\hat{P})\end{aligned}$      (5)

where, $\operatorname{tr}\{$. denotes the trace of a square matrix. Since the term $\operatorname{tr}\{\Sigma\}$ is a constant, minimizing the criterion $J_e$ is the same as maximizing the $J_v$ Given by [39]:

$\begin{aligned} J_v(\hat{P})=\mathbb{E}\left\{\hat{\mathbf{t}}^T \hat{\mathbf{t}}\right\} & =\mathbb{E}\left\{\sum_{j=1}^{\ell} t_j^2\right\}=\sum_{j=1}^{\ell} \mathbb{E}\left\{t_j^2\right\} \\ = & \sum_{j=1}^{\ell} \operatorname{Var}\left\{t_j\right\}\end{aligned}$           (6)

From the previous equation, maximizing the $J_v$ equivalent to maximizing the variance of the $t_j$. Thus, the optimization problem is reformulated as follows [37]:

$\hat{\boldsymbol{P}}=\arg \min _{\hat{P}}\left\{J_e(\hat{P})\right\}=\arg \max _{\hat{P}}\left\{J_v(\hat{P})\right\}$            (7)

To determine the column vectors of the matrix $\hat{P}$, we note by $t \in \mathbb{R}$ the projection of the data vector $\mathbf{x}$ along a direction represented by a unit vector $p \in \mathbb{R}^m$. The scalar product $t=$ $\mathbf{x}^T p=p^T \mathbf{x}, \mathbf{x}$ obtains the component $\mathrm{t}$ under the constraint $\|p\|^2=p^T p=1$. In particular, it represents a new variable with a mean and a variance that depend on the statistical properties of $\mathbf{x}$ as follows [39-41]:

$\mathbb{E}\{t\}=\mathbb{E}\left\{p^T \mathbf{x}\right\}=p^T \mathbb{E}\{\mathbf{x}\}=0$           (8)

$\begin{aligned} \operatorname{Var}\{t\} & =\mathbb{E}\left\{(t-\mathbb{E}\{t\})^2\right\}=\mathbb{E}\left\{t^2\right\} \\ & =\mathbb{E}\left\{\left(p^T \mathbf{x}\right)\left(\mathbf{x}^T p\right)\right\}=p^T \mathbb{E}\left\{\mathbf{x x}^T\right\} p \\ & =p^T \Sigma p\end{aligned}$                          (9)

The maximization of the projection variance, under the condition of a unit norm of the vector p, represents an equality-constrained optimization problem that the Lagrange function can formalize [40-43]:

$\mathcal{L}(p, \lambda)=J_v(p)-\lambda\left(p^T p-1\right)$ $=p^T \Sigma p-\lambda\left(p^T p-1\right)$          (10)

where, $\lambda \in \mathbb{R}$ is the Lagrange multiplier. Taking into account the symmetry of the matrix $\Sigma$, the vector $p$ maximizes the optimization criterion $J_v$ is a solution of the following system of equations [37]:

$\left\{\begin{array}{l}\partial \mathcal{L}(p, \lambda) / \partial p=\Sigma p-\lambda p=0 \\ \partial \mathcal{L}(p, \lambda) / \partial \lambda=p^T p-1=0\end{array}\right.$    (11)

Consequently, the solution of this system of equations is identified as a problem of estimating normalized eigenvalues and eigenvectors of the matrix $\Sigma$. Such a system of equations admits real solutions of the variables $\lambda$ obtained by solving the following characteristic equation [41]:

$\operatorname{Det}\left\{\Sigma-\lambda \mathbf{I}_m\right\}=0$    (12)

where, Det $\{$.$\} is the determinant of a square matrix. \mathbf{I}_m$ is the identity matrix of order $m$. The solutions of the previous equation represent the eigenvalues of $\Sigma$. To each eigenvalue, $\lambda$ is associated with an eigenvector $p$ verifying $\left(\Sigma-\lambda \mathbf{I}_m\right) p=0$.

This allows us to have $m$ eigenvectors $\mathbf{p}_i$ associated with the $m$ eigenvalues $\lambda_i$ of the matrix $\Sigma$, thus verifying the relation $\Sigma \mathbf{p}_i=\lambda_i \mathbf{p}_i$ with $i=\{1, \cdots, m\}$. In matrix form, such a relationship leads to writing the following [37, 44, 45]:

$\Sigma \mathbf{P}=\mathbf{P B}$     (13)

$\mathbf{P}=\left[\mathbf{p}_1, \cdots, \mathbf{p}_m\right] \in \mathbb{R}^{m \times m}$ represents the data projection matrix. It is orthonormal since its columns correspond to the eigenvectors of $\Sigma$:

$\mathbf{P}^T \mathbf{P}=\mathbf{P P}^T=\mathbf{I}_m \in \mathbb{R}^{m \times m}$     (14)

$\mathbf{B}=\operatorname{diag}\left\{\lambda_1, \cdots, \lambda_m\right\} \in \mathbb{R}^{m \times m}$ represents the diagonal matrix consisting of the diagonal elements of the eigenvalues of $\Sigma$.

From Eqns. (13) and (14), we can deduce that $\mathbf{P}^T \Sigma \mathbf{P}=\mathbf{B}$. This allows us to conclude that the first direction, having a maximum projection variance of $\mathbf{x}$, is carried by the eigenvector $\mathbf{p}_1$ associated with the largest eigenvalue $\lambda_1$. The latter represents the variance of such a direction. The second factorial axis also renders the maximum variance while orthogonal to the first. Its variance $\lambda_2$ is less important than that corresponding to the first direction. Therefore, the diagonal elements of $\mathbf{B}$ are arranged in descending order [46-49]: $\lambda_1 \geq \cdots \geq \lambda_m$

Considering the matrix $\mathbf{P}$, the data vector $\boldsymbol{x}(k)$ can be transformed without any loss of information into a vector of principal components (PC) [36, 37, 39]:

$\mathbf{t}(k)=\left[t_1(k), \cdots, t_m(k)\right]^T=\mathbf{P}^T \mathbf{x}(k) \in \mathbb{R}^m$        (15)

where, the $\mathrm{CP} t_j$ with $j=\{1, \cdots, m\}$ is defined by:

$t_j(k)=\mathbf{p}_j^T \mathbf{x}(k)=\mathbf{x}^T(k) \mathbf{p}_j$            (16)

These are statistically uncorrelated [37, 39]:

$\mathbb{E}\left\{t_i t_j\right\}=\mathbb{E}\left\{\mathbf{p}_i^T \mathbf{x x}^T \mathbf{p}_j\right\}=\mathbf{p}_i^T \Sigma \mathbf{p}_j=0 i \neq j$            (17)

The notation in matrix form allows us to define the matrix of CP as follows [40, 41]:

$\mathbf{T}=[\mathbf{t}(1), \cdots, \mathbf{t}(N)]^T=\mathbf{X P} \in \mathbb{R}^{N \times m}$        (18)

The determination of the data vector $x(k)$ à from the associated vector of CP $\mathbf{t}(k)$ is given by:

$\mathbf{x}(k)=\operatorname{Pt}(k)=\sum_{j=1}^m \mathbf{p}_j t_j(k)$          (19)

The data reduction is performed through the $\ell$ first $\mathrm{CP}$ with the largest variances. As a result, the $\ell$ first eigenvectors form the reduced vector subspace for the initial data. The estimation of $\hat{\mathbf{x}}(k)$ of the data vector $\boldsymbol{x}(k)$ in this reduced subspace (often called representation or principal subspace and denoted $\hat{S})$ is given by references [40, 42, 43, 50]:

$\hat{\mathbf{x}}(k)=\hat{\mathbf{P}} \hat{\mathbf{t}}(k)=\hat{\mathbf{P}} \hat{\mathbf{P}}^T \mathbf{x}(k)=\hat{\mathbf{C}} \mathbf{x}(k)$           (20)

where, the optimal representation matrix expressed in Eq. (7) is defined as follows:

$\hat{\mathbf{P}}=\left[\mathbf{p}_1, \cdots, \mathbf{p}_{\ell}\right] \in \mathbb{R}^{m \times \ell}$     (21)

$\hat{\mathbf{t}}(k)=\hat{\mathbf{P}}^T \mathbf{x}(k) \in \mathbb{R}^{\ell}$ represents the vector of the $\ell$ first $C P$. The matrix $\hat{\mathbf{C}} \in \mathbb{R}^{m \times m}$ thus characterizes the PCA model.

However, dimension reduction usually results in a loss of information that is recovered in a residual vector $\tilde{\mathbf{x}}(k)$. The latter is expressed in a residual subspace $\tilde{\mathcal{S}}$ consisting of the remainder of the CPs associated with the $(m-\ell)$ last eigenvectors [37, 39, 41]:

$\tilde{\mathbf{x}}(k)=\tilde{\mathbf{P}} \tilde{\mathbf{t}}(k)=\tilde{\mathbf{P}} \tilde{\mathbf{P}}^T \mathbf{x}(k)=\tilde{\mathbf{C}} \mathbf{x}(k)$      (22)

with

$\tilde{\mathbf{P}}=\left[\mathbf{p}_{\ell+1}, \cdots, \mathbf{p}_m\right] \in \mathbb{R}^{m \times(m-\ell)}$               (23)

and

$\tilde{\mathbf{C}}=\tilde{\mathbf{P}} \tilde{\mathbf{P}}^T=\mathbf{I}_m-\hat{\mathbf{C}}$       (24)

The matrix $\tilde{\mathbf{C}} \in \mathbb{R}^{m \times m}$ Describes the residual model. We can see here that PCA is a modelling approach that allows us to obtain a PCA model of a studied system.

The interpretation of the principle of PCA modelling represents a partitioning of the space $\mathbb{R}^m$ of measurements $\mathbf{x}(k)$ into the main subspace $\hat{\mathcal{S}}$ and a residual subspace $\tilde{\mathcal{S}}$. Therefore, the vector of measures $\mathbf{x}(k)$ is decomposed as follows:

$\mathbf{x}(k)=\hat{\mathbf{x}}(k)+\tilde{\mathbf{x}}(k)$          (25)

In particular, a geometric property of orthogonality between the estimated and the residual vector is always verified since [39, 41, 42]:

$\tilde{\mathbf{C}} \hat{\mathbf{C}}=\hat{\mathbf{C}} \tilde{\mathbf{C}}=\mathbf{0}_m \in \mathbb{R}^{m \times m}$     (26)

This implies that the principal subspace and the residual subspace are orthogonal for all values of $\ell$, thus [37, 40, 43],

$\tilde{\mathbf{x}}^T(k) \hat{\mathbf{x}}(k)=0$     (27)