Commit e733ab31 by Eric Coissac

### Send a first version of the main module source file

parent 949e2285
 require(matlib) #' Compute the variance, covariance matrix of K coordinate matrix #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export mvar = function(...) { ctrace <- function(MAT) sum(MAT^2) Xs <- list(...) Xnames=names(Xs) nXs= sapply(Xs, nrow) if (any(nXs!=nXs[1])) { stop("Matrices have different number of rows: ", cat(nXs)) } Xmeans<- sapply(Xs,colMeans) Xs <- lapply(Xs,scale,scale = FALSE) nX = length(Xs) Xx <- rep(1:nX,nX) Xy <- rep(1:nX,rep(nX,nX)) XXs <- mapply(function(x,y) crossprod(Xs[[x]], Xs[[y]]),Xx,Xy) sol_xxs <- lapply(XXs,svd) CovXXs = sapply(sol_xxs, function(sol) sum(sol$d)) dim(CovXXs)=c(nX,nX) colnames(CovXXs)=Xnames rownames(CovXXs)=Xnames return(CovXXs) } #' Compute the correlation matrix of K coordinate matrix #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export mcor = function(...) { cov = mvar(...) s = sqrt(diag(cov)) vv= outer(s,s) return(cov/vv) } #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export pm = function (model,data) { cl=match.call() terms = terms(model) if (missing(data)) data=NULL vars = eval(attr(terms(model),"variables"), envir = data, enclos = globalenv()) names(vars)=rownames(attr(terms,"factors")) cvars= sapply(vars, class) if (! all(cvars %in% c("matrix","data.frame"))) { stop("All the variables are not matrices", cat(cvars)) } vars = lapply(vars, as.matrix) nrvars= sapply(vars, nrow) if (any(nrvars!=nrvars[1])) { stop("Matrices have different number of rows: ", cat(nrvars)) } attr(vars,"terms")=terms nvar = length(vars) - 1 data.cov = do.call(mcov,vars) xy.cov = data.cov[2:(nvar+1),1,drop=FALSE] xx.cov = data.cov[2:(nvar+1),2:(nvar+1)] scales=inv(xx.cov) %*% xy.cov dim(scales)=NULL names(scales)=attr(terms,"term.labels") results = list(coefficients = scales, call=cl, model=vars, terms = terms ) class(results)="pm" return(results) } #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export print.pm = function(m,...) { cat(sprintf("\nCall:\n%s\n\nCoefficients:\n", deparse(m$call))) print(m$coefficients) invisible(m) } #' Estimate a procuste regression of several table on a reference one. #' #' @author Eric Coissac #' @author Christelle Gonindard-Melodelima #' @export mprocuste = function (Y, ...) { require(matlib) ctrace <- function(MAT) sum(MAT^2) Xs <- get_list_from_ellipsis(...) Xs <- list(...) nY = nrow(Y) nXs= sapply(Xs, nY) nXs= sapply(Xs, nrow) if (any(nXs!=nY)) { stop("Matrices have different number of rows: ", nY, ... ... @@ -26,15 +144,19 @@ mprocuste = function (Y, ...) 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) nX = length(Xs) Xx <- rep(1:nX,nX) Xy <- rep(1:nX,rep(nX,nX)) XXs <- mapply(function(x,y) crossprod(Xs[x], Xs[y]),Xx,Xy) sol_xxs <- lapply(XXs,svd) A_xxs <- lapply(sol_xxs, function(x) x$v %*% t(x$u)) Xrots <- mapply(function(x,a) x %*% a,Xs,A_xys) X1rot <- X1 %*% A_yx1 X2rot <- X2 %*% A_yx2 CovYXs = lapply(sol_yxs, function(sol) sum(sol$d)) CovYX1 = sum(sol_yx1$d) CovYX2 = sum(sol_yx2$d) CovX1X2 = sum(sol_x1x2\$d) VarY = ctrace(Y) ... ...
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!