Commit ccb4d408 by Eric Coissac

### We test with the new unbiased RLS

parent 574f654b
Bricolage.Rmd 0 → 100644
 --- 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)) 
 ... ... @@ -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'
 ... ... @@ -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)
 #' @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 0 → 100644
 #' @include internals.R #' @import permute #' #' @author Christelle Gonindard-Melodelima #' @author Eric Coissac #' NULL ..getPermuteMatrix = function (perm, N, strata = NULL) { if (length(perm) == 1) { perm <- how(nperm = perm) } if (!missing(strata) && !is.null(strata)) { if (inherits(perm, "how") && is.null(getBlocks(perm))) setBlocks(perm) <- strata } 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 }
 ... ... @@ -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 0 → 100644
 #' @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
 % 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 ... ...
 % 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 ... ...
 % 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 0 → 100644
 % 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} } \value{ a numeric matrix } \description{ 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. } \examples{ # Renerate a random matrix of size 10 x 15 m1 = simulate_matrix(10,15) m2 = simulate_matrix(10,20) mr = protate(m1,m2) }
 % Generated by roxygen2: do not edit by hand % Please edit documentation in R/simulate.R \name{simulate_correlation} \alias{simulate_correlation} \title{Simulate n points of dimension p correlated with a reference matrix.} \usage{ simulate_correlation(reference, p, r2, method = "chol") }