Commit e1365710 by Eric Coissac

Version of the module related to the manuscript on IRLs

parent 142cdcb8
......@@ -13,17 +13,22 @@ RoxygenNote: 6.1.1
Imports: MASS,
permute,
expm,
Matrix,
mvtnorm,
stats,
roxygen2
doParallel,
foreach
Suggests: knitr,
rmarkdown,
roxygen2,
vegan
VignetteBuilder: knitr
Collate:
'IR.R'
'internals.R'
'procmod_frame.R'
'multivariate.R'
'covls.R'
'corls_test.R'
'procuste.R'
'simulate.R'
# Generated by roxygen2: do not edit by hand
S3method("$",procmod.corls)
S3method("$",procmod.varls)
S3method("$<-",procmod.frame)
S3method("[",procmod.frame)
S3method("[[<-",procmod.frame)
......@@ -7,9 +9,11 @@ S3method(as.data.frame,dist)
S3method(as.list,procmod.frame)
S3method(as.procmod.frame,array)
S3method(as.procmod.frame,list)
S3method(as.procmod.frame,matrix)
S3method(as.procmod.frame,pm)
S3method(as.procmod.frame,procmod.frame)
S3method(dim,procmod.frame)
S3method(names,procmod.varls)
S3method(ortho,data.frame)
S3method(ortho,dist)
S3method(ortho,matrix)
......@@ -21,6 +25,8 @@ export(as.procmod.frame)
export(bicenter)
export(corls)
export(corls.partial)
export(corls.test)
export(icor)
export(is.euclid)
export(is.procmod.frame)
export(nmds)
......@@ -33,4 +39,5 @@ export(simulate_correlation)
export(simulate_matrix)
export(varls)
import(MASS)
import(expm)
import(doParallel)
import(foreach)
#' @author Christelle Gonindard-Melodelima
#' @author Eric Coissac
NULL
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
icor <- function(x, y = NULL) {
if (is.data.frame(y)) {
y <- as.matrix(y)
}
if (is.data.frame(x)) {
x <- as.matrix(x)
}
if (!is.matrix(x) && is.null(y)) {
stop("supply both 'x' and 'y' or a matrix-like 'x'")
}
if (!(is.numeric(x) || is.logical(x))) {
stop("'x' must be numeric")
}
stopifnot(is.atomic(x))
if (!is.null(y)) {
if (!(is.numeric(y) || is.logical(y))) {
stop("'y' must be numeric")
}
stopifnot(is.atomic(y))
}
if (!is.matrix(x)) {
x <- t(t(x))
}
if (is.null(y)) {
y <- x
}
if (!is.matrix(y)) {
y <- t(t(y))
}
xc <- scale(x, scale = FALSE, center = TRUE)
yc <- scale(y, scale = FALSE, center = TRUE)
n <- nrow(x)
cov <- crossprod(xc, yc) / (n - 1)
print(cov)
sdx <- apply(x, MARGIN = 2, sd)
sdy <- apply(y, MARGIN = 2, sd)
rcov <- sqrt(1 / (n - 1)) * (sdx %o% sdy)
print(rcov)
s <- sign(cov)
icov <- (s * cov - rcov) * s
print(icov)
isdx <- sqrt(1 - sqrt(1 / (n - 1))) * sdx
isdy <- sqrt(1 - sqrt(1 / (n - 1))) * sdy
ipearson <- icov / (isdx %o% isdy)
ipearson
}
#' @include covls.R
#'
#' @author Christelle Gonindard-Melodelima
#' @author Eric Coissac
NULL
#' Generate permutation matrix according to a schema.
#'
#' @param perm
#' @param n
#' @param strata
#'
#'
#' The permutation schema is defined using the `how` function.
#' The implementation of this function is inspired
#' from the VEGAN package and reproduced here to avoid an extra
#' dependency on an hidden vegan function.
#'
getPermuteMatrix = function(perm, n, strata = NULL)
{
if (length(perm) == 1) {
perm <- permute::how(nperm = perm)
}
if (!missing(strata) && !is.null(strata)) {
if (inherits(perm, "how") && is.null(permute::getBlocks(perm)))
permute::setBlocks(perm) <- strata
}
if (inherits(perm, "how"))
perm <- permute::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
}
#' Compute the person correlation matrix of K coordinate matrices
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
corls.test <- function(...,
permutations = permute::how(nperm = 999),
p.adjust.method="holm") {
eps <- sqrt(sqrt(.Machine$double.eps))
xs <- list(...)
if (length(xs) == 1) {
x <- xs[[1]]
if (is.procmod.frame(x)) {
xs <- x
} else {
xs <- procmod.frame(x)
}
}
else {
xs <- as.procmod.frame(xs)
}
x_names <- names(xs)
xs <- ortho(xs)
cov <- varls(xs, nperm = 0)
lcov <- cov - eps
ngreater <- array(0,dim = dim(cov))
n <- nrow(xs)
nx <- length(xs)
pmatrix <- getPermuteMatrix(permutations, n)
if (ncol(pmatrix) != n) {
stop(gettextf(
"'permutations' have %d columns, but data have %d observations",
ncol(pmatrix), n
))
}
npermutation <- nrow(pmatrix)
for (i in seq_len(npermutation)) {
ps <- sample(1:npermutation,
size = nx,
replace = FALSE
)
rcov = varls(as.procmod.frame(
lapply(
1:nx,
function(j) xs[[j]][pmatrix[ps[j], ], ]
)
),
nperm = 0
)
ngreater <- ngreater + (
rcov >= lcov)
}
p_values <- ngreater / npermutation
diag(p_values) <- 0
c_p_values <- p.adjust(p_values[upper.tri(p_values,diag = FALSE)],
method = p.adjust.method,
n = (nx - 1) * nx / 2)
p_values[upper.tri(p_values,diag = FALSE)] <- c_p_values
p_values <- as.matrix(Matrix::forceSymmetric(p_values, uplo = "U"))
colnames(p_values) <- x_names
rownames(p_values) <- x_names
p_values
}
#' @include procmod_frame.R
#' @include multivariate.R
#'
#' @import expm
#' @import doParallel
#' @import foreach
#'
#' @author Christelle Gonindard-Melodelima
#' @author Eric Coissac
NULL
library(doParallel)
registerDoParallel(1)
#' Compute the trace of a square matrix.
#'
#' @param X a square matrix
......@@ -15,7 +18,6 @@ NULL
#' @examples
#' m <- matrix(1:16, nrow = 4)
#' ProcMod:::.Trace(m)
#'
#' @note Internal function do not use.
#'
#' @rdname internal.getPermuteMatrix
......@@ -24,6 +26,17 @@ NULL
#'
.Trace <- function(X) sum(diag(X))
.estimate_mode <- function(x) {
d <- density(x)
d$x[which.max(d$y)]
}
.var2cor <- function(c) {
v <- sqrt(diag(c))
vv <- v %o% v
c / vv
}
#' Compute the variance, covariance matrix of K coordinate matrices.
#'
#' Covariance between two matrices is defined as the sum of the
......@@ -31,8 +44,13 @@ NULL
#' the same number of rows.
#'
#' @param ... the set of matrices
#' @param permutation
#' @param rcovls
#' @param nrand number of randomisation used to estimate the mean
#' covariance observed between two random matrix.
#' @param p.adjust.method the multiple test correction method used
#' to adjust p values. \code{p.adjust.method} belongs
#' one of the folowing values: "holm", "hochberg", "hommel",
#' "bonferroni", "BH", "BY", "fdr", "none". The default is
#' set to "holm".
#'
#' @examples
#' # Build Three matrices of 3 rows.
......@@ -40,109 +58,150 @@ NULL
#' 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)
#' varls(A, B, C)
#' varls(A = A, B = B, C = C)
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
varls <- function(..., nperm = 100, rcovls = FALSE) {
if (is.numeric(nperm) &&
length(nperm) == 1 &&
nperm == 0) {
nperm <- NULL
varls <- function(...,
nrand = 100,
p.adjust.method = "holm") {
if (!is.null(nrand) && length(nrand) > 1) {
stop("nrand must be a single numeric value or the NULL value")
}
if (is.null(nperm)) {
rcovls <- FALSE
}
if (nrand == 0) nrand <- NULL
xs <- list(...)
Xs <- list(...)
if (length(Xs) == 1) {
x <- Xs[[1]]
if (length(xs) == 1) {
x <- xs[[1]]
if (is.procmod.frame(x)) {
Xs <- x
} else if (is.pm(x)) {
return(x$cov)
xs <- x
} else {
Xs <- procmod.frame(x)
xs <- procmod.frame(x)
}
}
else {
Xs <- as.procmod.frame(Xs)
xs <- as.procmod.frame(xs)
}
x_names <- names(Xs)
x_names <- names(xs)
Xs <- ortho(Xs)
xs <- ortho(xs)
XXs <- as.procmod.frame(mapply(tcrossprod,
Xs,
SIMPLIFY = FALSE
))
# xxs <- as.procmod.frame(mapply(tcrossprod,
# xs,
# SIMPLIFY = FALSE
# ))
nx <- length(Xs)
N <- nrow(Xs) - 1
nx <- length(xs)
n <- nrow(xs)
n1 <- n - 1
cov_xxs <- matrix(0, nrow = nx, ncol = nx)
cov_xxs <- Matrix::Matrix(0, nrow = nx, ncol = nx)
# Computes the standard Covls covariance matrix
for (i in seq_len(nx))
for (j in i:nx) {
cov_xxs[i, j] <- Re(.Trace(sqrtm(XXs[[i]] %*% XXs[[j]])))
}
for (j in i:nx)
#cov_xxs[i, j] <- Re(.Trace(expm::sqrtm(xxs[[i]] %*% xxs[[j]])))
cov_xxs[i, j] <- sum(svd(crossprod(xs[[i]],xs[[j]]))$d) / n1
if (!is.null(nperm)) {
r_cov_xxs <- array(0, dim = c(nperm, nx, nx))
v_xs <- vector(mode = "list", nx)
l_xs <- vapply(Xs, FUN = function(x) prod(dim(x)), 0)
d_xs <- vapply(Xs, FUN = function(x) dim(ortho(x)), c(0, 0))
cov_xxs <- as.matrix(Matrix::forceSymmetric(cov_xxs, uplo = "U"))
for (i in seq_len(nx))
v_xs[i] <- list(t(sqrt(colMeans(Xs[i]^2) - colMeans(Xs[i])^2)
%*% t(rep(1, nrow(Xs[i])))))
# Computes Covls under null hypothesis
if (!is.null(nrand)) {
n_r_greater <- array(0, dim = dim(cov_xxs))
for (k in seq_len(nperm)) {
# Generate the random matrix for the kth simulation
v_xs <- vector(mode = "list", nx)
for (i in seq_len(nx))
v_xs[[i]] <- var(xs[[i]])
# s_cov_xxs <- array(0, dim = c(nx, nx, nrand))
# for(k in seq_len(nrand)) {
# #s1_cov_xxs <- matrix(0, nrow = nx, ncol = nx)
# r_xs <- vector(mode = "list", nx)
# r_ys <- vector(mode = "list", nx)
# for (i in seq_len(nx)) {
# r_xs[[i]] <- MASS::mvrnorm(n,
# rep(0,nrow(v_xs[[i]])),
# Sigma = v_xs[[i]],
# empirical = TRUE)
# r_ys[[i]] <- MASS::mvrnorm(n,
# rep(0,nrow(v_xs[[i]])),
# Sigma = v_xs[[i]],
# empirical = TRUE)
# }
#
# for (i in seq_len(nx))
# for (j in i:nx)
# s_cov_xxs[i, j, k] <- sum(svd(crossprod(r_xs[[i]],r_ys[[j]]))$d)
# }
s_cov_xxs <- foreach(k = seq_len(nrand),
.combine = cbind) %dopar% {
s1_cov_xxs <- matrix(0, nrow = nx, ncol = nx)
r_xs <- vector(mode = "list", nx)
for (i in seq_len(nx))
r_xs[i] <- list(tcrossprod(scale(array(rnorm(l_xs[i]),
dim = d_xs[, i]
)) * v_xs[[i]]))
r_ys <- vector(mode = "list", nx)
for (i in seq_len(nx)) {
r_xs[[i]] <- MASS::mvrnorm(n,
rep(0,nrow(v_xs[[i]])),
Sigma = v_xs[[i]],
empirical = TRUE)
r_ys[[i]] <- MASS::mvrnorm(n,
rep(0,nrow(v_xs[[i]])),
Sigma = v_xs[[i]],
empirical = TRUE)
}
for (i in seq_len(nx))
for (j in i:nx)
r_cov_xxs[k, i, j] <- .Trace(sqrtm(r_xs[[i]] %*% r_xs[[j]]))
}
}
s1_cov_xxs[i, j] <- sum(svd(crossprod(r_xs[[i]],r_ys[[j]]))$d)
r_cov_xxs <- apply(r_cov_xxs,
MARGIN = c(2, 3),
FUN = function(i) MASS::fitdistr(
Re(i),
"normal"
)$estimate["mean"]
)
s1_cov_xxs / n1
for (i in seq_len(nx))
for (j in i:nx) {
cov_xxs[i, j] <- max(cov_xxs[i, j] - r_cov_xxs[i, j], 0)
cov_xxs[j, i] <- cov_xxs[i, j]
r_cov_xxs[j, i] <- r_cov_xxs[i, j]
}
}
cov_xxs <- cov_xxs / N
colnames(cov_xxs) <- x_names
rownames(cov_xxs) <- x_names
dim(s_cov_xxs) = c(nx, nx, nrand)
s_cov_xxs <- s_cov_xxs
r_cov_xxs <- apply(s_cov_xxs,
MARGIN = c(1,2),
FUN = mean
)
for (k in seq_len(nrand))
n_r_greater <- n_r_greater + (s_cov_xxs[, , k] >= cov_xxs)
r_cov_xxs <- as.matrix(Matrix::forceSymmetric(r_cov_xxs, uplo = "U"))
cov_xxs <- pmax(cov_xxs - r_cov_xxs, 0)
colnames(r_cov_xxs) <- x_names
rownames(r_cov_xxs) <- x_names
cov_xxs <- Re(cov_xxs)
p_values <- n_r_greater / nrand
if (rcovls) {
attr(cov_xxs, "rcovls") <- r_cov_xxs / N
c_p_values <- p.adjust(p_values[upper.tri(p_values, diag = FALSE)],
method = p.adjust.method
)
p_values[upper.tri(p_values, diag = FALSE)] <- c_p_values
p_values <- as.matrix(Matrix::forceSymmetric(p_values, uplo = "U"))
colnames(p_values) <- x_names
rownames(p_values) <- x_names
#attr(cov_xxs, "scovls") <- s_cov_xxs
attr(cov_xxs, "rcovls") <- as.matrix(r_cov_xxs)
attr(cov_xxs, "p.value") <- p_values
attr(cov_xxs, "nrand") <- nrand
}
attr(cov_xxs, "nperm") <- nperm
colnames(cov_xxs) <- x_names
rownames(cov_xxs) <- x_names
return(make_subS3Class(Re(cov_xxs), "procmod.varls"))
make_subS3Class(cov_xxs, "procmod.varls")
}
......@@ -151,19 +210,23 @@ varls <- function(..., nperm = 100, rcovls = FALSE) {
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
corls <- function(..., nperm = 100, rcovls = FALSE) {
cov <- varls(..., nperm = nperm, rcovls = rcovls)
corls <- function(..., nrand = 100,
p.adjust.method = "holm") {
cov <- varls(...,
nrand = nrand,
p.adjust.method = p.adjust.method
)
s <- sqrt(diag(cov))
vv <- outer(s, s)
vv <- s %o% s
rls <- cov / vv
rls <- make_subS3Class(cov / vv, "procmod.corls")
attr(rls, "nperm") <- attr(cov, "nperm")
if (rcovls) {
attr(rls, "rcovls") <- attr(rls, "rcovls") / vv
if (!is.null(attr(cov, "rcovls"))) {
attr(rls, "nrand") <- attr(cov, "nrand")
attr(rls, "rcorls") <- attr(cov, "rcovls") / vv
}
return(rls)
make_subS3Class(rls, "procmod.corls")
}
#' Compute the person partial correlation matrix of K coordinate matrices
......@@ -171,19 +234,17 @@ corls <- function(..., nperm = 100, rcovls = FALSE) {
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
corls.partial <- function(..., nperm = 100, rcovls = FALSE) {
rls <- corls(..., nperm = nperm, rcovls = rcovls)
corls.partial <- function(..., nrand = 100) {
rls <- corls(..., nrand = nrand)
C <- solve(rls)
D <- sqrt(diag(C) %o% diag(C))
rp <- make_subS3Class(C / D, "procmod.corls")
S <- sqrt(diag(C))
D <- S %o% S
rp <- C / D
attr(rp, "nperm") <- attr(rls, "nperm")
attr(rp, "nrand") <- attr(rls, "nrand")
attr(rp, "rcovls") <- attr(rls, "rcovls")
if (rcovls) {
attr(rp, "rcovls") <- attr(rls, "rcovls")
}
return(rp)
make_subS3Class(rp, "procmod.corls")
}
#' Compute the person partial correlation matrix of K coordinate matrices
......@@ -193,11 +254,28 @@ corls.partial <- function(..., nperm = 100, rcovls = FALSE) {
#' @export
print.procmod.varls <- function(x, ...) {
class(x) <- "matrix"
attr(x, "permutations") <- NULL
attr(x, "nrand") <- NULL
attr(x, "rcovls") <- NULL
attr(x, "p.value") <- NULL
print(x)
}
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
`$.procmod.varls` <- function(x, name) {
attr(x,name)
}
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
names.procmod.varls <- function(x) {
names(attributes(x))
}
#' Compute the person partial correlation matrix of K coordinate matrices
#'
#' @author Eric Coissac
......@@ -205,7 +283,17 @@ print.procmod.varls <- function(x, ...) {
#' @export
print.procmod.corls <- function(x, ...) {
class(x) <- "matrix"
attr(x, "permutations") <- NULL
attr(x, "nrand") <- NULL
attr(x, "rcovls") <- NULL
attr(x, "rcorls") <- NULL
attr(x, "p.value") <- NULL
print(x)
}
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
`$.procmod.corls` <- function(x, name) {
attr(x,name)
}
......@@ -310,6 +310,24 @@ as.procmod.frame.pm <- function(data, ...) {
vars.procmod(terms(data), data$model)
}
#' Coerce to a ProcMod Frame.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
as.procmod.frame.matrix <- function(data, ...) {
l <- vector(mode = "list", length = ncol(data))
for (i in seq_len(ncol(data))) {
l[[i]] <- data[, i]
}
if (!is.null(colnames(data))) {
names(l) <- colnames(data)
}
as.procmod.frame(l)
}
#' Dimensions of a Matrix Frame.
#'
#' @author Eric Coissac
......
......@@ -8,7 +8,7 @@ NULL
#' the variances to 1 with all the covariances set to 0.
#'
#' @export
simulate_matrix <- function(n, p, equal.var = FALSE) {
simulate_matrix <- function(n, p, equal.var = TRUE) {
new <- rnorm(n * p, mean = 0, sd = 1)
dim(new) <- c(n, p)
......@@ -25,36 +25,35 @@ simulate_matrix <- function(n, p, equal.var = FALSE) {
#' Simulate n points of dimension p correlated with a reference matrix.
#'
#' @export
simulate_correlation <- function(reference, p, r2, equal.var = FALSE) {
simulate_correlation <- function(reference, p, r2, equal.var = TRUE) {
n <- nrow(reference)
maxdim <- max(ncol(reference), p)
noise <- simulate_matrix(n, p, equal.var = equal.var)
if (maxdim == p && maxdim > ncol(reference)) {