ProcMod.Rmd 9.51 KB
Newer Older
1
---
Eric Coissac committed
2 3
title: "ProcMod"
author: "Christelle Melodelima & Eric Coissac"
4 5
date: "`r Sys.Date()`"
output: 
Eric Coissac committed
6 7
  rmarkdown::html_vignette: default
  pdf_document: default
8
vignette: >
Eric Coissac committed
9
  %\VignetteIndexEntry{ProcMod}
10
  %\VignetteEncoding{UTF-8}
Eric Coissac committed
11
  %\VignetteEngine{knitr::rmarkdown}
12 13 14
---

```{r setup, include = FALSE}
Eric Coissac committed
15
knitr::opts_chunk$set(
16
  collapse = TRUE,
Eric Coissac committed
17
  comment = "#>"
18 19 20
)
```

Eric Coissac committed
21
## Aims of the module
22

Eric Coissac committed
23
Expliquer un tableau multavariées decrivant des individus ou des sites (ci dessous dénommés individus) par un ensemble de tableaux eux-même multivariés décrivant les mêmes individus. Par exemple, expliquer les changement de communauté d'espèces entre différents sites géographiques à partir de tableaux de données climatiques, chimique, d'espèces d'autres groupes taxinomiques... Chaque tableau est considéré comme une variable explicative sans chercher à donner un rôle à chacune des variables qui le composent.
24

Eric Coissac committed
25
## Model principes
26

Eric Coissac committed
27
L'idée est de s'appuyer sur les analyses procustéennes en les généralisant à plusieurs ($k$) tableaux. Pour mémoire, l'analyse procustéenne consiste à superposer deux nuages de points dans un espace de dimensions quelconque en réalisant trois operations: 
28

Eric Coissac committed
29 30 31
1. une translation (centrage des données) 
1. une rotation 
1. une mise à l'échelle. 
32

Eric Coissac committed
33
Dans notre cas, nous considérerons que les deux premières opérations ont pour seul but de projeter l'ensemble des tableaux réponse et explicatifs, dans un espace commun. La troisième operation d'homothétie servira de base à l'analyse de partition de la variance du tableau réponse. Cette approche a de forts liens avec les analyses de co-inerties développées par Chessel et al 
34

Eric Coissac committed
35
## Données en entrée
36

Eric Coissac committed
37
Les tableaux utilisés dans cette analyse doivent se projeter dans un espace orthogonal, ce qui implique aucune corrélation entre les colonnes des tableaux. Ils doivent tous décrire les mêmes individus et doivent donc tous avoir le même nombre de lignes $n$. 
38

Eric Coissac committed
39
### Tableaux de variables quelconques
40

Eric Coissac committed
41
Comme posé en préambule, l'analyse vise à mesurer l'effet global de $k$ tableaux de variables sur un tableau réponse sans s'interresser à l'effet individuel de chacune des variables des différents tableaux explicatifs. Les tableaux utilisés peuvent donc sans perte d'information être projetés dans un espace orthogonal par une simple PCA.
42

Eric Coissac committed
43
### Tableaux de distances
44

Eric Coissac committed
45
Il est possible de caractériser la dissimilarité entre les individus par une autre mesure que la corrélation (*cf* PCA ci-dessus) en estimant un tableau de distances par une méthode appropriée au type de données étudiées. Par exemple une distance de Bray-Curtis ou de Jaccard, lorsqu'il s'agit de comparer les communautés d'espèces présentes dans plusieurs sites.
46

Eric Coissac committed
47
Si la distance utilisée est une métrique, le tableau de distance pourra être projeter dans un espace orthogonal à $n-1$ dimensions par une PCoA.
48

Eric Coissac committed
49
Si la distance utilisée n'est pas une métrique, il faudra recourrir à une méthode non paramétrique telque la NMDS ou alterer le tableau de distance pour le rendre diagonalisable.
50

Eric Coissac committed
51
## Méthode de calculs du modèle procustéen
52

Eric Coissac committed
53
L'analyse procustéenne peut être assimilée à un modèle linéaire entre deux tableaux:
54

Eric Coissac committed
55 56
- Le tableau réponse : $Y$
- Le tableau explicatif : $X$ 
57

Eric Coissac committed
58
Si l'ensemble des variables de chaque tableau appartiennent à $\mathbb{R}$:
59

Eric Coissac committed
60 61
- L'opération de translation est réalisée via le centrage des variables des tableaux $X$ et $Y$. Dorénavant, la notation $X$ et $Y$ représentera les tableau après centrage.
- La rotation pour projeter le tableau $X$ sur $Y$ est calculer à partir de la décomposition en valeurs singulières de la matrice $Y'X$ soit $Y'X = U\Lambda V'$. On définit la rotation de $X$ sur $Y$ de la manière suivante:
62

Eric Coissac committed
63
$$ Rot(X|Y) = XVU' $$
64

Eric Coissac committed
65
- Le facteur d'homothétie $a$ se calcule selon :
66

Eric Coissac committed
67
$$ a = \frac{\sum diag(\Lambda)}{\sum diag(X'X)}$$
68

Eric Coissac committed
69
$\sum diag(\Lambda)$ est la co-inertie entre les matrices $X$ et $Y$ et $\sum diag(X'X)$ est l'inertie du tableau $X$. Co-inertie et inertie pouvant être assimilées à la covariance et à la variance de vecteurs. Le facteur $a$ est donc l'équivalent en dimention 1 de la pente de la droite de régression linéaire de $Y$ par $X$:
Christelle Melodelima committed
70 71

$$
Eric Coissac committed
72
a=\frac{Cov_{XY}}{Var_X}
Christelle Melodelima committed
73 74
$$

Eric Coissac committed
75
Nous proposons ici de généraliser l'analyse procustéen à $k$ tableaux en résolvant la régression multiple du tableau de réponse $Y$ par $k$ matrices explicatives $X_1,\,X_2,\,...\,X_k$.
Christelle Melodelima committed
76

Eric Coissac committed
77
Le calcul des coefficients d'échelles $a_i$ est réalisé par la même approche que pour ceux d'une régression linéaire multiple, mais ici à partir de la matrice d'inertie et de co-inertie des tableaux $Y$ et $X_i$:
Christelle Melodelima committed
78

Eric Coissac committed
79
### Calcul de l'inertie d'une matrice et de co-inertie entre deux matrices
Christelle Melodelima committed
80

Eric Coissac committed
81
Soit $CoI_{XY}$ la cointertie entre les tableaux $X$ et $Y$ et $I_X$ l'inertie du tableau $X$. $CoI_{XY}$ se calcule à partir de la décomposition en valeur singulière de la matrice $Y'X = U\Lambda V'$.
Christelle Melodelima committed
82

Eric Coissac committed
83 84 85
$$
CoI_{XY} = \frac{\sum{diag(\Lambda)}}{n-1}
$$
Christelle Melodelima committed
86

Eric Coissac committed
87
Comme pour la Covariance, $CoI_{XY}=CoI{YX}$ et $CoI_{XX}=I_X$
Christelle Melodelima committed
88

Eric Coissac committed
89
### Calcul du coefficient de corrélation entre deux matrices
Christelle Melodelima committed
90

Eric Coissac committed
91
Par analogie on définit $R_{XY}$ comme le coefficient de corrélation entre deux matrices $X$ et $Y$ de la manière suivante :
Christelle Melodelima committed
92

Eric Coissac committed
93 94 95
$$
R_{XY} = \frac{CoI_{XY}}{\sqrt{I_{X}I_{Y}}}
$$
Christelle Melodelima committed
96

Eric Coissac committed
97
### Calcul des facteurs d'échelle en procuste multiple à $k$ tableaux
Christelle Melodelima committed
98

Eric Coissac committed
99
On note à partir de maintenant :
Christelle Melodelima committed
100

Eric Coissac committed
101 102 103 104 105 106
- $X = \{X_1,\,X_2,\,...\,X_k\}$ l'ensemble des $k$ tableaux explicatifs centrés
- $Y$ le tableau réponse centré
- $M \in \{Y\} \cup X$ tel que $M$ est la matrice de plus grande dimension.
- $CoI_{YX}$ la matrice colonne des co-inerties entre $Y$ et chacun des éléments de $X$
- $CoI_{XX}$ la matrice d'inertie, co-inerties entre tous les éléments de $X$
- $a = \{a_1\,a_2,\,...\,a_k\}$ l'ensemble des coefficients d'échelle associés à chacun des éléments de $X$ dans le modèle procustéen multiple.
Christelle Melodelima committed
107

Eric Coissac committed
108
$a$ se calcule par :
Christelle Melodelima committed
109

Eric Coissac committed
110 111
$$
  a = (CoI_{XX})^{-1} CoI_{YX}
Christelle Melodelima committed
112 113
$$

Eric Coissac committed
114
Les prédictions du modèle peuvent donc s'écrirent:
Christelle Melodelima committed
115

116
$$
Eric Coissac committed
117
\widehat{Rot(Y|M)} = a_1 Rot(X_1|M) + a_2 Rot(X_2|M)+\,...\,+ a_k Rot(X_k|M)
Christelle Melodelima committed
118 119
$$

Eric Coissac committed
120
### Interaction entre deux tableaux explicatifs
121

Eric Coissac committed
122
L'interaction entre deux tableaux explicatifs $X_i$ et $X_j$ est estimée de manière analogue à l'interaction dans un modèle linéaire multiple. Le principe est de poser que $a_i$ le facteur d'échelle associé à $X_i$ est une fonction affine de $X_j$ et symétriquement $a_j$ est une fonction affine de $X_i$. Pour réaliser ce calcul il est nécessaire de projeter les deux matrices explicatives sur $M$.
Christelle Melodelima committed
123 124

$$
Eric Coissac committed
125 126
a_i = (b_i Rot(X_j|M) + c_i) \\
a_j = (b_j Rot(X_i|M) + c_j)
Christelle Melodelima committed
127 128
$$

Eric Coissac committed
129 130
Dans le cas de deux tableaux, le modèle procustéen expliquant $Y$ par $X_1$, $X_2$ et l'interaction de $X_1$ et $X_2$
peut donc s'écrire : 
Christelle Melodelima committed
131 132

$$
Eric Coissac committed
133 134
\begin{aligned}
\widehat{Rot(Y|M)} = & a_1 \cdot Rot(X_1|M) +  a_2 \cdot Rot(X_2|M) \\
Christelle Melodelima committed
135
\\
Eric Coissac committed
136 137
  = & (b_1 Rot(X_2|M) + c_1) \cdot Rot(X_1|M) + \\ 
& (b_2 Rot(X_1|M) + c_2) \cdot Rot(X_2|M) \\
Christelle Melodelima committed
138
\\
Eric Coissac committed
139 140
  = & c_1 Rot(X_1|M) + c_2 Rot(X_2|M) + (b_1+b_2) Rot(X_1|M) \cdot Rot(X_2|M)
\end{aligned}
Christelle Melodelima committed
141 142
$$

Eric Coissac committed
143 144
En renommant les facteurs d'échelle les prédictions s'écrivent:

Christelle Melodelima committed
145
$$
Eric Coissac committed
146
\widehat{Rot(Y|M)} =  a_1 Rot(X_1|M) + a_2 Rot(X_2|M) + a_{1,2} Rot(X_1|M) \cdot Rot(X_2|M)
Christelle Melodelima committed
147 148
$$

Eric Coissac committed
149
Ce qui revient à construire un nouveau modèle sans interaction incluant un tableau explicatif supplémentaire 
Christelle Melodelima committed
150 151

$$
Eric Coissac committed
152
X_{1,2}=Rot(X_1|M) \cdot Rot(X_2|M)
Christelle Melodelima committed
153 154
$$

Eric Coissac committed
155
Ce principe peut être généralisé à l'interaction entre plus de deux tableaux.
Christelle Melodelima committed
156

Eric Coissac committed
157
## Partition de l'inertie de $Y$ selon le modèle procustéen
Christelle Melodelima committed
158

Eric Coissac committed
159
La variation total $SCT$ est définie comme suit:
160

161
$$
Eric Coissac committed
162
SCT = I_Y * (n-1)
163
$$
Eric Coissac committed
164
La variation résiduelle non expliquée par le modèle est :
165

166
$$
Eric Coissac committed
167
SCR = \sum{(Rot(Y|M)-\widehat{Rot(Y|M)})^2}
168 169
$$

Eric Coissac committed
170
Selon l'approche de Scherrer on propose de définir $SCE_i$ la contribution de $X_i$ à la variation de $Y$ peut s'écrire :
171 172

$$
Eric Coissac committed
173
SCE_i = a_i \; R_{X_iY} \; SCT
174 175
$$

Eric Coissac committed
176
On propose de tester la signignificativité de l'effet de chacun des tableaux explicatifs par une méthode de permutation des lignes de chacune des matrices explicatives. 
Christelle Melodelima committed
177

Eric Coissac committed
178
## How to build a procrustean model
179 180

```{r}
Eric Coissac committed
181
library(ProcMod)
182 183 184
```


Eric Coissac committed
185
### Loading of the demo data
186

Eric Coissac committed
187
They consist in two MOTUs tables describing bacterial and eukaryiote communities at 21 sites accros the eastern coast of Australia.
188 189

```{r}
Eric Coissac committed
190 191
data("eukaryotes")
data("bacteria")
192 193
```

Eric Coissac committed
194
At each site environmental data are also available to describe soil chemistry, climat and site location.
195 196

```{r}
Eric Coissac committed
197 198 199
data("soil")
data("climat")
data("geography")
200 201
```

Eric Coissac committed
202
### Processing of the MOTUs data
203

Eric Coissac committed
204
Using the *vegan* package MOTUs frequencies are transformed according to Hellinger by take the square root of the MOTUs relative frequencies at each site
205 206 207 208

```{r}
library(vegan)

Eric Coissac committed
209 210
bac.hellinger = decostand(bacteria,method ="hellinger")
euk.hellinger = decostand(eukaryotes,method ="hellinger")
211 212
```

Eric Coissac committed
213
Bray Curtis distances among sites are computed according to both the communities
214 215

```{r}
Eric Coissac committed
216 217
bac.bray = vegdist(bac.hellinger,method = "bray")
euk.bray = vegdist(euk.hellinger,method = "bray")
218 219
```

Eric Coissac committed
220
### Processing the environmental data
221

Eric Coissac committed
222
Soil and climatic data are centered and rescaled for having the same influence
223

224
```{r}
Eric Coissac committed
225 226
soil.rescaled   = scale(soil,   center = TRUE, scale = TRUE)
climat.rescaled = scale(climat, center = TRUE, scale = TRUE)
227 228
```

Eric Coissac committed
229
### Assembling the data for the model
230

231
```{r}
Eric Coissac committed
232 233 234 235 236
data = procmod.frame(euk       = euk.bray,
                     bac       = bac.bray,
                     climat    = climat,
                     soil      = soil,
                     geography = geography)
237
```
238

239
```{r}
Eric Coissac committed
240
euk.pm = pm(formula = euk ~ soil + climat + geography, data = data)
241 242
euk.pm
```
243

244 245 246
```{r}
plot(euk.pm)
```
247

Eric Coissac committed
248 249

```{r}
Eric Coissac committed
250
anova(euk.pm)
Eric Coissac committed
251 252
```

253