Commit d5063c61 by Eric Coissac

Produce now a R package ProcMod

parent eebc9dc9
......@@ -2,10 +2,11 @@ Package: ProcMod
Type: Package
Title: What the Package Does (Title Case)
Version: 0.1.0
Author: Who wrote it
Author: Eric Coissac, Christelle Gonindard-Melodelima
Maintainer: The package maintainer <yourself@somewhere.net>
Description: More about what it does (maybe more than one line)
Use four spaces when indenting paragraphs within the Description.
License: What license is it under?
Encoding: UTF-8
LazyData: true
\ No newline at end of file
LazyData: true
RoxygenNote: 6.0.1
exportPattern("^[[:alpha:]]+")
# Generated by roxygen2: do not edit by hand
S3method(print,pm)
export(mcor)
export(mprocuste)
export(mvar)
export(pm)
......@@ -18,3 +18,4 @@ StripTrailingWhitespace: Yes
BuildType: Package
PackageUseDevtools: Yes
PackageInstallArgs: --no-multiarch --with-keep.source
PackageRoxygenize: rd,collate,namespace,vignette
mprocuste = function (Y, ...)
{
require(matlib)
ctrace <- function(MAT) sum(MAT^2)
Xs <- get_list_from_ellipsis(...)
nY = nrow(Y)
nXs= sapply(Xs, nY)
if (any(nXs!=nY)) {
stop("Matrices have different number of rows: ", nY,
" and ", cat(nXs))
}
Ymean <- colMeans(Y)
Xmeans<- sapply(Xs,colMeans)
Y <- scale(Y, scale = FALSE)
Xs <- lapply(Xs,scale,scale = FALSE)
XYs <- lapply(Xs,function(x) crossprod(Y, x))
sol_yxs <- lapply(XYs,svd)
A_xys <- lapply(sol_yxs, function(x) x$v %*% t(x$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)
VarY = ctrace(Y)
VarX1 = ctrace(X1)
VarX2 = ctrace(X2)
CovEx = matrix(c(VarX1,CovX1X2,CovX1X2,VarX2),nrow=2)
Cov2 = matrix(c(CovYX1,CovYX2),nrow=2)
pentes = inv(CovEx) %*% Cov2
print(pentes)
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
SCX1 = sum((X1rot * a1)^2)
SCX2 = sum((X2rot * a2)^2)
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) <- "mprocuste"
return(res)
}
......@@ -151,34 +151,32 @@ myproc = function (X, Y, scale = TRUE, symmetric = FALSE, scores = "sites",
```
```{r}
myproc3 = function (Y, X1, X2)
myproc3 = function (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))
require(matlib)
ctrace <- function(MAT) sum(MAT^2)
Xs <- get_list_from_ellipsis(...)
nY = nrow(Y)
nXs= sapply(Xs, nY)
c <- c(1,1)
if (any(nXs!=nY)) {
stop("Matrices have different number of rows: ", nY,
" and ", cat(nXs))
}
Ymean <- colMeans(Y)
X1mean <- colMeans(X1)
X2mean <- colMeans(X2)
Xmeans<- sapply(Xs,colMeans)
Y <- scale(Y, scale = FALSE)
X1 <- scale(X1, scale = FALSE)
X2 <- scale(X2, scale = FALSE)
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)
Xs <- lapply(Xs,scale,scale = FALSE)
XYs <- lapply(Xs,function(x) crossprod(Y, x))
sol_yxs <- lapply(XYs,svd)
A_xys <- lapply(sol_yxs, function(x) x$v %*% t(x$u))
X1X2 <- crossprod(X1, X2)
sol_x1x2 <- svd(X1X2)
......@@ -195,6 +193,13 @@ myproc3 = function (Y, X1, X2)
VarX1 = ctrace(X1)
VarX2 = ctrace(X2)
CovEx = matrix(c(VarX1,CovX1X2,CovX1X2,VarX2),nrow=2)
Cov2 = matrix(c(CovYX1,CovYX2),nrow=2)
pentes = inv(CovEx) %*% Cov2
print(pentes)
SdX1 = sqrt(VarX1)
SdX2 = sqrt(VarX2)
......@@ -272,6 +277,12 @@ xxx
xxx = myproc3(bac.pco.li,euk.pco.li,env.pca.li)
xxx
```
```{r}
xxx = myproc3(env.pca.li,bac.pco.li,euk.pco.li)
xxx
```
## Estimation des coefficients lors d'une rgression linaire multiple (2 variables)
......
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