We test with the new unbiased RLS

Bricolage.Rmd

+---
+title: "R Notebook"
+output:
+  html_document:
+    df_print: paged
+---
+
+```{r}
+Xs = test2
+nperm=100
+```
+
+
+```{r}
+ix = c(4,3,6,5,8,7)
+iy = 1
+
+
+r = corls(Xs,permutations = how(nperm = nperm))
+rx = r[ix,ix]
+ry = r[iy,ix]
+rxi = r[ix[(ix %% 2)>0],ix[(ix %% 2)>0]]
+ryi = r[iy,ix[(ix %% 2)>0]]
+```
+
+```{r}
+r2 = (solve(rx) %*% ry) * ry
+r2
+```
+
+```{r}
+sum(r2)
+```
+
+```{r}
+r2 = (solve(rxi) %*% ryi) * ryi
+r2
+```
+
+```{r}
+sum(r2)
+```
+
+```{r}
+v = varls(Xs,permutations = how(nperm = nperm))
+vt = v[c(iy,ix[(ix %% 2)>0]),c(iy,ix[(ix %% 2)>0])] * (nrow(Xs)-1)
+vr = v[c(iy+1,ix[(ix %% 2)==0]),c(iy,ix[(ix %% 2)>0])] * (nrow(Xs)-1)
+
+vt
+vr
+```
+
+```{r}
+#diag(vr)=0
+vp = vt - vr
+vp
+```
+
+```{r}
+d = sqrt(diag(vp))
+n = d %o% d
+
+r = vp / n
+r
+```
+
+```{r}
+rx = r[ix[(ix %% 2)==0]/2,ix[(ix %% 2)==0]/2]
+ry = r[iy,ix[(ix %% 2)==0]/2]
+r2 = (solve(rx) %*% ry) * ry
+r2
+```
+```{r}
+sum(r2)
+```
+
+```{r}
+rp = solve(r)
+rp = rp / sqrt(diag(rp) %o% diag(rp))
+rp
+rp^2
+```
+
+
+```{r}
+library(ProcMod)
+r2s = seq(from=0.1, to = 0.9, by = 0.1)
+n.sim=20
+r2.obs = matrix(0,nrow = n.sim,ncol = length(r2s))
+d1 = simulate_matrix(10,2000)
+for (j in seq_along(r2s)) {
+  print(c("r2",as.character(r2s[j])))
+  for (i in seq_len(n.sim)) {
+    print(c("round",as.character(i)))
+    d2 = simulate_correlation(10,3000,d1,r2s[j])
+    v = varls2(d1,d2,permutations = 10)
+    s = sqrt(diag(v))
+    r2.obs[i,j] = (v / outer(s,s))[1,2]^2
+  }
+}
+
+boxplot(r2.obs)
+cbind(R2_theo=r2s,
+      R2_obs=colMeans(r2.obs),
+      R2_sd=apply(r2.obs,MARGIN = 2,FUN = sd))
+```
+

DESCRIPTION

@@ -11,6 +11,9 @@ Encoding: UTF-8
 LazyData: true
 RoxygenNote: 6.1.1
 Imports: MASS,
+    permute,
+    expm,
+    mvtnorm,
     stats,
     roxygen2
 Suggests: knitr,
@@ -27,6 +30,8 @@ Collate:
     'mprocuste.R'
     'anova.pm.R'
     'mcov.R'
+    'permute.R'
     'plot.pm.R'
     'procuste.R'
+    'simulate.R'
     'zzzz.R'

NAMESPACE

@@ -28,13 +28,12 @@ S3method(residuals,pm)
 S3method(subset,procmod.frame)
 export(as.procmod.frame)
 export(bicenter)
+export(corls)
+export(corls.partial)
 export(is.euclid)
 export(is.pm)
 export(is.procmod.frame)
-export(mcor)
-export(mcor.partial)
 export(model.procmod.default)
-export(mvar)
 export(nmds)
 export(ortho)
 export(pca)
@@ -42,6 +41,13 @@ export(pcoa)
 export(pm)
 export(procmod.frame)
 export(protate)
+export(simulate_correlation)
+export(simulate_matrix)
+export(varls)
+export(varls2)
 export(vars.procmod)
 export(weighted.residuals)
 import(MASS)
+import(expm)
+import(mvtnorm)
+import(permute)

R/mcov.R

 #' @include procmod.frame.R
 #' @include multivariate.R
 #'
+#' @import expm
+#'
+#' @author Christelle Gonindard-Melodelima
+#' @author Eric Coissac
 NULL
 
+#' Compute the trace of a square matrix.
+#'
+#'
+#' @param X a square matrix
+#' @return the trace of X
+#'
+#' @examples
+#'     m = matrix(1:16,nrow=4)
+#'     ProcMod:::.Trace(m)
+#'
+#' @note Internal function do not use.
+#'
+#' @rdname internal.getPermuteMatrix
+#' @author Eric Coissac
+#' @author Christelle Gonindard-Melodelima
+#'
+.Trace = function(X) sum(diag(X))
+
 #' Compute the variance, covariance matrix of K coordinate matrices.
 #'
 #' Covariance between two matrices is defined as the sum of the
@@ -20,9 +42,10 @@ NULL
 #'    mvar(A,B,C)
 #'    mvar(A=A,B=B,C=C)
 #'
-#' @author Eric Coissac & Christelle Gonindard-Melodelima
+#' @author Eric Coissac
+#' @author Christelle Gonindard-Melodelima
 #' @export
-mvar = function(...) {
+varls = function(...,permutations = how(nperm = 999)) {
 
   Xs <- list(...)
   if (length(Xs)==1) {
@@ -37,34 +60,170 @@ mvar = function(...) {
   else
       Xs=as.procmod.frame(Xs)
 
-    Xnames=names(Xs)
+  Xnames=names(Xs)
 
-    Xs <- ortho(Xs)
+  Xs <- ortho(Xs)
 
-    nX = length(Xs)
+  XXs = as.procmod.frame(mapply(tcrossprod, Xs,
+                        SIMPLIFY = FALSE))
 
-    Xx <- rep(1:nX,nX)
-    Xy <- rep(1:nX,rep(nX,nX))
+  nX = length(Xs)
+  N  = nrow(Xs)-1
 
-    CovXXs <- mapply(function(x,y) sum(svd(crossprod(Xs[[x]], Xs[[y]]))$d),
-                  Xx,Xy,
-                  SIMPLIFY = TRUE) / (nrow(Xs)-1)
+  if (! is.null(permutations)) {
+    pmatrix = .getPermuteMatrix(perm=permutations,N=nrow(Xs))
+    rCovXXs = matrix(0, nrow = nX * 2, ncol = nX * 2)
+    pval = matrix(0, nrow = nX * 2, ncol = nX * 2)
+    for (i in 1:(2*nX))
+      for (j in 1:(2*nX)) {
+        if (i %% 2 && j %% 2) {
+          rCovXXs[i,j]=.Trace(sqrtm(XXs[[ceiling(i / 2)]] %*% XXs[[ceiling(j/2)]])) / N
+          pval[i,j]=-1
+        }
+        else if (i ==j) {
+          vv = numeric(nrow(pmatrix))
+          for (k in seq_len(nrow(pmatrix))) {
+            d = Xs[[ceiling(i / 2)]][pmatrix[k,],]
+            dd = tcrossprod(d)
+            vv[k] = Re(.Trace(sqrtm(dd %*% dd)))
+          }
+          rCovXXs[i,j]=mean(vv) / N
+          pval[i,j]=shapiro.test(vv)$p.value
+        }
+        else if (i <=j) {
+          vv = numeric(nrow(pmatrix))
+          for (k in seq_len(nrow(pmatrix))) {
+            d = Xs[[ceiling(i / 2)]][pmatrix[k,],]
+            dd = tcrossprod(d)
+            vv[k] = Re(.Trace(sqrtm(dd %*% XXs[[ceiling(j/2)]])))
+          }
+          rCovXXs[i,j]=mean(vv) / N
+          rCovXXs[j,i]=rCovXXs[i,j]
+          pval[i,j]=shapiro.test(vv)$p.value
+          pval[j,i]=shapiro.test(vv)$p.value
+        }
+      }
 
+    CovXXs=rCovXXs
+    Xnames = rep(Xnames,rep(2,length(Xnames)))
+    Xsuff  = rep("",length(Xnames))
+    Xsuff[seq(from=2,to=length(Xnames),by=2)]="r_"
+    Xnames = mapply(paste0, Xsuff,Xnames,collapse="")
+    colnames(CovXXs)=Xnames
+    rownames(CovXXs)=Xnames
+    colnames(pval)=Xnames
+    rownames(pval)=Xnames
+
+    print(pval)
+  }
+  else {
+    Xx <- rep(1:nX,nX)
+    Xy <- rep(1:nX,rep(nX,nX))
+    CovXXs <- mapply(function(x,y) .Trace(sqrtm(XXs[[x]] %*% XXs[[y]])),
+                     Xx,Xy,
+                     SIMPLIFY = TRUE) / N
 
     dim(CovXXs)=c(nX,nX)
     colnames(CovXXs)=Xnames
     rownames(CovXXs)=Xnames
+  }
+
+  return(Re(CovXXs))
+}
+
+
+#' Compute the variance, covariance matrix of K coordinate matrices.
+#'
+#' Covariance between two matrices is defined as the sum of the
+#' sigular values of the X'Y matrix. All the matrices must have
+#' the same number of rows.
+#'
+#' @param ... the set of matrices
+#'
+#' @examples
+#'    # Build Three matrices of 3 rows.
+#'    A <- matrix(1:9,nrow=3)
+#'    B <- matrix(10:15,nrow=3)
+#'    C <- matrix(20:31,nrow=3)
+#'    # compute the variance covariance matrix
+#'    varls2(A,B,C)
+#'    varls2(A=A,B=B,C=C)
+#'
+#' @author Eric Coissac
+#' @author Christelle Gonindard-Melodelima
+#' @export
+varls2 = function(...,permutations = how(nperm = 999)) {
 
-    return(CovXXs)
+  Xs <- list(...)
+  if (length(Xs)==1) {
+    x = Xs[[1]]
+    if (is.procmod.frame(x))
+      Xs=x
+    else if (is.pm(x))
+      return(x$cov)
+    else
+      Xs=procmod.frame(x)
+  }
+  else
+      Xs=as.procmod.frame(Xs)
+
+  Xnames=names(Xs)
+
+  Xs <- ortho(Xs)
+
+  XXs = as.procmod.frame(mapply(tcrossprod, Xs,
+                        SIMPLIFY = FALSE))
+
+  nX = length(Xs)
+  N  = nrow(Xs)-1
+
+  CovXXs = matrix(0, nrow = nX, ncol = nX)
+
+  for (i in seq_len(nX))
+    for (j in i:nX) {
+        CovXXs[i,j] = .Trace(sqrtm(XXs[[i]] %*% XXs[[j]]))
+      }
+
+  if (! is.null(permutations)) {
+    pmatrix = .getPermuteMatrix(perm=permutations,N=nrow(Xs))
+    nP = nrow(pmatrix)
+    rCovXXs = array(0,dim=c(nX,nX,nP))
+
+    for (k in seq_len(nP)) {
+      Xp = Xs[pmatrix[k,],]
+      for (i in seq_len(nX)) {
+        dd = tcrossprod(Xp[[i]])
+        for (j in i:nX)
+          rCovXXs[i,j,k] = Re(.Trace(sqrtm(dd %*% XXs[[j]])))
+      }
+    }
+
+    for (i in seq_len(nX))
+      for (j in i:nX)
+        CovXXs[i,j] = CovXXs[i,j] - mean(rCovXXs[i,j,])
+
+
+  }
+
+  for (i in seq_len(nX))
+    for (j in i:nX)
+      CovXXs[j,i] = CovXXs[i,j]
+
+  CovXXs = CovXXs / N
+  colnames(CovXXs)=Xnames
+  rownames(CovXXs)=Xnames
+
+  return(Re(CovXXs))
 }
 
+
 #' Compute the person correlation matrix of K coordinate matrices
 #'
 #' @author Eric Coissac
 #' @author Christelle Gonindard-Melodelima
 #' @export
-mcor = function(...) {
-  cov = mvar(...)
+corls = function(...,permutations = how(nperm = 999)) {
+  cov = varls(...,permutations = permutations)
   s = sqrt(diag(cov))
   vv= outer(s,s)
   return(cov/vv)
@@ -75,8 +234,8 @@ mcor = function(...) {
 #' @author Eric Coissac
 #' @author Christelle Gonindard-Melodelima
 #' @export
-mcor.partial = function(...) {
-  C = solve(mcor(...))
+corls.partial = function(...,permutations = how(nperm = 999)) {
+  C = solve(corls(...,permutations = permutations))
   D = sqrt(diag(C) %o% diag(C))
   return(C/D)
 }

R/permute.R

+#' @include internals.R
+#' @import permute
+#'
+#' @author Christelle Gonindard-Melodelima
+#' @author Eric Coissac
+#'
+NULL
+
+#' Generate a permutation matrix using the permute package
+#'
+#' The code is copied from the \code{\link[vegan]{getPermuteMatrix}} internal function
+#' of the \code{vegan} pakage. This is mainly a wrappoer around the functions of
+#' the \code{permute} package.
+#'
+#' The permutation design is controled through the \code{\link[permute]{how}} function
+#' of the \code{permute} package.
+#'
+#' @param perm a list of control values for the permutations
+#'             as returned by the function how, or the number
+#'             of permutations required, or a permutation
+#'             matrix where each row gives the permuted indices.
+#' @param N    number of elements to permute.
+#'
+#' @return a matrix of integer of size nperm x n. Each line corresponds to a
+#'         permutation of n element
+#'
+#' @note Internal function do not use.
+#'
+#' @examples
+#'     ProcMod:::.getPermuteMatrix(perm=how(nperm=5),N=10)
+#'
+#' @rdname internal.getPermuteMatrix
+#'
+#' @author Christelle Gonindard-Melodelima
+#' @author Eric Coissac
+#'
+.getPermuteMatrix = function (perm, N)
+{
+    if (length(perm) == 1) {
+        perm <- how(nperm = perm)
+    }
+
+    if (inherits(perm, "how"))
+        perm <- shuffleSet(N, control = perm)
+    else {
+        if (!is.integer(perm) && !all(perm == round(perm)))
+            stop("permutation matrix must be strictly integers: use round()")
+    }
+    if (is.null(attr(perm, "control")))
+        attr(perm, "control") <- structure(list(within = list(type = "supplied matrix"),
+            nperm = nrow(perm)), class = "how")
+    perm
+}

R/procuste.R

@@ -7,6 +7,23 @@
 #'
 NULL
 
+#' Rotate the \code{src} matrix to fit into the space of the \code{dest} matrix.
+#'
+#' The optimal rotation is computed according to the procruste methode.
+#' Rotation is based on singular value decomposition (SVD).
+#' No scaling is done, only the rotation.
+#'
+#' @param src a numeric matrix to be rotated
+#' @param dest a numeric matrix used as reference space
+#'
+#' @return a numeric matrix
+#'
+#' @examples
+#'    # Renerate a random matrix of size 10 x 15
+#'    m1 = simulate_matrix(10,15)
+#'    m2 = simulate_matrix(10,20)
+#'    mr = protate(m1,m2)
+#'
 #' @export
 protate = function(src,dest) {
   YX = crossprod(dest, src)

R/simulate.R

+#' @import mvtnorm
+#'
+NULL
+
+#' Simulate n points of dimension p.
+#'
+#' Points are simulated using the \code{\link[mvtnorm]{rmvnorm}} from
+#' the \code{mvtnorm}. The mean of each dimension is set to 0 and all
+#' the variances to 1 with all the covariances set to 0.
+#'
+#' @export
+simulate_matrix = function(n,p,method = "chol") {
+  new = rmvnorm(n,rep(0,p),method = method)
+
+  new.sd = as.numeric(sqrt(varls(new,permutations=NULL)))
+
+  new = new/new.sd
+
+  return(new)
+}
+
+#' Simulate n points of dimension p correlated with a reference matrix.
+#'
+#' @export
+simulate_correlation = function(reference,p,r2,method = "chol") {
+  n = nrow(reference)
+  maxdim = max(ncol(reference),p)
+
+  noise = simulate_matrix(n,p,method = method)
+
+  if (maxdim == p && maxdim > ncol(reference)) {
+    temp = reference
+    reference = noise
+    noise = temp
+    switched = TRUE
+  }
+  else
+    switched = FALSE
+
+  noise.rotate = protate(noise,reference)
+
+  inflate = sqrt(r2/(1-r2))
+
+  if (switched)
+    new = protate(noise.rotate * inflate + reference,
+                  reference)
+  else
+    new = protate(reference * inflate + noise.rotate,
+                  noise)
+
+  new.sd = as.numeric(sqrt(varls(new,permutations=NULL)))
+
+  new = new/new.sd
+
+  return(new)
+}
+
+#simulate_matrix_tree

man/mcor.Rd → man/corls.Rd

 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/mcov.R
-\name{mcor}
-\alias{mcor}
+\name{corls}
+\alias{corls}
 \title{Compute the person correlation matrix of K coordinate matrices}
 \usage{
-mcor(...)
+corls(..., permutations = how(nperm = 999))
 }
 \description{
 Compute the person correlation matrix of K coordinate matrices

man/mcor.partial.Rd → man/corls.partial.Rd

 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/mcov.R
-\name{mcor.partial}
-\alias{mcor.partial}
+\name{corls.partial}
+\alias{corls.partial}
 \title{Compute the person partial correlation matrix of K coordinate matrices}
 \usage{
-mcor.partial(...)
+corls.partial(..., permutations = how(nperm = 999))
 }
 \description{
 Compute the person partial correlation matrix of K coordinate matrices

man/internal.getPermuteMatrix.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mcov.R, R/permute.R
+\name{.Trace}
+\alias{.Trace}
+\alias{.getPermuteMatrix}
+\title{Compute the trace of a square matrix.}
+\usage{
+.Trace(X)
+
+.getPermuteMatrix(perm, N)
+}
+\arguments{
+\item{X}{a square matrix}
+
+\item{perm}{a list of control values for the permutations
+as returned by the function how, or the number
+of permutations required, or a permutation
+matrix where each row gives the permuted indices.}
+
+\item{N}{number of elements to permute.}
+}
+\value{
+the trace of X
+
+a matrix of integer of size nperm x n. Each line corresponds to a
+        permutation of n element
+}
+\description{
+The code is copied from the \code{\link[vegan]{getPermuteMatrix}} internal function
+of the \code{vegan} pakage. This is mainly a wrappoer around the functions of
+the \code{permute} package.
+}
+\details{
+The permutation design is controled through the \code{\link[permute]{how}} function
+of the \code{permute} package.
+}
+\note{
+Internal function do not use.
+
+Internal function do not use.
+}
+\examples{
+    m = matrix(1:16,nrow=4)
+    ProcMod:::.Trace(m)
+
+    ProcMod:::.getPermuteMatrix(perm=how(nperm=5),N=10)
+
+}
+\author{
+Eric Coissac
+
+Christelle Gonindard-Melodelima
+
+Christelle Gonindard-Melodelima
+
+Eric Coissac
+}

man/protate.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/procuste.R
+\name{protate}
+\alias{protate}
+\title{Rotate the \code{src} matrix to fit into the space of the \code{dest} matrix.}
+\usage{
+protate(src, dest)
+}
+\arguments{
+\item{src}{a numeric matrix to be rotated}
+
+\item{dest}{a numeric matrix used as reference space}
+}