diff --git a/debroussaillage.html b/debroussaillage.html deleted file mode 100644 index 09eff2e..0000000 --- a/debroussaillage.html +++ /dev/null @@ -1,600 +0,0 @@ - - - - - - - - - - - - - - -Untitled - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - -
library(ade4)
-
## Warning: package 'ade4' was built under R version 3.4.3
-
library(vegan)
-
## Warning: package 'vegan' was built under R version 3.4.3
-
## Loading required package: permute
-
## Loading required package: lattice
-
## This is vegan 2.4-6
-
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)
-
env=env[,lapply(env,class)=="numeric"]
-euk=euk[rownames(env),]
-bac=bac[rownames(env),]
-
env = sweep(env,MARGIN = 2,colMeans(env),'-')
-env = sweep(env,MARGIN = 2,apply(env,2,sd),'/')
-
euk=sqrt(euk)
-#euk = sweep(euk,MARGIN = 2,colMeans(euk),'-')
-#euk = sweep(euk,MARGIN = 2,apply(euk,2,sd),'/')
-
bac=sqrt(bac)
-#bac = sweep(bac,MARGIN = 2,colMeans(bac),'-')
-#bac = sweep(bac,MARGIN = 2,apply(bac,2,sd),'/')
-
euk.dist=vegdist(euk,method = "jaccard")
-euk.pco =dudi.pco(euk.dist,full = TRUE)
-
-euk.pco.li = euk.pco$li
-dim(euk.pco.li)
-
## [1] 62 61
-
plot(euk.pco.li[,1:2],cex=0)
-text(euk.pco.li[,1:2],labels = rownames(euk.pco.li),cex=0.6)
-

-
bac.dist=vegdist(bac,method = "jaccard")
-bac.pco =dudi.pco(bac.dist,full = TRUE)
-
-bac.pco.li = bac.pco$li
-dim(bac.pco.li)
-
## [1] 62 61
-
plot(bac.pco.li[,1:2],cex=0)
-text(bac.pco.li[,1:2],labels = rownames(bac.pco.li),cex=0.6)
-

-
euk.sct = sum(euk.pco.li**2)
-bac.sct = sum(bac.pco.li**2)
-
x  = procrustes(euk.pco.li,bac.pco.li,scale = TRUE)
-
y  = procrustes(bac.pco.li,euk.pco.li,scale = TRUE)
-
SCR.E.B = sum((euk.pco.li-x$Yrot)^2)
-SCT.E   = sum((euk.pco.li- colMeans(euk.pco.li))^2)
-1-SCR.E.B/SCT.E
-
## [1] 0.8496279
-
SCR.B.E = sum((bac.pco.li-y$Yrot)^2)
-SCT.B   = sum((bac.pco.li- colMeans(bac.pco.li))^2)
-1-SCR.B.E/SCT.B
-
## [1] 0.8496279
-
myproc = function (X, Y, scale = TRUE, symmetric = FALSE, scores = "sites", 
-    ...) 
-{
-    X <- scores(X, display = scores, ...)
-    Y <- scores(Y, display = scores, ...)
-    if (nrow(X) != nrow(Y)) 
-        stop("Matrices have different number of rows: ", nrow(X), 
-            " and ", nrow(Y))
-    if (ncol(X) < ncol(Y)) {
-        warning("X has fewer axes than Y: X adjusted to comform Y\n")
-        addcols <- ncol(Y) - ncol(X)
-        for (i in 1:addcols) X <- cbind(X, 0)
-    }
-    ctrace <- function(MAT) sum(MAT^2)
-    c <- 1
-    if (symmetric) {
-        X <- scale(X, scale = FALSE)
-        Y <- scale(Y, scale = FALSE)
-        X <- X/sqrt(ctrace(X))
-        Y <- Y/sqrt(ctrace(Y))
-    }
-    xmean <- apply(X, 2, mean)
-    ymean <- apply(Y, 2, mean)
-    if (!symmetric) {
-        X <- scale(X, scale = FALSE)
-        Y <- scale(Y, scale = FALSE)
-    }
-    XY <- crossprod(X, Y)
-    sol <- svd(XY)
-    A <- sol$v %*% t(sol$u)
-    
-    print(A)
-    
-    if (scale) {
-        c <- sum(sol$d)/ctrace(Y)
-    }
-    Yrot <- c * Y %*% A
-    b <- xmean - c * ymean %*% A
-    R2 <- ctrace(X) + c * c * ctrace(Y) - 2 * c * sum(sol$d)
-    reslt <- list(Yrot = Yrot, X = X, ss = R2, rotation = A, 
-        translation = b, scale = c, xmean = xmean, symmetric = symmetric, 
-        call = match.call())
-    reslt$svd <- sol
-    class(reslt) <- "procrustes"
-    reslt
-}
-
myproc3 = function (Y, X1, X2) 
-{
-    if (nrow(Y) != nrow(X1)) 
-        stop("Matrices have different number of rows: ", nrow(Y), 
-            " and ", nrow(X1))
-    if (nrow(Y) != nrow(X2)) 
-        stop("Matrices have different number of rows: ", nrow(Y), 
-            " and ", nrow(X1))
-
-    ctrace <- function(MAT) sum(MAT^2)
-    
-    c <- c(1,1)
-        
-    Y  <- scale(Y,  scale = FALSE)
-    X1 <- scale(X1, scale = FALSE)
-    X2 <- scale(X2, scale = FALSE)
-
-    Ymean <- colMeans(Y)
-    X1mean <- colMeans(X1)
-    X2mean <- colMeans(X2)
-
-    YX1 <- crossprod(Y, X1)
-    sol_yx1 <- svd(YX1)
-    A_yx1 <- sol_yx1$v %*% t(sol_yx1$u)
-    
-    YX2 <- crossprod(Y, X2)
-    sol_yx2 <- svd(YX2)
-    A_yx2 <- sol_yx2$v %*% t(sol_yx2$u)
-
-    X1X2 <- crossprod(X1, X2)
-    sol_x1x2 <- svd(X1X2)
-    A_x1x2 <- sol_x1x2$v %*% t(sol_x1x2$u)
-
-    X1rot <- X1 %*% A_yx1
-    X2rot <- X2 %*% A_yx2
-    
-    CovYX1  = sum(sol_yx1$d)
-    CovYX2  = sum(sol_yx2$d)
-    CovX1X2 = sum(sol_x1x2$d)
-    
-    
-}
-
-

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 :

- -

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} -\]

-
- - - - -
- - - - - - - -