Commit 31363933 by Eric Coissac

Merge branch 'master' of git.metabarcoding.org:lecasofts/ProcMod

# Conflicts:
#	Bricolage.Rmd
#	DESCRIPTION
#	NAMESPACE
#	R/mcov.R
#	R/permute.R
#	R/simulate.R
#	man/corls.Rd
#	man/corls.partial.Rd
#	man/internal.getPermuteMatrix.Rd
#	man/simulate_correlation.Rd
#	man/simulate_matrix.Rd
#	man/varls.Rd
parents 8a115d95 ccb4d408
......@@ -84,6 +84,7 @@ rp^2
```{r}
library(ProcMod)
<<<<<<< HEAD
r2s = seq(from=0.0, to = 0.9, by = 0.1)
n.sim=1000
r2.obs = matrix(0,nrow = n.sim,ncol = length(r2s))
......@@ -98,6 +99,20 @@ for (j in seq_along(r2s)) {
# v = varls(d1,d2,permutations = 10)
# s = sqrt(diag(v))
# r2.obs[i,j] = (v / outer(s,s))[1,2]^2
=======
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
>>>>>>> ccb4d408c4936c553b8d097d94022285a7bc89a7
}
}
......@@ -106,6 +121,7 @@ cbind(R2_theo=r2s,
R2_obs=colMeans(r2.obs),
R2_sd=apply(r2.obs,MARGIN = 2,FUN = sd))
```
<<<<<<< HEAD
```{r}
a= simulate_matrix(40,1)
......@@ -409,4 +425,6 @@ abline(h=1)
legend("topright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5)
```
=======
>>>>>>> ccb4d408c4936c553b8d097d94022285a7bc89a7
......@@ -30,7 +30,7 @@ Collate:
'mprocuste.R'
'anova.pm.R'
'covls.R'
'permute.R'
'mcov.R'
'plot.pm.R'
'procuste.R'
'simnull.R'
......
......@@ -47,9 +47,8 @@ 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)
#' @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
#' 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
#' mvar(A,B,C)
#' mvar(A=A,B=B,C=C)
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
varls = function(...,permutations = how(nperm = 999)) {
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
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)) {
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
corls = function(...,permutations = how(nperm = 999)) {
cov = varls(...,permutations = permutations)
s = sqrt(diag(cov))
vv= outer(s,s)
return(cov/vv)
}
#' Compute the person partial correlation matrix of K coordinate matrices
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
corls.partial = function(...,permutations = how(nperm = 999)) {
C = solve(corls(...,permutations = permutations))
D = sqrt(diag(C) %o% diag(C))
return(C/D)
}
#' @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
}
#' @import mvtnorm
#'
NULL
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/covls.R
% Please edit documentation in R/covls.R, R/mcov.R
\name{corls}
\alias{corls}
\title{Compute the person correlation matrix of K coordinate matrices}
\usage{
corls(..., nperm = 100, rcovls = FALSE)
corls(..., permutations = how(nperm = 999))
corls(..., permutations = how(nperm = 999))
}
\description{
Compute the person correlation matrix of K coordinate matrices
Compute the person correlation matrix of K coordinate matrices
}
\author{
Eric Coissac
Christelle Gonindard-Melodelima
Eric Coissac
Christelle Gonindard-Melodelima
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/covls.R
% Please edit documentation in R/covls.R, R/mcov.R
\name{corls.partial}
\alias{corls.partial}
\title{Compute the person partial correlation matrix of K coordinate matrices}
\usage{
corls.partial(..., nperm = 100, rcovls = FALSE)
corls.partial(..., permutations = how(nperm = 999))
corls.partial(..., permutations = how(nperm = 999))
}
\description{
Compute the person partial correlation matrix of K coordinate matrices
Compute the person partial correlation matrix of K coordinate matrices
}
\author{
Eric Coissac
Christelle Gonindard-Melodelima
Eric Coissac
Christelle Gonindard-Melodelima
}
\name{hello}
\alias{hello}
\title{Hello, World!}
\usage{
hello()
}
\description{
Prints 'Hello, world!'.
}
\examples{
hello()
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/covls.R, R/permute.R
% Please edit documentation in R/covls.R, R/mcov.R
\name{.Trace}
\alias{.Trace}
\alias{.getPermuteMatrix}
\title{Compute the trace of a square matrix.}
\usage{
.Trace(X)
.getPermuteMatrix(perm, N)
.Trace(X)
}
\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.}
\item{X}{a square matrix}
}
\value{
the trace of X
a matrix of integer of size nperm x n. Each line corresponds to a
permutation of n element
the trace of X
}
\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.
Compute the trace of a square matrix.
Compute the trace of a square matrix.
}
\note{
Internal function do not use.
......@@ -43,7 +32,8 @@ Internal function do not use.
m = matrix(1:16,nrow=4)
ProcMod:::.Trace(m)
ProcMod:::.getPermuteMatrix(perm=how(nperm=5),N=10)
m = matrix(1:16,nrow=4)
ProcMod:::.Trace(m)
}
\author{
......@@ -51,7 +41,7 @@ Eric Coissac
Christelle Gonindard-Melodelima
Christelle Gonindard-Melodelima
Eric Coissac
Christelle Gonindard-Melodelima
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/covls.R
% Please edit documentation in R/covls.R, R/mcov.R
\name{varls}
\alias{varls}
\title{Compute the variance, covariance matrix of K coordinate matrices.}
\usage{
varls(..., nperm = 100, rcovls = FALSE)
varls(..., permutations = how(nperm = 999))
varls(..., permutations = how(nperm = 999))
}
\arguments{
\item{...}{the set of matrices}
\item{rcovls}{}
\item{...}{the set of matrices}
}
\description{
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.
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.
}
\examples{
# Build Three matrices of 3 rows.
......@@ -25,9 +33,21 @@ the same number of rows.
varls2(A,B,C)
varls2(A=A,B=B,C=C)
# 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
mvar(A,B,C)
mvar(A=A,B=B,C=C)
}
\author{
Eric Coissac
Christelle Gonindard-Melodelima
Eric Coissac
Christelle Gonindard-Melodelima
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/mcov.R
\name{varls2}
\alias{varls2}
\title{Compute the variance, covariance matrix of K coordinate matrices.}
\usage{
varls2(..., permutations = how(nperm = 999))
}
\arguments{
\item{...}{the set of matrices}
}
\description{
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.
}
\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
Christelle Gonindard-Melodelima
}
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