-
Notifications
You must be signed in to change notification settings - Fork 8
/
AMV_chapter.qmd
1066 lines (656 loc) · 51.5 KB
/
AMV_chapter.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
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).