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