Commit 55aee53a by Eric Coissac

A first version of the vignette describing the package usage.

parent 97d1ad1c
---
title: "ProcMod Package"
author: "Eric Coissac & Christelle Gonindard-Melodelima"
date: "`r Sys.Date()`"
output:
rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Vignette Title}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
library(ProcMod)
knitr::opts_knit$set(
collapse = TRUE,
comment = "#>",
root.dir = system.file("extdata", package="ProcMod")
)
```
## Mathematic priciples
## Estimation des coefficients lors d'une régression linéaire multiple (2 variables)
$$ \hat{y}=a_1 x_1 + a_2 x_2 + b$$
On cherche à minimiser l'erreur commise entre la prédiction et la vraie valeur (méthode des moindres carrés)
$$
\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 :
$$
\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.
$$
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 :
$$
\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 :
$$
\hat{a}_2=\frac{S_{x_2y}-a_1S_{x_1x_2}}{S^{2}_{x_{2}}}
$$
$$
\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}}}
$$
$$
\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}}
$$
$$
\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}}
$$
$$
\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}}
$$
## Calcul matriciel :
Soit X la matrice à expliquer et Y la matrice explicative.
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|X-(1c'+\rho Y A)| $$
On se retrouve à minimiser comme en relation linéaire l'écart entre X notre tableau de données explicative et la transformation de la matrice Y (qui représente l'approximation en régression linéaire) :
$$1c'+\rho Y A $$
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) :
$$Yrot=\hat{\rho} Y A
\\
b=\overline{X}-\hat{\rho}\overline{Y} A
$$
? j'aurai plus tôt écrit, la translation a déjà du être faite dans procruste avant cette étape??
$$Yrot=1c'+ \hat{\rho} Y A
\\
b=\overline{X}-\hat{\rho}\overline{Y} A
$$
? 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
$$
\sum(X-Y_{rot})^2=SRC_{residuel}
$$
## ProcMod usage
### 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"]
```
Eukaryote and bacterial data are arranged in the same order than environmental data.
```{r}
euk=euk[rownames(env),]
bac=bac[rownames(env),]
```
Environmental data are centered and reduced
```{r}
env = scale(env,scale = TRUE)
```
Relative frequency tables for Eukaryota and Bacteria are root square transforms which corresponds to Hellinger transformation.
```{r}
euk=sqrt(euk)
bac=sqrt(bac)
```
Using the `vegan` package the Jaccard distance is computed among the sites according to the Eukaryota and Bacteria markers.
```{r}
library(vegan)
euk.dist=vegdist(euk,method = "jaccard")
bac.dist=vegdist(bac,method = "jaccard")
```
Using the `ade4` package a principal coordinate analysis is done to place the sites in a cartesian space reflecting the distances among the sites.
```{r}
library(ade4)
euk.pco =dudi.pco(euk.dist,full = TRUE)
bac.pco =dudi.pco(bac.dist,full = TRUE)
```
Coordinates of the sites are extracted from the PCoA analysis
```{r}
euk.pco.li = euk.pco$li
bac.pco.li = bac.pco$li
```
```{r, fig.show='hold'}
plot(euk.pco.li[,1:2],cex=0,
main="Eukariota")
text(euk.pco.li[,1:2],
labels = rownames(euk.pco.li),
cex=0.4)
plot(bac.pco.li[,1:2],cex=0,
main="Bacteria")
text(bac.pco.li[,1:2],
labels = rownames(bac.pco.li),
cex=0.4)
```
The environmental variable are transformed using a Principal Component Analysis.
```{r}
env.pca = dudi.pca(env,scannf = FALSE,nf=nrow(env)-1)
env.pca.li = env.pca$li
dim(env.pca.li)
plot(env.pca.li[,1:2],cex=0,
main="Environmental data")
text(env.pca.li[,1:2],
labels = rownames(env.pca.li),
cex=0.4)
```
```{r}
```
## Vignette Info
Note the various macros within the `vignette` section of the metadata block above. These are required in order to instruct R how to build the vignette. Note that you should change the `title` field and the `\VignetteIndexEntry` to match the title of your vignette.
## Styles
The `html_vignette` template includes a basic CSS theme. To override this theme you can specify your own CSS in the document metadata as follows:
output:
rmarkdown::html_vignette:
css: mystyles.css
## Figures
The figure sizes have been customised so that you can easily put two images side-by-side.
```{r, fig.show='hold'}
plot(1:10)
plot(10:1)
```
You can enable figure captions by `fig_caption: yes` in YAML:
output:
rmarkdown::html_vignette:
fig_caption: yes
Then you can use the chunk option `fig.cap = "Your figure caption."` in **knitr**.
## More Examples
You can write math expressions, e.g. $Y = X\beta + \epsilon$, footnotes^[A footnote here.], and tables, e.g. using `knitr::kable()`.
```{r, echo=FALSE, results='asis'}
knitr::kable(head(mtcars, 10))
```
Also a quote using `>`:
> "He who gives up [code] safety for [code] speed deserves neither."
([via](https://twitter.com/hadleywickham/status/504368538874703872))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment