Principes-francais.Rmd 9.52 KB
Newer Older
1 2 3 4
---
title: "ProcMod"
author: "Christelle Melodelima & Eric Coissac"
date: "`r Sys.Date()`"
5
output: 
Eric Coissac committed
6
  rmarkdown::html_vignette: default
Eric Coissac committed
7
  pdf_document: default
8 9 10
vignette: >
  %\VignetteIndexEntry{Vignette Title}
  %\VignetteEncoding{UTF-8}
11
  %\VignetteEngine{knitr::rmarkdown}
12 13 14 15 16 17 18 19 20 21 22
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Aims of the module

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 25 26

## Model principes

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

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

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

35
## Données en entrée
36

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 39 40

### Tableaux de variables quelconques

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 43 44

### Tableaux de distances

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

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

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

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

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

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

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

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 63 64

$$ Rot(X|Y) = XVU' $$

65
- Le facteur d'homothétie $a$ se calcule selon :
66 67 68

$$ a = \frac{\sum diag(\Lambda)}{\sum diag(X'X)}$$

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$:
70 71 72 73 74

$$
a=\frac{Cov_{XY}}{Var_X}
$$

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$.
76

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$:
78 79 80

### Calcul de l'inertie d'une matrice et de co-inertie entre deux matrices

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'$.
82 83 84 85 86 87 88

$$
CoI_{XY} = \frac{\sum{diag(\Lambda)}}{n-1}
$$

Comme pour la Covariance, $CoI_{XY}=CoI{YX}$ et $CoI_{XX}=I_X$

89
### Calcul du coefficient de corrélation entre deux matrices
90

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 :
92 93 94 95 96

$$
R_{XY} = \frac{CoI_{XY}}{\sqrt{I_{X}I_{Y}}}
$$

97
### Calcul des facteurs d'échelle en procuste multiple à $k$ tableaux
98

99
On note à partir de maintenant :
100

101 102
- $X = \{X_1,\,X_2,\,...\,X_k\}$ l'ensemble des $k$ tableaux explicatifs centrés
- $Y$ le tableau réponse centré
103
- $M \in \{Y\} \cup X$ tel que $M$ est la matrice de plus grande dimension.
104 105 106
- $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.
107 108 109 110 111 112 113

$a$ se calcule par :

$$
  a = (CoI_{XX})^{-1} CoI_{YX}
$$

114
Les prédictions du modèle peuvent donc s'écrirent:
115 116 117 118 119 120 121

$$
\widehat{Rot(Y|M)} = a_1 Rot(X_1|M) + a_2 Rot(X_2|M)+\,...\,+ a_k Rot(X_k|M)
$$

### Interaction entre deux tableaux explicatifs

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$.
123 124 125 126 127 128

$$
a_i = (b_i Rot(X_j|M) + c_i) \\
a_j = (b_j Rot(X_i|M) + c_j)
$$

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 : 
131 132

$$
133
\begin{aligned}
134 135 136 137 138 139
\widehat{Rot(Y|M)} = & a_1 \cdot Rot(X_1|M) +  a_2 \cdot Rot(X_2|M) \\
\\
  = & (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) \\
\\
  = & 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)
140
\end{aligned}
141 142
$$

143
En renommant les facteurs d'échelle les prédictions s'écrivent:
144 145 146 147 148

$$
\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)
$$

149
Ce qui revient à construire un nouveau modèle sans interaction incluant un tableau explicatif supplémentaire 
150 151 152 153 154

$$
X_{1,2}=Rot(X_1|M) \cdot Rot(X_2|M)
$$

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

157
## Partition de l'inertie de $Y$ selon le modèle procustéen
158

159
La variation total $SCT$ est définie comme suit:
160 161 162 163

$$
SCT = I_Y * (n-1)
$$
164
La variation résiduelle non expliquée par le modèle est :
165 166 167 168 169

$$
SCR = \sum{(Rot(Y|M)-\widehat{Rot(Y|M)})^2}
$$

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 173 174 175

$$
SCE_i = a_i \; R_{X_iY} \; SCT
$$

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
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. 

## How to build a procrustean model

```{r}
library(ProcMod)
```


### Loading of the demo data

They consist in two MOTUs tables describing bacterial and eukaryiote communities at 21 sites accros the eastern coast of Australia.

```{r}
data("eukaryotes")
data("bacteria")
```

At each site environmental data are also available to describe soil chemistry, climat and site location.

```{r}
data("soil")
data("climat")
data("geography")
```

### Processing of the MOTUs data

Using the *vegan* package MOTUs frequencies are transformed according to Hellinger by take the square root of the MOTUs relative frequencies at each site

```{r}
library(vegan)

bac.hellinger = decostand(bacteria,method ="hellinger")
euk.hellinger = decostand(eukaryotes,method ="hellinger")
```

Bray Curtis distances among sites are computed according to both the communities

```{r}
bac.bray = vegdist(bac.hellinger,method = "bray")
euk.bray = vegdist(euk.hellinger,method = "bray")
```

### Processing the environmental data

Soil and climatic data are centered and rescaled for having the same influence

```{r}
soil.rescaled   = scale(soil,   center = TRUE, scale = TRUE)
climat.rescaled = scale(climat, center = TRUE, scale = TRUE)
```

### Assembling the data for the model

```{r}
data = procmod.frame(euk       = euk.bray,
                     bac       = bac.bray,
                     climat    = climat,
                     soil      = soil,
                     geography = geography)
```

```{r}
euk.pm = pm(formula = euk ~ soil + climat + geography, data = data)
euk.pm
```

```{r}
plot(euk.pm)
```


```{r}
anova(euk.pm)
```