AMV_chapter.qmd

---
title: "Multivariate Analysis"
format: html
editor: visual
author: ANSALDI Lucile, BARREAU Julie, BONGIBAULT Baptiste, GUINEBRETIERE Aldric, MOREAU Cyrielle
---

# Introduction to Multivariate Analysis

Biologists' data, whether on individuals, populations or communities, are usually presented in the form of rectangular tables, with observations (n) in the rows and variables (p) in the columns.

The graphical representation of these n observations and p variables is easily achieved when there are only 2 or 3 variables (dimension). However, when the number of variables increases, the graphical representation becomes complicated, and Multivariate Analyses come into their own!

The aim of these analyses is to reduce the number of dimensions and examine the structure of the data by answering the following questions:

-   Which observations are similar?

-   Are there observations that stand out? Subgroups?

-   Which variables are correlated?

-   Are there particular links between certain observations and/or variables?

Multivariate analysis methods are therefore used to describe the data and generate hypotheses that can then be tested.

# PART 1: Principal Component Analysis (PCA)

## Introduction

PCA can be used to process a measurement table with :

-   In row: the **n observations**.

-   In column: the **p quantitative variables**.

The 2 sets (observations and variables) are **totally distinct and non-interchangeable**. In other words, if you interchange the rows and columns (with the variables in the rows and the observations in the columns), the table no longer has the same meaning.

In PCA tables, the mean of a column has a meaning, while the mean of a row does not.

The aim of PCA is to achieve the best geometric representation of individuals and variables. To achieve this, we seek to reduce the dimension by finding the best projection plane (subspace) for "best" visualization of the point cloud in reduced space.

However, this reduction must :

-   Two individuals who resemble each other must be close in the representation space.

-   Preserve correlations between variables:

    -   Two variables that are correlated must be represented by vectors forming an acute angle.

    -   Two independent variables are represented by orthogonal vectors.

## Mathematics

Principal Component Analysis (PCA) is a powerful method for simplifying the complexity of highly correlated data. This approach aims to reduce dimensionality in a linear manner, thereby providing a clearer perspective on the underlying structure of the data.

In essence, PCA performs a linear transformation of the original variables to new variables called principal components. These components are carefully selected to capture the majority of the variance present in the dataset, thereby maximizing the separation between different observations.

The space of principal components can be expressed as :

$$
Z_p=\Sigma^p_{j=1}ø_j\times X_j
$$Where :

-   $Z_p$ is the principal component $p$

-   $Ø_j$ is the weight vector comprising the $j$ weights for principal component $p$, i.e., the coefficients of the linear combination of the original variables from which the principal components are constructed.

-   $X_j$ is the standardized predictor, meaning it has a mean equal to 0 and a standard deviation of 1.

The process begins with a matrix $Y$, where each row represents an observation and each column represents a continuous variable.

![](table-inds-vars.png){fig-align="center" width="230"}

First, we normalize the observations as follows :

$$
Y_{std} = \frac{y_i-\overline{y}}{\sigma_y}
$$ *; which is equivalent to centering as in:*

$$
y_c= [y_i-\overline{y}]
$$ *; then scaling as in:*

$$
 y_s = \frac{y_i}{\sigma_y}
$$

We can then calculate the variance-covariance matrix:

$$
R = cov(Y_{std})
$$The decomposition of this matrix provides the matrix $U$ of eigenvectors, which contains the principal components.

The eigenvectors, defined as the principal components, have associated eigenvalues. The graphical representation of the distances of observations to these eigenvectors illustrates how each principal component captures the variation in the data.

![](Capture%20d’écran%202023-11-15%20à%2016.45.12.png){fig-align="center" width="190"}

![](Capture%20d’écran%202023-11-15%20à%2016.45.20.png){fig-align="center" width="191"}

The orthogonal aspect of the principal components becomes evident when observing the second principal component, which, in turn, maximizes the variance in the data. These principal components, each aligned along an orthogonal axis, form a new coordinate system in which the information is more clearly structured.

In addition to representing the variance of the principal components, the eigenvectors also provide a relative understanding of their influence. This influence is calculated by dividing the eigenvalues by the total sum of these values, thus providing context on the contribution of each principal component $v_k$ to the overall information.

$$
Explained\ variance\ of\ v_k = \frac{\lambda_{v_k}}{\Sigma_{i=1}^p\lambda_v}
$$The conclusion of the PCA involves obtaining the coordinate matrix $F$. This matrix results from the multiplication of matrix $U$ by the standardized matrix $Y_{std}$, enabling a rotation of the new data space. This process aligns the data along the principal components, providing a simplified yet informative view of the underlying structure of the initial data.

## Interpretation

Here's an example of how PCA can be applied. *dataset presentation* *package needed*

### a) Importing the dataset :

```{r}
# Loading the "iris" dataset available on R :
library(factoextra)
data("decathlon2") 
data=decathlon2[decathlon2$Competition=="OlympicG",c(1:10,12)]
summary(data)
```

In this data set, we have 14 individuals (athletes) on whom 10 quantitative variables (performances in different sports disciplines) have been recorded:

-   X100m

-   Long.jump

-   Shot.put

-   X400m

-   X110m.hurdle

-   Discus

-   Pole.vault

-   Javelin

-   X1500m

If NAs are present in the data. In order to carry out the PCA without any problems with NAs, they can be omitted from the analysis with the following command:

```{r}
# data=na.omit(object = data)
```

We want to find out whether particular links between certain observations and/or variables can be observed in this data set.

### b) Study of correlations :

Obtain the correlation matrix between variables using the **cor()** function:

```{r}
cor(data)
```

Variables are correlated if their value is greater than 0.9 (as a general rule). Here, we can see that no variables are correlated.

Here's another way of visualizing correlations between variables: (useful when you have a lot of variables, as here):

```{r}
abs(cor(data))>0.9
```

If TRUE (outside the diagonal), then both variables are correlated. If two variables are correlated, one must be removed. The choice of deleting one of the two correct variables is arbitrary and depends on the question being asked.

To remove a variable from the data set, use the following function: **data=data\[,-(column_of_the_variable_to_remove)}\]**

### c) Performing PCA :

To run the PCA, you need to load the following two packages: **factoextra** and **FactoMineR**

Then, to perform the PCA on R, you can use the function : **PCA()**. Remember to store the result of this PCA in a new variable, so that you can easily retrieve the PCA information you need for subsequent interpretation.

```{r, echo=FALSE, hide=TRUE}
library(FactoMineR)
library(factoextra)
PCA=PCA(data)
```

### d) Interpretation of outputs :

#### 1) Inertia and choice of axes :

```{r}
# Recover the eigenvalues of each axis and the inertia carried by the principal components
PCA$eig

# Graphical display of the inertia of each axis : 
fviz_eig(PCA, addlabels = TRUE)
```

To determine the number of axes used in the PCA analysis, we need to identify the jump in variance explained by the different axes. Here we can see that the first axis represents 43.2% of the variance, the second axis 22%, the third axis 14.1%, ...

We can therefore see that the jump in variance explained by the different axes is between the first and second axes. The difference between the variance explained by the axes other than the first is negligible compared to the difference between the first axis and the others.

For the PCA interpretation, we therefore retain the first two axes, which together explain 65.19% of the variance.

#### 2) Interpretation of the biological meaning of the axes :

In PCA, the axes are interpreted according to the columns (i.e. variables) via the correlation circle and the absolute contributions of the variables.

To obtain this information in R, use the following commands:

```{r}
# Obtain the absolute contributions of variables :
PCA$var$contrib

# Obtain the circle of correlations :
fviz_pca_var(PCA)

# Graphical representation of absolute variable contribution values for an axis: 
# Axis 1 case : 
fviz_contrib(PCA, choice = "var", axes = 1)

# Axis 2 case: 
fviz_contrib(PCA, choice = "var", axes = 2)
```

The absolute contribution shows how the initial variable contributes to the formation of the axis. To find out from these results which variables contribute to the formation of the synthetic axis, we use the threshold of $\frac{1}{nb\_variables}\times100$ (in general). For a variable to contribute, its contribution value must be greater than this threshold. This threshold value corresponds to the red dotted line on the graphs showing the absolute contributions of the variables along the axes.

From these results, we can see that axis 1 is mainly represented by the following variables:

-   Points, Discus, Shotput, Long.jump, X100m, High.jump

Axis 2 is represented by the variables:

-   X1500m, Pole.vault, X110m.hurdle, X100m

The correlation circle also shows that some variables are related and others independent. Two related variables are represented by vectors forming an acute angle. Two independent variables are represented by orthogonal vectors. Two negatively related variables are represented by two opposite vectors.

#### 3) Projection of individuals to create a typology based on the axes:

For this, we use the relative contributions of the individuals. These relative contributions can be used to create a typology of individuals based on the axes. The relative contribution of an individual corresponds to the share of information on the axis explained by the point and supported on the axis (in percentage).

To obtain the relative contributions of individuals on R, use the following commands:

```{r}
# Relative contributions of individuals :
PCA$ind$cos2

# Obtaining the projection of individuals : 
fviz_pca_ind(PCA)

```

The individuals who most explain Axis 1 are :

-   Sebrle, Clay, Karpov, Schwarzl, Drews, Schoenbeck

Those who most explain axis 2 are :

-   Macey, Warners, Zsivoczky, Barras

### e) Biological conclusion :

Thus, thanks to the results of the absolute contributions of the variables and the relative contributions of the individuals, we can represent the projection of the individuals and variables as follows:

```{r}
fviz_pca_biplot(PCA)
```

Knowing the biological interpretation of the axes :

![**PCA interpretation**](PCA_interpretation.png)

We can therefore see that the individuals with the highest scores at this tournament are those who had the best results in the discus, shot put, long jump and high jump events. In addition, having good results in the 100m race doesn't give you a good ranking.

Thanks to PCA, we can see that the final score is more linked to the score of certain disciplines than others. And therefore, that links between certain variables and individuals can be observed.

# PART 2: Correspondence Analysis (CA)

## Introduction

CA deals with **tables that cross two categorical variables with several modalities**, such as contingency tables, frequency tables, ...

In these tables, rows and columns are **interchangeable** and the **variables are all in the same units**.

The aim of CA is to **reveal structures and associations between these qualitative variables** by representing them in a multidimensional space. It **allows us to visualize relationships between categories of different variables, and to highlight trends, groupings or disparities**. It is particularly useful for analyzing categorical data, when you want to understand the underlying structure of the data and identify associations between categories.

Unlike PCA, for a CA, we're not interested in the distances between points, but in the distances between profiles of different modalities. A CA is therefore simply a PCA on the profiles, by equipping the space with a suitable distance: **the** $\chi^2$ **distance**. With this distance, the weight of rows (or columns) is relativized, but not cancelled out, and the symmetry between rows and columns is preserved.

## Mathematics

Correspondence Analysis (CA) represents a statistical method aimed at exploring the relationship between two sets of categories. In contrast to PCA, which applies to continuous variables, CA is designed for categorical variables, allowing the analysis of the structure of a contingency table.

![](Table_de_contingence.png){fig-align="center" width="261"}

In essence, Correspondence Analysis (CA) seeks to graphically represent associations between rows and columns in a contingency table. It helps reveal trends, groupings, or oppositions within categories, providing insight into the structure of dependence between variables.

The process begins with a contingency table showing the observed frequencies of co-occurrence between different categories. These data are then transformed to obtain a probability table where the frequency of each co-occurrence is calculated as follows:

$$
f_{ij} = \frac{X_{ij}}{n}
$$With $f_{ij}$ representing the probability that events $i$ and $j$ occur jointly, $X_{ij}$ the number of times this co-occurrence has been observed, and $n$ the total number of observations.The marginal probabilities of the columns $(f_{.j})$ and rows $(f_{i.})$ are calculated as follows:

$$
f_{.j} = \sum_{i=1}^{I}f_{ij}\\\\\\  
f_{i.}=\sum_{j=1}^{J}f_{ij}
$$

### - Profiles

-   Row profiles

    The row profile is $L_i = (\frac{f_{i1}}{f_{i.}},...,\frac{f_{ij}}{f_{i.}})$. In matrix terms, we will set up $D_n=diag(f_{i.})$ such that $L = D_n^{-1}F$. We can then define the average row profile, obtained by summing the columns of our frequency table :

    $$
    (\sum_{i=1}^nf_{i.}\frac{f_{i1}}{f_{i.}},...,\sum_{i=1}^nf_{i.}\frac{f_{ij}}{f_{i.}})=(f_{.1},...,f_{.j})
    $$

-   **Column profiles**

    The column profile is $C_j=(\frac{f_{1j}}{f_{.j}},...,\frac{f_{ij}}{f_{.j}})$. As with row profiles, it can be expressed and calculated in matrix terms by setting: $D_p=diag(f_{.j})$, such that $C=D_p^{-1}F^T$. We can then define the average column profile by summing the rows of our frequency table:

    $$
    (\sum_{j=1}^pf_{.j}\frac{f_{1j}}{f_{.j}},...,\sum_{j=1}^pf_{.j}\frac{f_{ij}}{f_{i.}})=(f_{1.},...f_{i.})
    $$

### - Distances between profiles

One can measure the distance between two row profiles by$$d^2(i,i')=\sum_{j=1}^p(\frac{f_{i j}}{f_{i.}}-\frac{f_{i'j}}{f_{i'.}})^2$$

But this distance does not take into account the importance of each column. A more sensible choice is to use the $\chi^2$ distance, which does take the importance of each column into consideration.

$d^2(i,i')=\sum_{j=1}^p\frac{1}{f_{.j}}(\frac{f_{i j}}{f_{i.}}-\frac{f_{i'j}}{f_{i'}})^2$

Let's use the notation $x$ to refer to the vector of distances between two rows $L_{11}-L_{21}$ or between two columns $C_{11}-C_{12}$. Adopting matrix notation, the square distance of the $\chi^2$ is of the form:

$$
x^TD_p^{-1}x = \sum_{j=1}^p\frac{x^2_j}{f_{.j}}
$$for a row point $x\ \epsilon\ \mathbb{R}^n$ .

### - Independence

In the theory of tests, the objective is to determine whether there is a link between the two qualitative variables. The null hypothesis $H_0$ is the independence between the two variables (no effect). The alternative hypothesis $H_1$ is the presence of a link between the two variables.

The association between the two categorical variables is assessed through the discrepancy between the observed data and the independence model.This independence model is calculated based on the following equation:

$$
P(A\ \cap\ B) = P(A) \times P(B)
$$Therefore, the joint probability is the product of the marginal probabilities:

$$
f_{ij} = f_{i.} \times f_{.j}
$$One can then calculate the discrepancy between the observed data $(f_{ij})$ and the independence model:

-   The significance of the association is measured with a $\chi^2$ test

    $$
    \chi_{obs}^2 = \sum_{i=1}^I\sum_{j=1}^J\frac{(obs.headcount - th.headcount)^2}{th.headcount}=\sum_{i=1}^I\sum_{j=1}^J\frac{(n\times f_{ij} - n\times f_{i.}\times f_{.j})^2}{n\times f_{i.}\times f_{.j}}
    $$

-   The strength of the association $Ø^2$ is measured by the discrepancy between theoretical and observed probabilities.

$$
Ø^2 = \sum_{i=1}^I\sum_{j=1}^J\frac{(obs.prob - th.prob)^2}{th.prob}
$$

### - Binary Correspondence Analysis

By conducting a correspondence analysis, we aim to simultaneously represent the row profiles belonging to $\mathbb{R}^p$ and the column profiles belonging to $\mathbb{R}^n$ of a frequency table. As this representation is done in the Cartesian plane, it will be obtained through a double principal component analysis.

The two principal component analyses are not directly performed on the variables but rather on the row profiles (direct analysis) and the column profiles (dual analysis) presented earlier. Moreover, binary correspondence analysis incorporates the notion of column (or row) weights and the $\chi^2$ distance.

#### - Direct Analysis (Rows)

The direct analysis is performed on the row profiles. The rows of $L = D_n^{-1}F$ are elements of $\mathbb{R}^p$. We aim to represent them in this space using the distance function $x^TD_p^{-1}x$.

In the direct analysis, we seek the vector $u\ \epsilon\ \mathbb{R}^p$ such as :

$$
(u^TD_p^{-1}F^TD_n^{-1})D_n(D_n^{-1}FD_p^{-1}u)
$$

*With* $F$ *: the observed frequency matrix*

is maximal, given that $u^TD_p^{-1}u=1$.

The solution is given by the principal eigenvector of

$$
D_pD_p^{-1}F^TD_n^{-1}FD_p^{-1}=F^TD_n^{-1}FD_p^{-1}\equiv S
$$

#### - Dual Analysis (Columns)

Dual analysis is performed on the column profiles. The rows of $C = D_p^{-1}F^T$ are elements of $\mathbb{R}^n$. We aim to represent them in this space using the distance function $x^TD-n^{-1}x$.

In dual analysis, we seek the vector $v\ \epsilon\ \mathbb{R}^n$ such as

$$
(v^TD_n^{-1}FD_p^{-1})D_p(D_p^{-1}F^TD-n^{-1}v)
$$is maximal, given that $v^TD_n^{-1}v=1$.

The solution is given by the principal eigenvector of

$$
D_n(D_n^{-1}FD_p^{-1}F^TD_n^{-1})=FD_p^{-1}F^TD_n^{-1}\equiv T
$$

### - Projection of both analyses onto the same plane

The usual goal of correspondence analysis is to produce a 2-dimensional graph that summarizes the information contained in the frequency table and highlights different interesting associations. In other words, we want to somehow overlay the figures resulting from dual and direct analyses.

To present both analyses in the same plane, it is essential to ensure that both PCAs project the data into the same dimensions. The transition relationships and the concept of the center of gravity presented in the following sections allow us to ensure that we are indeed projecting the "individuals" from dual and direct analyses into the same dimensions.

#### - Transition Relationships

The following relationships

$$
D_n^{-1}F\varphi_j=\sqrt{\lambda_j}\Psi_j\equiv \hat{\Psi}_j
$$

$$
D_p^{-1}F^T\Psi_j=\sqrt{\lambda_j}\varphi_j\equiv \hat{\varphi}_j
$$

*With* $\varphi_j$ : *the* eigenvector *associated with the eigenvalue* $\lambda_j$ *of the matrix F;* $\Psi_j$ : *the* modality *coordinates of the row variable in the row space defined by the eigenvalues;* $\hat{\Psi}_j$ *the* contrinution *of the component* $j$ *to the row variable in the row space and* $\hat{\varphi}_j$ *the* contrinution *of the component* $j$ *to the column variable in the column space*

establish an almost barycentric link between the two types of analysis, meaning that the coordinates of points in one space are proportional to the components of the factor from the other space corresponding to the same eigenvalue.

Verification: The first relationship is confirmed as follows.

$$
D_n^{-1}F\varphi_j=D_n^{-1}(FD_p^{-1}u_j) = D_n^{-1}(\sqrt{\lambda_j}v_j)=\sqrt{\lambda_j}\Psi_j
$$The second identity is demonstrated in a similar manner.

The consequences of the preceding transition relationships are that...

1.  $1=\lambda_1 \geq \lambda_2 \geq … \geq \lambda_p \geq 0$

2.  $(f_1,…,f_p)^T$ is an eigenvector associated with the eigenvalue $\lambda_1=1$ of $S=F^TD_n^{-1}FD_p^{-1}$.

    Note that since the first eigenvalue is always equal to one, many software packages only mention the other $p-1$ eigenvalues.

#### - Center of Gravity

To project the observations onto the same plane, we need to determine the origin of the graph.

The center of gravity of the rows is defined by $G_L=(g_1,…,g_p)^T$, where

$$
g_j=\sum_{i=1}^nf_i\frac{f_{ij}}{f_i}=\sum_{i=1}^nf_{ij}=f_j, 1\leq j \leq p
$$Similarly, the center of gravity of the columns is defined by

$$
G_C=(f_1,...,f_n)^T
$$Centering of the rows is obtained by calculating

$$
\frac{f_{ij}}{f_i}-g_j=\frac{f_{ij}}{f_i}-f_j=\frac{f_{ij}-f_if_j}{f_i}
$$so that $\sum_{j=1}^p\frac{f_{ij}-f_if_j}{f_i}=0$

for all $i\ \epsilon$ ${1,…,n}$.

The consequence of centering the rows is that the analysis is no longer conducted on $S=F^TD_n^{-1}FD_p^{-1}$ but rather on $S*=(s^*_{jj'})$, where

$$
s^*_{jj'}=\sum_{i=1}^n\frac{(f_{ij}-f_if_j)(f_{ij'}-f_if_j)}{f_if_{j'}}
$$To detect associations between rows and columns, one must establish a connection with the statistic of the $\chi^2$ test.

By definition, $trace(S*)=\sum_{j=1}^p\sum_{i=1}^n\frac{(f_{ij}-f_if_j)^2}{f_if_j}$

## Interpretation

Here's an example of how CA can be applied.

You can import R package using the code

```{r}
# Creating the dataset :
## Step 1: Create eye and hair color vectors
colors_of_eyes <- c("Blue", "Brown", "Green", "Gray", "Black", "Yellow")
hair_colors <- c("Black", "Brown", "Blond", "Redhead", "Chestnut", "Blues", "White", "Pink")

## Step 2: Create a data array with fixed values
fixed_values <- matrix(c(
  50, 35, 32, 24, 0, 0,
  15, 84, 34, 18, 16, 0,
  7, 4, 1, 20, 10, 0,
  0, 18, 0, 21, 34, 6,
  5, 3, 1, 0, 0, 3,
  8, 4, 2, 5, 13, 24,
  9, 10, 0, 12, 0, 10,
  1, 20, 5, 46, 2, 5
), nrow = length(hair_colors), ncol = length(colors_of_eyes))

## Store the matrix in an array variable : 
data <- as.data.frame(fixed_values)

## Name rows and columns
colnames(data) <- colors_of_eyes
rownames(data) <- hair_colors


data
```

This table shows the number of people with each eye color (in columns) and hair color (in rows), for a total of 617 people sampled.

By performing a CA on this dataset, we seek to show the correspondences/oppositions between the different modalities of a variable.

### a) Performing CA :

To run the CA, such as the PCA, you need to load the following two packages: **factoextra** and **FactoMineR**.

Then, to perform the CA on R, you can use the function : **CA()**. Remember to store the result of this PCA in a new variable, so that you can easily retrieve the PCA information you need for subsequent interpretation.

```{r, echo=FALSE, hide=TRUE}
library(FactoMineR)
library(factoextra)
CA=CA(data)
```

In a CA, rows are called **rows** and columns called **col**.

### b) Interpretation of outputs :

#### 1) Inertia and choice of axes :

```{r}
# Recover the eigenvalues of each axis and the inertia carried by the principal components
CA$eig

# Graphical display of the inertia of each axis : 
fviz_eig(CA, addlabels = TRUE)
```

To determine the number of axes used in the CA analysis, we need to identify the jump in variance explained by the different axes. It's the same as for PCA. Here we can see that the first axis represents 55.4% of the variance, the second axis 28.7%, the third axis 9.8%, ...

We can therefore see that the jump in variance explained by the different axes is between the first and second axes. The difference between the variance explained by the axes other than the first is negligible compared to the difference between the first axis and the others.

For the CA interpretation, we therefore retain the first two axes, which together explain 65.19% of the variance.

#### 2) Interpretation of the biological meaning of the axes :

In CA, in contrast to PCA, axes are interpreted either column-wise or row-wise. The choice depends on the initial biological question.

Here we'll interpret the axes according to the lines, i.e. according to the hair colors. To find out which variables contribute to the synthetic axis, we use the threshold of : $\frac{1}{nb\_variables}$

The threshold in our case is: $\frac{1}{8}$.

To obtain this information in R, use the following commands:

```{r}
# Obtaining absolute line contribution values : 
Row_contrib=CA$row$contrib
Row_contrib>(1/8)*100


```

According to these results, the lines explaining the first axis are :

-   Black, Blues, Pink

The lines explaining the second axis are :

-   Chestnut, White

#### 3) Interpretation of relative contributions :

Relative contributions represent the quality of representation of the column/row by the axis.

Either the rows or the columns are interpreted, depending on the choice made in the next step. In our case, we chose the columns, since we interpreted the axes with the absolute contributions of the rows.

```{r}
# Obtain relative column contribution values :
CA$col$cos2

```

Blue and yellow eye colors are well represented by axis 1. Grey eye color is represented by axis 2.

```{r}
# Relative line contribution values : 
CA$row$cos2

```

This verifies that the modalities that contribute most to the formation of the first two axes all have a good representation quality (not always the case).

### c) Biological conclusion :

For the question of links between modalities in rows or columns, you need to know that :

-   If we have a positive scalar product, i.e. less than 90 degree difference between the arrows of two different modalities. Then we have a conjunction representing an affinity between these two modalities.

-   If we have a negative scalar product, i.e. more than 90? difference between the arrows of two different modalities. Then we have an opposition of these two modalities. They reject each other.

-   If we have a zero scalar product, i.e. the arrows between two different modalities are perpendicular. Then there is no link between these two modalities.

```{r}
fviz_ca_biplot(CA)
```

On this graph, the rows are shown in blue and the columns in red. In other words, hair color in blue and eye color in red.

From this analysis, we can see :

-   A conjunction between black, pink, brown or blond hair and blue or brown eyes.

-   A conjunction between yellow eyes and blue hair

-   An opposition between these two groups on axis 1.

In an analysis where the data are not created randomly, but come from different individuals from different populations. This could reveal a particular typology of one population in relation to another.

# PART 3: Multiple Correspondance Analysis (MCA)

## Introduction

MCA is a statistical method adapted to table of type "individuals x quatitative variable". Eingen Values correspond to means of squared correlation ratios. It could be used a pre-processing before a classification or a coinertia analysis on tables with quantitative.

## Mathematics

The MCA stands out from other multivariate analyzes with a single table (PCA, CA or mix between MCA and PCA, etc.) because it ***only takes qualitative variables*** as input. Its goal is to find the relationships between modalities by visualizing the possible associations and producing a quantitative indicator of their relationships.

The data taken by an MCA is a table comprising $J$ qualitative variables to describe $I$ statistical individuals. All individuals have the same weight $\frac{1}{I}$. The value $V_{i,j}$ in the table corresponds to the modality taken by individual $i$ for variable $j$.

To obtain a ***complete disjunctive table (CDT)*** and carry out the statistical analysis, it is necessary to subdivide each variable $j$ into $n_{j}$ modalities corresponding to all of the responses $V_{.,j}$ of the $I$ individuals for this variable. The table will always have $I$ rows but $\sum_{j=1}^J {n_j}$ columns (we will denote K columns), in which each of the $V_{i,k}$ values is either 1 or 0 depending on the presence or absence of the modality.

The value $V_{i,k}$ is transformed again : the old value of $V_{i,k}$ (1 or 0) is divided by the probability $p_{k}$ of having this modality (i.e. $p_k=\frac{\sum_{i=0}^I V_{i,k}}{I}$, result from which we subtract 1 to center the result.

$$V'_{i,k}=\frac{\sum_{i=1}^I V_{i,k}}{I}-1$$

We define the ***distance D*** between 2 individuals $i$ and $i’$ as : $$D=\frac{1}{J}\sum_{k=1}^K \frac{1}{p_{k}}(V'_{i,k}-V'_{i',k})^2$$

So 2 individuals taking the all same modalities are at a distance D = 0, the more similar the responses are over a large number of modalities, the lower the distance D will be and conversely. If two individuals share a rare modality, the distance will be reduced to take into account the common specificity.

The inertia of the point cloud : $N = \frac{K}{J} – 1$. $η^2(x,y)$ is the correlation ratio between 2 variables $x$ and $y$ and $η^2(F_s,k)$ is the ***correlation ratio*** between the modality $k$ and l'axe $F_s$. The axis $F_1$ is the one which maximize its correlation ratio with all the modalities $k$ : $$F_1=\max(\sum_{k=1}^K η^2(F_s,V_{.,k}))$$

## Interpretation

This part is a simple example of MCA using R. You can import R package using the code. Let's have a look to fictive the data set we will be working on :

```{r}
library(tidyverse)
hippo <- read.table(file = "https://github.com/raw/AnsaldiL/MODE_reproduciblescience/master/hp.csv", sep=';', header=TRUE)
hippo=hippo[,-1]
str(hippo)
```

Water consumption behaviour of Hippopotamus was observed in Penjari National Park, Benin. Drinking frequency was evaluated by a technician and rated "rarely" or "regularly". Sex is indicated with F for female and M for male. As there are 3 lakes in the park which are noted "G1", "G2" and "G3".

We are performing MCA with the package FactoMineR.

```{r, include=FALSE}
library(FactoMineR)
res.mca = MCA(hippo, quanti.sup=1)
```

Let's have a look to Eigen Values :

```{r}
res.mca$eig
fviz_eig(res.mca, addlabels = TRUE)
```

The total variance of the data set is divided between four dimensions. The first dimension concentrates 32.6% of the total variance. The first plan gathers 58% of the variance. You can better visualize the repartition of variance with the plot.

Here, you can see the results for the individuals, their coordinates ($coord$), absolute contributions ($contrib$), and the quality of their projection ($cos2$).

```{r}
head(res.mca$ind$coord)
head(res.mca$ind$cos2)
head(res.mca$ind$contrib)
```

Coordinates of the individuals are their position on the first plan. $cos2$ represents the quality of the representation of the individuals on the first plan. Contribution is how the point contributes to the creation of different axis.

You can access the same information for the variable (instead of the individuals):

```{r}
res.mca$var$coord
res.mca$var$contrib
res.mca$var$cos2
```

Let's have a look to the plot of this MCA analysis :

```{r}
plot(res.mca)
```

As this is a factice data set, some individuals are overlapping. We will not focus on this artefact of the data set.

**How to interprete these graphs?**

***General description*** : The individuals are in black. The variables are in red. The percentage of variance of the first two axis is written on them. Individuals in the center of the cloud are individuals taking a mean value for all of their caracteristics, unlikely individuals far from the middle which are specific individuals, very different from others. Close individuals present close characteristics.

***Interpretation***: Here, we can see that the first axis separates individuals drinking regularly from individuals drinking rarely. The second axis separates individuals drinking at G2 from individuals drinking at G3. Male and female seem to be separated by the first axis. As we can gather individuals sharing close properties, females seems to drink regularly and male more rarely. It is a bit less clear but male are drinking preferably in pound 1.

Note that supplementary quantitaive data could be added to the analysis to help interpretation. As they are quantitative, they would not participate to the creation of axis but they may be ploted on the final graph.

# PART 4: Canonical Correspondence Analysis (CCA)

## Introduction

Two table analysis is used when the scientist possesses two data set with expecting the second to explain the first. The first one is called Y, the response variable and the second one is called X, the explanatory variable. For example :

| Table Y \| Table X
| a data set with temperature data \| a data set with cities caracteristics such as heigth of buildings, concreted area, number of cars... \|
| a data set with the quantity of chemicals in soil \| a data set with % of pesticide applied \|
| a data set with functional traits of an animal \| a data set with the quantities of different food given to the animal \|

RDA is used when expecting a linear response from $X$ to $Y$. Only the variables of $X$ explaining $Y$ would be kept. The canonical axis are a linear combinaison of the explicative variable ($Y$). The $Y$ table will be ordered with a PCA. Then, a multiple regression of $Y$ on $X$ is done and stored in a new table, $Y'$. The residual matrix correspond to $Y-Y'$, what was not explained by $X$.

## Interpretation

Here are the library you will need to perfor the RDA, make sure your import them at the begining of your working document.

```{r, include=FALSE}
require(ade4)
require(vegan)
require(tidyverse)
```

Let's import the data, they are from @thilenius1963synecology. However, the data of abundance have been changed to facilitated the interpretation. For access to real data and metadata, please see the JDBakker Git page: [appliedmultivariatestatistics/Oak_data_47x216.csv at main · jon-bakker/appliedmultivariatestatistics (github.com)](https://github.com/jon-bakker/appliedmultivariatestatistics/blob/main/Oak_data_47x216.csv)

```{r}
Oak <- read.csv("https://github.com/raw/AnsaldiL/MODE_reproduciblescience/master/Oak_data.csv", sep=';', header = TRUE) 

Topo <- Oak %>%
select(LatAppx, LongAppx, Elev.m, Slope.deg) %>%
decostand("range")

Oak_sp=Oak[,28:83]
```

We can explore the tables

```{r}
head(Oak_sp)
head(Topo)
```

So here, the first table ("Oak_sp") represents a table of plant abundance for 47 sites (47 rows). We have 56 differents species. The second table ("Topo") corresponds to 4 topographics variables (Latitute, Longitute of sites, the elevation and the slope). As it is not the same unit for each variable we have centred-reduced these four variables.

As explained above, the CCA objective is to link the abundance of plant by the topographic variables. Here, we have *a priori* hypothesis on the effect of these environmental variables on the plants.

We have seen in introduction of the part 3 that CCA was used when the respond have a non-linear distribution but rather a gaussian relationship.

So we begin by performing a CA on the species table:

```{r}
caoak <- dudi.coa(Oak_sp,scannf=F,nf=2)
```

Then, we can perform the PCA on the CA that corresponds to the CCA with the **ADE4** package:

```{r}
ccaoak <- pcaiv(caoak,Topo,scannf=F,nf=2)
```

You can also realize the CCA with one code row only with the **vegan** package:

```{r}
ccavegan=cca(Oak_sp,Topo,scan=F)
```

It is important to use the two methods because each allows us to have different information during the analyse of the result

And we already have completed our analyse, it was no difficult! But now, we need to analyse the result and test the significant.

To test the significant, we use a bootstrapping test thank you the **ADE4** package.It is important to test if the result of the RLQ is not only due to random combination of values but that we have a real correlation between are different tables. To produce this, we perform a permutation test with the function *randtest*.

```{r}
randtest(ccaoak)#with ADE4
plot(randtest(ccaoak))
```

The outputs above correspond to the permutation test. We see that the number of permutations of columns and rows was to 999 (default value).

The results of this test show that our result is significatively different to the permutation result with a threshold of 5% (p-value = 0.02)

Finally, the results of the CCA are plot in the distribution law calculated by the permutation. So we can conclude that our environmental variables explain a part of the plant distribution.

We have now to look the collinearity between the environmental variables. It corresponds to the correlation between them. For that, we use the Variance Inflation Factor (VIF).

```{r}
test=vif.cca(ccavegan) 
test
```

If all the variables have a number lower than 10, you can conclude that you do not not have collinearity. If you have, you need to select some variables to remove with, for example, an *ordistep* function.

Now that we have tested all the conditions, we can look at the results.

```{r}
ccaoak
print(paste("The pourcentage of variance explained by topographic variables is: ",sum(ccaoak$eig)/sum(caoak$eig),"%"))
```

The first output (by **ADE4**) allows us to see each part of the analyse that could be use to describe the result. We will use especially the last part afterward. You also have, at the beginning, the number of rank that corresponds to the number of eigenvalue with their score below. We can calculate the percentage of variance explained by topographic variables.

```{r}
ccavegan
```

Here, it is the output with the **vegan** package. We can see the total inertia in the CCA analyse (5.55) and the inertia explained by the topographic variables (0.71) and by the residuals (4.85). We can therefore conclude that our topographic variables explain only 13% of the plant variances (same result that before). Below that, we have the eigenvalues for the constrained axes (variance part explained by topographic variables) and for the unconstrained axes (variance part explained by the residuals)

We can now, observe the complet plot of the CCA result:

```{r}
plot(ccaoak)
```

Here, we can have a complete plot with the **ADE4** package.

-   The "Loadings" part represents the canonical coefficients. For each axis, the arrows explain the relative weight of the topographic variables in the multiple regression calculation.

-   The "Correlation" part shows the correlation between topographic variables themselves and with CCA axes

-   The "Scores and Predictions" part allows us to understand, for each site, the real plant abundance (base of arrows) and the abundance predicted by multiple regressions with the topographic variables.

-   The "Species" part corresponds to the coordinates of each plant species in the CCA analyse

But this can also be decomposed in different plots and output.We can detail each part and observe the biological results.

We can begin by the absolute contribution of species:

```{r}
contrib=inertia.dudi(ccaoak, col = TRUE, row = FALSE)
contrib$col.abs
```

The output, here, is the absolute contribution of each species for the axis 1 and the axis 2. This allows us to observe what species contributes the most of the axes.

```{r}
ccaoak$fa
s.corcircle(ccaoak$fa)
```

This line gives you the canonical coefficients on each axis. The circle is explain in the "Loadings" part of the main plot of the RDA. On axis 1, the variable with the main contribution of the predictive power is $LatAppx$. On axis 2, $Slope.deg$ and $Elev.m$ contributes the most of the predictive power.

```{r}
ccaoak$co
s.label(ccaoak$co,boxes=F)
```

$co$ gives you access to the coordinates of each species on the first plan. The table explicits the coordinates and the plot help you to find close species. On the first axis, two species have high negative value : Adbi (-1.97) & Abgr.t (-1.50). They are opposed at species as Coco.t (0.86) & Cato (0.96). You can realize the same for the second axis.

Here, we can highlight that there is a group formed by Agse (coord = 0.31 ; 1.90 respectively for the axis 1 and 2), Agha (coord = -0.38 ; 1.25) Avfa (coord = 0.33 ; 0.88) and ALL (coord = 0.20 ; 1.14). An other on with two outliers, Abgr.t and Adbi. These groups will be correlated with the environment data further in the analysis.

```{r}
ccaoak$cor 
s.corcircle(ccaoak$cor)
```

Here, we have the correlation between axes and the topographic variables. We can conclude that the elevation (corr = -0.15\[axis1\] ; -0.45\[axis2\]) and the slope (corr = 0.12\[axis1\] ; 0.65\[axis2\]) are strongly negatively correlated. The Longitude seems to be positively correlated with the slope and negatively correlated with the elevation.

```{r}
s.match(ccaoak$ls, ccaoak$li)
plot(ccavegan,scaling=1) 
plot(ccavegan,scaling=2) 
```

The first plot with the arrows represents the real abundances for each site (base of the arrow) and the predictions of abundance calculated by a multiple regression thanks to the topographic variables. We see here that there are some large differences between reality and prediction. Most of the sites are in a center group but some sites differ. It is the case with the sites n°7,40,41,42. We have another group with the sites n°43,44,45,46.

The second plot, with *scaling = 1* preserves euclidian distances between stations.\
The third plot, with *scaling = 2* preserves correlations between species.

To conclude on the CCA analyse, we have seen that the topographic variables explain 13% of the plant variances. The first axis represents 34% of the constrained inertia and is explained by the Latitude (corr = 0.99). We can see that two species have a high coordinate in this axis (Adbi & Abgr.t).\
The second axis (18% of the constrained inertia) is explained by the Elevation (corr = -0.45), the Slope (corr = 0.65) and the Longitude (corr = 0.33).

Finally, the sites n°43,44,45,46 are sites with a small latitude value and with a high abundance of species as Abgr.t.

# PART 5: RLQ

## Introduction

The RLQ analysis is a generalized co-inertia with 3 tables : environmental variables for each sites (R), species abundances for each sites (L) and trait variables for each species (Q). We will have to perform a CA on the table L and weighted analysis (MCA or PCA) on tablesR and Q. Then a function named *rlq()* perform 3 co-inertia analysis between R and L and L and Q and then between the results of each co-inertia.

## Interpretation

The use of RLQ analysis is important in ecology to integrate the traits of species with the environmental variables. So here, we don't have 2 tables (environment & specie) as RDA part but 3 tables:

-   environmental variables by sites (R)

-   abundance of species by sites (L)

-   trait values by species (Q)

To perform the RLQ, we need to decompose the analyse by three types of analyses already done in this chapter:

-   we will use a CA analyse on the abundance of species

-   we will use a MCA on the environmental table by taking the sites weight on the CA

-   we will use a MCA on the trait table by taking the species weight on the CA

Here, we use MCA for R and Q because our variables are factors but you can perform a PCA if your variables are quantitative. Warning, for R and Q you have the obligation to weight by the L table (see below).

You can import R package using the code

```{r}
library(tidyverse)
library(ade4)
library(vegan)
library(ggplot2)
library(factoextra)
library(corrplot)
library(RVAideMemoire)
library(PerformanceAnalytics)
```

Here, we work with a dataset of **ADE4** package. The data comes from @dray2008testing, @legendre1997relating

```{r}
#import the dataset
data(aviurba)

#create the three tables
summary(aviurba$mil)    #(R)
R<-aviurba$mil

summary(aviurba$fau)    #(L)
L<-aviurba$fau

summary(aviurba$traits) #(Q)
Q<-aviurba$traits
```

and explore the tables

```{r}
head(R)  
head(L)
head(Q)
```

The first part is to perform our CA on specie table

```{r}
afcL <- dudi.coa(log(L+1), scannf = FALSE) 
afcL  

```

The first CA is done. We use log transformation because the abundance of species has a large range and we add "+1" to avoid the log(0) for some species. You must adapt the presence of transformation (or not) to your data.

Now, we can perform the two MCA analysis on the trait table and environmental table

```{r}
acmR <- dudi.acm(R, row.w = afcL$lw, scannf = FALSE,nf = 4)
scatter(acmR)

acmQ <- dudi.acm(Q, row.w = afcL$cw, scannf = FALSE,nf = 4)
scatter(acmQ)

```

The scatterplot allows to see the ordination of each table and the repartition of factor on the simple axe of MCA.

But now, we will use the RLQ analyse that creates two co-inertia (R-L, L-Q), assembles and compares the co-inertia. We use the rlq function for that.

```{r}
rlq <- rlq(acmR, afcL, acmQ, scannf = FALSE)
rlq
axe=c(1:8)
print(paste("The contribution of axe n°",axe, "are", rlq$eig/(sum(rlq$eig))*100,"%"))
#randtest(rlq)
#summary(rlq)
#plot(rlq)

```

Here, the output of the RLQ is complex but only few information are, at this point important. We see that we have 8 eigenvalues and we have their values. All the different compounds of the output will be used after in representations or analysis.

Nevertheless, we can calculate the contribution of each axis of the RLQ by performing the formule below: $$\frac{rlq$eig}{\sum{rlq$eig}}.100$$

After that, and before to plot and analyse the result, it is important to test if the result of the RLQ is not only due to random combination of values but that we have a real correlation between the different tables. To produce this, we perform a permutation test with the function *randtest*.

```{r}
randtest(rlq)
plot(randtest(rlq))


```

The outputs above correspond to the permutation test. We see that the number of permutations of columns and rows was to 999 (default value).

The results of this test shows that the permutation of sites (rows) is the first result (p_value=0.1%) and the permutation of species (columns) is the second result (p_value=1.9%) So each result is significant (5% threshold) and we can conclude that our RLQ result is not linked to random effect.

Finally, the results of the RLQ are plot in the distribution law calculated by the permutation.

We can know analyse the result by using different types of plot:

```{r}
plot(rlq)
```

The plot above is the complet plot, difficult to understand and that will be resumed step by step afterward.

However, we have here, up the score (the position) of the sites and species in this analyse. Below in the centre, we have the correlation of the environmental variables between them (R Canonical weights) and the correlation of the traits between them (Q Canonical weight). The plots of the R axes and the Q axes represent the reprojection of axis in the RLQ analyse in order to the simple analyses (CA or MCA here).

Now, we can decompose the result and start with the environmental variables:

```{r}
#analyse the environmental variables:
round(rlq$l1,dig=2)
s.label(rlq$l1, boxes=FALSE)

#calcul of the absolute contribution
iner=inertia.dudi(rlq,col.inertia=T,row.inertia=T)

abscoiE=iner$row.abs
abscoiE

```

The table corresponds to the position of each modality for each environmental variables on the two first axis. You can link this table to the plot. For example, "veg.cover.R100" seems to stand out from the other modalities. In the table, we can see that the position on the first axis is 4.35 and -4.95 on the second axis. Closer to the center of the cloud, we can talk about "veg.cover.R98" that have position of 1.74 on the first axis and -1.64 on the second axis. In addition, we can observe the contribution of this two variables. On the second table we have the absolute contribution of the modalities for each variable. The contribution of the "veg.cover.R100" of the first axis is 17% and 22% on the second axis. This is a high value if we take the threshold to 1/N with N the number of modalities (28 here so 3%). The contribution of "veg.cover.R98" is 4.6% for the first axis and 4.1% for the second axis.

We can therefore conclude that these two variables are correlated. This seems to be biologically coherent since these two variables are part of a plant cover gradient and are the two highest values.

We can also see that the variable negatively correlated with the first axis are modalities explaining the urbanization ()

Now we can realise the same analyse for the traits.

```{r}
#analyse the traits:
round(rlq$c1,dig=2)
s.arrow(rlq$c1, boxes=FALSE)

#absolute contribution
abscoiV=iner$col.abs
abscoiV

```

As seen before, the first table is the position of the modalities for the two first axes of the RLQ analyses. The plot is the representation of that and the second table is the absolute contributions.

Here, we have larger differences between the different variable and modalities. We're going to focus on a few examples, but you can take the time to go into detail about each modality. "feed.strat.foliage" (coord= 2.60; -1.97 for first and second axis respectively) and "breeding.scrub" (coord= 1.71: -2.22) break away from the centre of the cloud. In addition their contributions are all greater than 10%. We can conclude that these modalities are correlated and highly contribute of axes.

The last point is to observe the species coordinates.

```{r}
#analyse the species:
round(rlq$lQ,dig=2)
s.label(rlq$lQ,boxes=FALSE)
s.label(rlq$lQ, label=aviurba$species.names.fr,boxes=FALSE)
```

Here, is to see that we seem to have a gradient of response along the first axis (i.e. Loriot Jaune/Sp37 \[coord=1.18; -1.33\] and Alouette des champs/Sp9 \[coord=0.99; 1.17\]) and along the second axis (i.e. Loriot Jaune/Sp37 \[coord=1.18; -1.33\] and Fauvette des jardins/Sp21 \[coord=0.99; 1.17\]). We will see that we need to link this to the other variables (trait and environment) For species and sites, we do not have contribution because these variables do not contribute of the axis creation.

To conclude on this analysis, the two first axes explain 86.2% of the total variation with respectively 66.7%, and 19.5% of the total inertia.

The correlations between the environmental variables and the RLQ axes showed that the first axis is explained by Veg.cover.R100 (contribution of 16.95%), veg.cover.R22 (7.43%), noisy.yes (8.82%), small.bui.yes (11.22%), field.yes (8.17%), field.no (8.2%) and Veg.cover.R22 (7.43%). The traits variables that contribute to the first axis of the RLQ axes are breeding.ground (28.19%), feed.strat.foliage (19.8%), breeding.scrub (13.3%) and negatively with breeding.building (17.18%) and feed.strat.aerial (7.09%). In other words, breeding.ground, feed.strat.foliage, breeding.scrub, feed.strat.aerial and breeding.building were the higher explanatory attributes for this RLQ axis (explained by the length and the angle between the axes and the vectors).

# To go further with Multivariate Analysis

Two other methods of two tables analysis were not presented here :

-   **Co-inertia** is used when you have a table with you response variable $Y$ and a second one with *potentially* explaining variables, $X$, **but you do not have any *a priori*** about the effect of $X$ on $Y$.

-   **CCA** is used with two table, where both $Y$ **and** $X$ **are quantitative data**. It works the same as RDA.

-   PCCA is CCA like but it includes a spacial dependency... Spacial and Temporal dependencies are a great deal in data analysis !

Even if this presentation of multivariate analysis was not exhaustive, you can have a better view of what you can do when you are to analyse a data table. The method chosen will depend on the number of tables and the type of the variables (quantitative or qualitative). To help you decide what analysis to use regarding to the ones we presented here, you can have a look to this decision tree !

![**Decision tree multivariate analysis (from A. Pannard)**](decision_tree.jpg)