Commit cea08e2e by Eric Coissac

### Premiere version de myproc3

parent 74341c64
 ... ... @@ -69,6 +69,14 @@ text(bac.pco.li[,1:2],labels = rownames(bac.pco.li),cex=0.6)  {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) text(env.pca.li[,1:2],labels = rownames(env.pca.li),cex=0.6)  {r} euk.sct = sum(euk.pco.li**2) bac.sct = sum(bac.pco.li**2)  ... ... @@ -143,54 +151,127 @@ myproc = function (X, Y, scale = TRUE, symmetric = FALSE, scores = "sites",  {r} myproc3 = function (X, Y, Z) myproc3 = function (Y, X1, X2) { if (nrow(X) != nrow(Y)) stop("Matrices have different number of rows: ", nrow(X), " and ", nrow(Y)) if (nrow(X) != nrow(Z)) stop("Matrices have different number of rows: ", nrow(X), " and ", nrow(Y)) 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) X <- scale(X, scale = FALSE) Y <- scale(Y, scale = FALSE) Z <- scale(Z, scale = FALSE) Ymean <- colMeans(Y) X1mean <- colMeans(X1) X2mean <- colMeans(X2) xmean <- apply(X, 2, mean) ymean <- apply(Y, 2, mean) Y <- scale(Y, scale = FALSE) X1 <- scale(X1, scale = FALSE) X2 <- scale(X2, scale = FALSE) XY <- crossprod(X, Y) sol1 <- svd(XY) A1 <- sol1$v %*% t(sol1$u) YX1 <- crossprod(Y, X1) sol_yx1 <- svd(YX1) A_yx1 <- sol_yx1$v %*% t(sol_yx1$u) XZ <- crossprod(X, Z) sol2 <- svd(XZ) A2 <- sol2$v %*% t(sol2$u) YX2 <- crossprod(Y, X2) sol_yx2 <- svd(YX2) A_yx2 <- sol_yx2$v %*% t(sol_yx2$u) Yrot <- Y %*% A1 Zrot <- Z %*% A2 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) VarY = ctrace(Y) VarX1 = ctrace(X1) VarX2 = ctrace(X2) SdX1 = sqrt(VarX1) SdX2 = sqrt(VarX2) a1 = (CovYX1 * VarX2 - CovYX2 * CovX1X2)/(VarX1*VarX2-CovX1X2^2) a2 = (CovYX2 - a1 * CovX1X2) / VarX2 b = Ymean - a1 * X1mean %*% A_yx1 - a2 * X2mean %*% A_yx2 CovXY = crossprod(X, Yrot) CovXY = crossprod(X, Zrot) CovYZ = crossprod(Yrot, Zrot) SCX1 = sum((X1rot * a1)^2) SCX2 = sum((X2rot * a2)^2) 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 Yhat = X1rot * a1 + X2rot * a2 SCR = sum((Y - Yhat)^2) SCT = VarY SCI = SCT - SCX1 - SCX2 - SCR ddl.t = (nrow(Y)-1)^2 ddl.r = ddl.t - 3 vt = SCT/ddl.t vx1= SCX1 vx2= SCX2 vi = SCI vr = SCR/ddl.r fx1=vx1/vr fx2=vx2/vr fi =vi/vr pf.x1=1-pf(fx1,1,ddl.r) pf.x2=1-pf(fx2,1,ddl.r) pf.i =1-pf(fi ,1,ddl.r) ddl = (nrow(Y)-1)^2-1 sd.a1 = sqrt((VarY/VarX1 - a1^2)/ddl) sd.a2 = sqrt((VarY/VarX2 - a2^2)/ddl) t.a1 = a1/sd.a1 t.a2 = a2/sd.a2 p.a1 = 1 - pt(t.a1,ddl) p.a2 = 1 - pt(t.a2,ddl) res = list(coefficients = list(a1=a1,a2=a2,b=b), sd = list(a1=sd.a1,a2=sd.a2), t = list(a1=t.a1,a2=t.a2), p.value = list(a1=p.a1,a2=p.a2), anova = list(X1=SCX1/SCT, X2=SCX2/SCT, X1X2=SCI/SCT, Res=SCR/SCT), SunSq = list(X1=SCX1, X2=SCX2, X1X2=SCI, Res=SCR), MeanSq = list(X1=vx1,X2=vx2,X1X2=vi,Res=vr), Fvalue = list(X1=fx1,X2=fx2,X1X2=fi), Fvalue = list(X1=pf.x1,X2=pf.x2,X1X2=pf.i) ) class(res) <- "procustes3" return(res) }  {r} xxx = myproc3(euk.pco.li,bac.pco.li,env.pca.li) xxx  {r} xxx = myproc3(bac.pco.li,euk.pco.li,env.pca.li) xxx  ## Estimation des coefficients lors d'une rgression linaire multiple (2 variables) ... ... @@ -250,6 +331,54 @@ $$\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. ... ...
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