Commit a972db98 by Eric Coissac

### Patch bug for regression with a single explicative variable

parent e733ab31
 ... ... @@ -6,7 +6,6 @@ require(matlib) #' @author Christelle Gonindard-Melodelima #' @export mvar = function(...) { ctrace <- function(MAT) sum(MAT^2) Xs <- list(...) Xnames=names(Xs) ... ... @@ -16,7 +15,6 @@ mvar = function(...) { stop("Matrices have different number of rows: ", cat(nXs)) } Xmeans<- sapply(Xs,colMeans) Xs <- lapply(Xs,scale,scale = FALSE) nX = length(Xs) ... ... @@ -24,11 +22,15 @@ mvar = function(...) { Xx <- rep(1:nX,nX) Xy <- rep(1:nX,rep(nX,nX)) XXs <- mapply(function(x,y) crossprod(Xs[[x]], Xs[[y]]),Xx,Xy) XXs <- mapply(function(x,y) crossprod(Xs[[x]], Xs[[y]]), Xx,Xy, SIMPLIFY = FALSE) 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 ... ... @@ -71,6 +73,7 @@ pm = function (model,data) { vars = lapply(vars, as.matrix) nrvars= sapply(vars, nrow) if (any(nrvars!=nrvars[1])) { ... ... @@ -81,21 +84,63 @@ pm = function (model,data) { nvar = length(vars) - 1 data.cov = do.call(mcov,vars) data.cov = do.call(mvar,vars) xy.cov = data.cov[2:(nvar+1),1,drop=FALSE] xx.cov = data.cov[2:(nvar+1),2:(nvar+1)] Y = vars[[1]] Xs= vars[-1] Ymean <- colMeans(Y) Xmeans<- lapply(Xs,colMeans) YC <- scale(Y,scale=FALSE) XCs <- lapply(Xs,scale,scale = FALSE) XYs <- lapply(XCs,function(x) crossprod(YC, x)) sol_yxs <- lapply(XYs,svd) A_xys <- lapply(sol_yxs, function(x) x$v %*% t(x$u)) Xrots <- mapply(function(x,a) x %*% a, XCs,A_xys, SIMPLIFY = FALSE) if (length(xx.cov)>1) ixx.cov=inv(xx.cov) else ixx.cov=1/xx.cov scales=ixx.cov %*% xy.cov scales=inv(xx.cov) %*% xy.cov dim(scales)=NULL b = Ymean - Reduce(function(a,b) a+b, mapply(function(d,coef,rot) coef * (d %*% rot), Xmeans, scales, A_xys, SIMPLIFY = FALSE)) yhat = mapply(function(coef,xrot) coef * xrot, scales, Xrots, SIMPLIFY = FALSE) init=matrix(rep(0,nrow(yhat[[1]])*ncol(yhat[[1]])), nrow=nrow(yhat[[1]])) yhat = Reduce(function(a,b) a+b, yhat, init) # yhat = sweep(yhat,2,b,"+") names(scales)=attr(terms,"term.labels") results = list(coefficients = scales, fitted.values=yhat, residuals=Y - yhat, call=cl, model=vars, terms = terms # model=vars, model=c(list(Y=YC),XCs), terms = terms, A = A_xys ) class(results)="pm" ... ... @@ -113,6 +158,10 @@ print.pm = function(m,...) { invisible(m) } anova.pm = function(m,...) { } #' Estimate a procuste regression of several table on a reference one. #' #' @author Eric Coissac ... ... @@ -149,11 +198,15 @@ mprocuste = function (Y, ...) Xx <- rep(1:nX,nX) Xy <- rep(1:nX,rep(nX,nX)) XXs <- mapply(function(x,y) crossprod(Xs[x], Xs[y]),Xx,Xy) XXs <- mapply(function(x,y) crossprod(Xs[x], Xs[y]), Xx,Xy, SIMPLIFY = FALSE) 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) Xrots <- mapply(function(x,a) x %*% a, Xs,A_xys, SIMPLIFY = FALSE) CovYXs = lapply(sol_yxs, function(sol) sum(sol$d)) ... ...
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!