Commit 026d9e01 by Eric Coissac

Change vignette name

parent 41cb691a
ProcMod.html
ProcMod.R
ProcMod_files
Principes-english.R
---
title: "ProcMod"
author: "Christelle Melodelima & Eric Coissac"
date: "`r Sys.Date()`"
output:
rmarkdown::html_vignette: default
pdf_document: default
vignette: >
%\VignetteIndexEntry{Vignette Title}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
## Aims of the module
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.
## Model principes
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:
1. une translation (centrage des données)
1. une rotation
1. une mise à l'échelle.
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
## Données en entrée
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$.
### Tableaux de variables quelconques
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.
### Tableaux de distances
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.
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.
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.
## Méthode de calculs du modèle procustéen
L'analyse procustéenne peut être assimilée à un modèle linéaire entre deux tableaux:
- Le tableau réponse : $Y$
- Le tableau explicatif : $X$
Si l'ensemble des variables de chaque tableau appartiennent à $\mathbb{R}$:
- 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:
$$ Rot(X|Y) = XVU' $$
- Le facteur d'homothétie $a$ se calcule selon :
$$ a = \frac{\sum diag(\Lambda)}{\sum diag(X'X)}$$
$\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$:
$$
a=\frac{Cov_{XY}}{Var_X}
$$
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$.
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$:
### Calcul de l'inertie d'une matrice et de co-inertie entre deux matrices
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'$.
$$
CoI_{XY} = \frac{\sum{diag(\Lambda)}}{n-1}
$$
Comme pour la Covariance, $CoI_{XY}=CoI{YX}$ et $CoI_{XX}=I_X$
### Calcul du coefficient de corrélation entre deux matrices
Par analogie on définit $R_{XY}$ comme le coefficient de corrélation entre deux matrices $X$ et $Y$ de la manière suivante :
$$
R_{XY} = \frac{CoI_{XY}}{\sqrt{I_{X}I_{Y}}}
$$
### Calcul des facteurs d'échelle en procuste multiple à $k$ tableaux
On note à partir de maintenant :
- $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.
$a$ se calcule par :
$$
a = (CoI_{XX})^{-1} CoI_{YX}
$$
Les prédictions du modèle peuvent donc s'écrirent:
$$
\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
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$.
$$
a_i = (b_i Rot(X_j|M) + c_i) \\
a_j = (b_j Rot(X_i|M) + c_j)
$$
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 :
$$
\begin{aligned}
\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)
\end{aligned}
$$
En renommant les facteurs d'échelle les prédictions s'écrivent:
$$
\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)
$$
Ce qui revient à construire un nouveau modèle sans interaction incluant un tableau explicatif supplémentaire
$$
X_{1,2}=Rot(X_1|M) \cdot Rot(X_2|M)
$$
Ce principe peut être généralisé à l'interaction entre plus de deux tableaux.
## Partition de l'inertie de $Y$ selon le modèle procustéen
La variation total $SCT$ est définie comme suit:
$$
SCT = I_Y * (n-1)
$$
La variation résiduelle non expliquée par le modèle est :
$$
SCR = \sum{(Rot(Y|M)-\widehat{Rot(Y|M)})^2}
$$
Selon l'approche de Scherrer on propose de définir $SCE_i$ la contribution de $X_i$ à la variation de $Y$ peut s'écrire :
$$
SCE_i = a_i \; R_{X_iY} \; SCT
$$
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)
```
---
title: "ProcMod Package"
author: "Eric Coissac & Christelle Gonindard-Melodelima"
title: "ProcMod"
author: "Christelle Melodelima & Eric Coissac"
date: "`r Sys.Date()`"
output:
rmarkdown::html_vignette
rmarkdown::html_vignette: default
pdf_document: default
vignette: >
%\VignetteIndexEntry{Vignette Title}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{ProcMod}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
---
```{r setup, include = FALSE}
library(ProcMod)
knitr::opts_knit$set(
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
root.dir = system.file("extdata", package="ProcMod")
comment = "#>"
)
```
## TODO
## Aims of the module
- interactions
On fait la rotation des tableau sur Y puis on multiplie terme à terme les
deux tableau après rotation et on introduit le nouveau tableau comme une
variable suplementaire
- VIF selection tableau non corrélés
- pmm proc model mixte avec les facteurs
- virer l'intersept
- partitioner un tableau explicatif en sous-tableau
- connecter avec glasso
- connecter avec path analysis (lavaan)
- test randomisation anova
- regarder les functions : shuffle, how
- test rendomisation coeficients
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.
## interrogation :
## Model principes
Quid des pentes négative alors que les cov sont > 0
--> Forcer à 0 les pentes négatives avec un warning
pour prévenir l'utilisateur
## Idéee de présentation du papier
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:
-> partir du procruste, proposer de lénéralisé à K tableau et montrer que c'est un cas particulier de l'analyse de coinertie. Apport supplementaire avec la partition de variance....
-> tableau de var/cov donc de corrélation --> corrélation partiel --> réseaux
1. une translation (centrage des données)
1. une rotation
1. une mise à l'échelle.
## Mathematic principles
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
## Estimation des coefficients lors d'une régression linéaire multiple (2 variables)
## Données en entrée
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$.
$$ \hat{y}=a_1 x_1 + a_2 x_2 + b$$
### Tableaux de variables quelconques
On cherche à minimiser l'erreur commise entre la prédiction et la vraie valeur (méthode des moindres carrés)
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.
$$
\left\{
\begin{array}{rcr}
\sum_i{(y_i- \hat{y_i})^2}=0 \\
\sum_i{(y_i- (a_1 x_{1,i} + a_2 x_{2,i} + b))^2}=0
\end{array}
\right.
$$
Pour trouver les paramètres qui minimise cette fonction, on calcule les trois dérivées partielles et on les annule :
### Tableaux de distances
$$
\left\{
\begin{array}{rcr}
-2 \sum_i{x_{1,i}(y_i- a_1 x_{1,i} + a_2 x_{2,i} + b)}=0 \\
-2 \sum_i{x_{2,i}(y_i- a_1 x_{1,i} + a_2 x_{2,i} + b)}=0 \\
-2 \sum_i{(y_i- a_1 x_{1,i} + a_2 x_{2,i} + b)}=0
\end{array}
\right.
$$
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.
A partir de la dernière équation, on obtient :
$$
\overline{y}-a_1\overline{x_1}-a_2\overline{x_2}-b=0
\\
\Leftrightarrow \hat{b}=\overline{y}-a_1\overline{x_1}-a_2\overline{x_2}
$$
En remplaçant dans la première équation on obtient :
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.
$$
\frac{\sum_i{x_{1,i}y_i}}{n}-a_1\frac{\sum_i{x_{1,i}^2}}{n}-a_2\frac{\sum_i{x_{1,i}x_{2,i}}}{n}-b\frac{\sum_i{x_{1,i}}}{n}=0
\\
\Leftrightarrow
\frac{\sum_i{x_{1,i}y_i}}{n}-a_1\frac{\sum_i{x_{1,i}^2}}{n}-a_2\frac{\sum_i{x_{1,i}x_{2,i}}}{n}-(\overline{y}-a_1\overline{x_1}-a_2\overline{x_2}) \overline{x_1}=0
\\
\Leftrightarrow
\frac{\sum_i{x_{1,i}y_i}}{n}-a_1\frac{\sum_i{x_{1,i}^2}}{n}-a_2\frac{\sum_i{x_{1,i}x_{2,i}}}{n}-(\overline{y}-a_1\overline{x_1}-a_2\overline{x_2} ) \overline{x_1}=0
\\
\Leftrightarrow
S_{xy}-a_1S^{2}_{x_{1}}-a_2S_{x_1x_2}=0
\\
\Leftrightarrow
\hat{a}_1=\frac{S_{x_1y}-a_2S_{x_1x_2}}{S^{2}_{x_{1}}}
$$
De la même manière, on obtient :
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.
$$
\hat{a}_2=\frac{S_{x_2y}-a_1S_{x_1x_2}}{S^{2}_{x_{2}}}
$$
## Méthode de calculs du modèle procustéen
L'analyse procustéenne peut être assimilée à un modèle linéaire entre deux tableaux:
$$
\hat{a}_1=\frac{S_{x_1y}-\frac{S_{x_2y}-a_1S_{x_1x_2}}{S^{2}_{x_{2}}}S_{x_1x_2}}{S^{2}_{x_{1}}}
$$
- Le tableau réponse : $Y$
- Le tableau explicatif : $X$
Si l'ensemble des variables de chaque tableau appartiennent à $\mathbb{R}$:
$$
\hat{a}_1-\hat{a}_1\frac{S^2_{x_1x_2}}
{S^2_{x_2}S^2_{x_1}}=\frac{S_{x_1y}}
{S^2_{x_1}}
-\frac{S_{x_1x_2}S_{x_2y}}
{S^2_{x_2}S^2_{x_1}}
$$
- 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:
$$
\hat{a}_1(1-\frac{S^2_{x_1x_2}}
{S^2_{x_2}S^2_{x_1}})=\frac{S_{x_1y}}
{S^2_{x_1}}
-\frac{S_{x_1x_2}S_{x_2y}}
{S^2_{x_2}S^2_{x_1}}
$$
$$ Rot(X|Y) = XVU' $$
$$
\hat{a}_1\frac{S^2_{x_2}S^2_{x_1}-S^2_{x_1x_2}}
{S^2_{x_2}S^2_{x_1}}=
\frac{S_{x_1y}}
{S^2_{x_1}}
-\frac{S_{x_1x_2}S_{x_2y}}
{S^2_{x_2}S^2_{x_1}}
$$
$$
\hat{a}_1\frac{S^2_{x_2}S^2_{x_1}-S^2_{x_1x_2}}
{S^2_{x_2}S^2_{x_1}}=
\frac{S_{x_1y}S^2_{x_2}}
{S^2_{x_1}S^2_{x_2}}
-\frac{S_{x_1x_2}S_{x_2y}}
{S^2_{x_2}S^2_{x_1}}
$$
$$
\hat{a}_1
=
\frac{S_{x_1y}S^2_{x_2}
-S_{x_1x_2}S_{x_2y}}
{S^2_{x_2}S^2_{x_1}-S^2_{x_1x_2}}
$$
- Le facteur d'homothétie $a$ se calcule selon :
D'un point de vue matricielle, les coefficients des pentes peuvent s'écrire à partir des matrices de covariances des variables explicatives (X) et des covariances entre la matrice explicative (X) et et de la variable à expliquer Y:
$$ a = \frac{\sum diag(\Lambda)}{\sum diag(X'X)}$$
Soit $\tilde{\theta}(y)$ l'estimation des différentes paramètres de la régression multiple :
dans notre cas :
$$
\tilde{\theta}(y) = (\hat a_1, \hat a_2)
$$
$\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$:
En généralisant :
$$
\tilde{\theta}(y) = (X'X)^{-1}X'Y
a=\frac{Cov_{XY}}{Var_X}
$$
## Calcul matriciel : calculs de parts de variances expliquées d'un tableau Y à partir d'un tableau X explicatif
Soit Y la matrice à expliquer et X la matrice explicative.
Par procuste on a trois opérations :
* translation (revient à aligner les barycentres : centrer )
* échelle : homothétie
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$.
* rotation trouver l'angle qui minimise : méthode des moindres carrés
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$:
### Calcul de l'inertie d'une matrice et de co-inertie entre deux matrices
Seule contrainte de la méthode, la rotation doit être effectuée en dernier.
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'$.
$$
CoI_{XY} = \frac{\sum{diag(\Lambda)}}{n-1}
$$
On cherche : $$ Min|Y-(1c'+\rho X A)| $$
On se retrouve à minimiser comme en relation linéaire l'écart entre X notre tableau de données explicative et la matrice Y (qui représente l'approximation en régression linéaire) :
$$1c'+\rho X A $$
Comme pour la Covariance, $CoI_{XY}=CoI{YX}$ et $CoI_{XX}=I_X$
### Calcul du coefficient de corrélation entre deux matrices
avec $c'$ translation (centrage), $\rho$ échelle (normalisation), $A$ rotation
Par analogie on définit $R_{XY}$ comme le coefficient de corrélation entre deux matrices $X$ et $Y$ de la manière suivante :
D'un point de vue matriciel : il faut que A'A=I orthogonal, diag(A'A)=1
$$
R_{XY} = \frac{CoI_{XY}}{\sqrt{I_{X}I_{Y}}}
$$
Les matrices centrées de X et Y sont notées :
$$\tilde{Y}, \tilde{X}$$
### Calcul des facteurs d'échelle en procuste multiple à $k$ tableaux
Décomposition en valeur singulière : avec U'U=V'V=I
$$\tilde{Y}' \tilde{X} = U \Lambda V'$$
La matrice de rotation est donc estimée par : $$\hat{A}=VU'$$
On note à partir de maintenant :
Paramètre d'échelle qui correspond à c dans procruste :
$$\hat{\rho}=\frac{trace(\hat{A}\tilde{Y}'\tilde{X})}{trace(\tilde{Y}'\tilde{X})}$$
- $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.
Si on fait l'analogie avec la fonction programmée dans procrustre (équivalence : translation <=> b, échelle <=> c et rotation A) :
$a$ se calcule par :
$$
Xrot = \hat{\rho} X A \\
b = \overline{Y}-\hat{\rho}\overline{X} A
$$
a = (CoI_{XX})^{-1} CoI_{YX}
$$
Si la translation n'a pas encore été faite, on a :
Les prédictions du modèle peuvent donc s'écrirent:
$$
Xrot = 1c'+ \hat{\rho} X A
\\
b = \overline{Y}-\hat{\rho}\overline{X} A
\widehat{Rot(Y|M)} = a_1 Rot(X_1|M) + a_2 Rot(X_2|M)+\,...\,+ a_k Rot(X_k|M)
$$
Xrot est la nouvelle estimation de X qui a été transformé pour "être représenté dans la même base" que Y
### Interaction entre deux tableaux explicatifs
Et du coup la différence entre Xrot et Y représente la part non expliquée du tableau Y par X
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$.
$$
\sum(Y-X_{rot})^2 = SRC_{residuel}
a_i = (b_i Rot(X_j|M) + c_i) \\
a_j = (b_j Rot(X_i|M) + c_j)
$$
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 :
## Calcul matriciel : calculs de parts de variances expliquées d'un tableau Y à partir de deux tableaux X1 et X2 explicatifs
Soit Y la matrice à expliquer, X1 et X2 deux matrices explicatives.
Par procuste on a trois opérations :
* translation (revient à aligner les barycentres : centrer )
* échelle : homothétie
* rotation trouver l'angle qui minimise : méthode des moindres carrés
Seule contrainte de la méthode, la rotation doit être effectuée en dernier.
On cherche : $$ Min|Y-[(1c_1'+\rho_1 X_1 A_1)+(1c_2'+\rho_2 X_2 A_2)]| $$
On se retrouve à minimiser comme en relation linéaire l'écart entre X1 et X2 nos tableaux de données explicatives et la matrice Y (qui représente l'approximation en régression linéaire) :
$$(1c_1'+\rho_1 X_1 A_1)+(1c_2'+\rho_2 X_2 A_2) $$
avec $c_1'$ et $c_2'$ translation (centrage), $\rho1$ et $\rho_2$ échelle (normalisation), $A_1$ et $A_2$ rotation
D'un point de vue matriciel, les matrices de rotations doivent vérifier les conditions suivantes : A'A=I orthogonal, diag(A'A)=1
Les matrices centrées sont notées :
$$\tilde{Y}, \tilde{X_1}, \tilde{X_2}$$
Décomposition en valeur singulière avec :
$$
U_1'U_1 = V_1'V_1 = I
\\
U_2'U_2 = V_2'V_2 = I
$$
$$
\tilde{Y}' \tilde{X_1} = U_1 \Lambda V_1'
\begin{aligned}
\widehat{Rot(Y|M)} = & a_1 \cdot Rot(X_1|M) + a_2 \cdot Rot(X_2|M) \\
\\
\tilde{Y}' \tilde{X_2} = U_2 \Lambda V_2'
$$
La matrice de rotation est donc estimée par :
$$
\hat{A_1} = V_1 U_1'
= & (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) \\
\\
\hat{A_2} = V_2 U_2'
= & 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}
$$
Paramètre d'échelle qui correspond à c dans procruste :
En renommant les facteurs d'échelle les prédictions s'écrivent:
$$
\hat{\rho_1} = \frac{trace(\hat{A_1}\tilde{Y}'\tilde{X_1})}{trace(\tilde{Y}'\tilde{X_1})}
\\
\hat{\rho_2} = \frac{trace(\hat{A_2}\tilde{Y}'\tilde{X_2})}{trace(\tilde{Y}'\tilde{X_2})}
\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)
$$
Si on fait l'analogie avec la fonction programmée dans procrustre (équivalence : translation <=> b, échelle <=> c et rotation A) :
Ce qui revient à construire un nouveau modèle sans interaction incluant un tableau explicatif supplémentaire
$$
X_1rot = \hat{\rho_1} X_1 A_1
\\
X_2rot = \hat{\rho_2} X_2 A_2
\\
b = \overline{Y}-\hat{\rho_1}\overline{X_1} A_1 -\hat{\rho_2}\overline{X_2} A_2
$$
Si la translation n'a pas encore été faite, on a :
$$
X_1rot = 1c_1'+ \hat{\rho_1} X_1 A_1
\\
X_2rot = 1c_2'+ \hat{\rho_2} X_2 A_2
\\
b = \overline{Y}-\hat{\rho_1}\overline{X_1} A_1 -\hat{\rho_2}\overline{X_2} A_2
X_{1,2}=Rot(X_1|M) \cdot Rot(X_2|M)
$$
Ce principe peut être généralisé à l'interaction entre plus de deux tableaux.
X1rot et X2 rot sont les nouvelles estimations de X1 et X2 qui ont été transformé pour "être représenté dans la même base" que Y.
## Partition de l'inertie de $Y$ selon le modèle procustéen
avec $c'$ translation (centrage), $\rho$ échelle (normalisation), $A$ rotation
D'un point de vue matriciel : il faut que A'A=I orthogonal, diag(A'A)=1
Les matrices centrées de X et Y sont notées :
$$\tilde{Y}, \tilde{X}$$
Décomposition en valeur singulière : avec U'U=V'V=I
$$\tilde{X}' \tilde{Y}=U \Lambda V'$$
La matrice de rotation est donc estimée par : $$\hat{A}=VU'$$
Paramètre d'échelle qui correspond à c dans procruste :
$$\hat{\rho}=\frac{trace(\hat{A}\tilde{X}'\tilde{Y})}{trace(\tilde{Y}'\tilde{Y})}$$
Si on fait l'analogie avec la fonction programmée dans procrustre (équivalence : translation <=> b, échelle <=> c et rotation A) :
La variation total $SCT$ est définie comme suit:
$$
Yrot=\hat{\rho} Y A
\\
b=\overline{X}-\hat{\rho}\overline{Y} A
SCT = I_Y * (n-1)
$$
? j'aurai plus tôt écrit, la translation a déjà du être faite dans procruste avant cette étape??
La variation résiduelle non expliquée par le modèle est :
$$
Yrot=1c'+ \hat{\rho} Y A
\\
b=\overline{X}-\hat{\rho}\overline{Y} A
SCR = \sum{(Rot(Y|M)-\widehat{Rot(Y|M)})^2}
$$
? Yrot est la nouvelle estimation de Y qui a été transformé pour "être représenté dans la même base" que X. ???
Et du coup la différence entre X et Yrot représente la part non expliquée par Y du tableau de données X
Selon l'approche de Scherrer on propose de définir $SCE_i$ la contribution de $X_i$ à la variation de $Y$ peut s'écrire :
$$
\sum(X-Y_{rot})^2=SRC_{residuel}
SCE_i = a_i \; R_{X_iY} \; SCT
$$
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.
## ProcMod usage
### The procmd.frame data structure
## How to build a procrustean model
```{r}
A = matrix(1:6,nrow=3)
B = matrix(1:9,nrow=3)
C = matrix(1:12,nrow=3)
pf = procmod.frame(A,B,D=C,
row.names = c('x','y','z'))
pf
```
### Build a model
#### Prepare the data
Three datasets are used
```{r}
euk = read.csv("australia.euk.reads.plot.csv",
head=TRUE,
row.names = 1)
bac = read.csv("australia.bac.reads.plot.csv",
head=TRUE,
row.names = 1)
env = read.csv("biogegraphy_and_environment.csv",
head=TRUE,
row.names = 1)
```
Only numerical environmental variables are kept.
```{r}
env=env[,lapply(env,class)=="numeric"]
geovar = which(colnames(env) %in% c('Latitude','Longitude'))
soilvar= which(colnames(env) %in% c("KLg", "pH", "AlLg",
"FeLg", "PLg", "SLg",
"CaLg", "MgLg", "MnLg",
"CNratio", "CLg", "NLg"))
climvar= which(colnames(env) %in% c("Aspect", "TempSeasonality",
"MaxMonTemp", "Elevation",
"MeanMonTempRange", "AnnMeanTemp",
"Isothemality"))
geo = env[,geovar]
library(ProcMod)
```
Environmental data are centered and reduced
```{r}
env = scale(env,scale = TRUE)
soil = env[,soilvar]
climat = env[,climvar]
```
### Loading of the demo data
Eukaryote and bacterial data are arranged in the same order than environmental data.
They consist in two MOTUs tables describing bacterial and eukaryiote communities at 21 sites accros the eastern coast of Australia.
```{r}
euk=euk[rownames(env),]
bac=bac[rownames(env),]
data("eukaryotes")
data("bacteria")
```
Relative frequency tables for Eukaryota and Bacteria are root square transforms which corresponds to Hellinger transformation.