R/concorscano.R

In BMconcor: CONCOR for Structural- And Regular-Equivalence Blockmodeling

#' simultaneous concorgmcano
#'
#' concorgmcano with the set of r solutions simultaneously optimized
#'
#' This function uses the concors function
#'
#' @param x are the \code{n} times \code{p} and \code{n} times \code{q} matrices of \code{p} and \code{q} centered column
#' @param y See \code{x}
#' @param px A row vector which contains the numbers pi, i=1,...,kx, of the kx subsets xi of x : sum(pi)=sum(px)=p. px is the partition vector of x
#' @param py The partition vector of y. A row vector containing the numbers \code{qi} for \code{i = 1,...,ky} of the \code{ky} subsets \code{yi} of \code{y : sum(qi)=sum(py)=q}.
#' @param r The number of wanted successive solutions rmax <= min(min(px),min(py),n)
#'
#' @return A \code{list} with following components:
#' \item{cx}{a \code{n} times \code{r} matrix of the r canonical components of x}
#' \item{cy}{a \code{n.ky} times \code{r} matrix. The ky blocks cyi of the rows n*(i-1)+1 : n*i contain the r canonical components relative to Yi}
#' \item{cov2}{a \code{ky} times \code{r} matrix; each column \code{k} contains \code{ky} squared covariances \eqn{\mbox{cov}(x*u[,k],y_i*v_i[,k])^2}, the partial measures of link}
#'
#' @author \enc{Lafosse, R.}{R. Lafosse}
#'
#' @references Hanafi & Lafosse (2001) Generalisation de la regression lineaire simple pour analyser la dependance de K ensembles de variables avec un K+1 eme.  Revue de Statistique Appliquee vol.49, n.1
#'
#' @examples
#'
#' x <- matrix(runif(50),10,5);y <- matrix(runif(90),10,9)
#' x <- scale(x);y <- scale(y)
#' cca <- concorscano(x,c(2,3),y,c(3,2,4),2)
#' cca$rho2[1,1,]
#'
#' @export

concorscano <-
  function(x,px,y,py,r) {
    if (sum(px) != dim(x)[2] | sum(py) != dim(y)[2] ) stop("px or py IS NOT SUITABLE")
    n <- dim(x)[1]
    kx <- length(px)
    rx <- matrix(0,1,kx)
    Px <- NULL
    cux=c(0,cumsum(px))
    for (j in 1:kx) {
      s <- svd(x[,(cux[j]+1):cux[j+1]])
      rx[j]<-sum(s$d > max(c(n,px[j]))*s$d[1]*1e-8)
      Px <- cbind(Px,s$u[,1:rx[j]]*sqrt(n))
    }
    cux <- c(0,cumsum(rx))
    Px <- matrix(Px,nrow=n)
    ky <- length(py)
    ry <- matrix(0,1,ky)
    Py <- NULL
    cuy=c(0,cumsum(py))
    for (j in 1:ky) {
      s <- svd(y[,(cuy[j]+1):cuy[j+1]])
      ry[j]<-sum(s$d > max(c(n,py[j]))*s$d[1]*1e-8)
      Py <- cbind(Py,s$u[,1:ry[j]]*sqrt(n))
    }
    if (r > min(c(min(ry),min(rx),n))) stop("r IS TOO HIGH")
    cuy <- c(0,cumsum(ry))
    Py <- matrix(Py,nrow=n)
    s <- concors(Px,rx,Py,ry,r)
    cy <- matrix(0,n*ky,r)
    cx <- matrix(0,n*kx,r)

    for  (j in 1:kx) {
      cx[((j-1)*n+1):(j*n),]<-matrix(Px[,(cux[j]+1):cux[j+1]],nrow=n)%*%s$u[(cux[j]+1):cux[j+1],]
    }
    for  (j in 1:ky) cy[((j-1)*n+1):(j*n),]<-matrix(Py[,(cuy[j]+1):cuy[j+1]],nrow=n)%*%s$v[(cuy[j]+1):cuy[j+1],]

    list(cx=cx,cy=cy,rho2=s$cov2)
  }