Skip to content

P
ProcMod
Commit 142cdcb8 by Eric Coissac
Options

First round of source cleaning to focus only on correlation

parent 6f7e253c
Showing with 893 additions and 1626 deletions
.gitignore
......@@ -6,3 +6,9 @@ inst/doc
.DS_Store
vignettes.ooo
before_trash
R2_H0.csv
R2_H0.v2.csv
Simulation_R2_Null.Rmd
Simulation_R2_Null.html
demo.Rmd
demo.html
DESCRIPTION
......@@ -22,15 +22,8 @@ Suggests: knitr,
VignetteBuilder: knitr
Collate:
'internals.R'
'procmod.frame.R'
'procmod_frame.R'
'multivariate.R'
'formula.procmod.frame.R'
'omit.action.R'
'model.procmod.R'
'mprocuste.R'
'anova.pm.R'
'covls.R'
'plot.pm.R'
'procuste.R'
'simnull.R'
'simulate.R'
NAMESPACE
......@@ -3,51 +3,34 @@
S3method("$<-",procmod.frame)
S3method("[",procmod.frame)
S3method("[[<-",procmod.frame)
S3method(AIC,pm)
S3method(BIC,pm)
S3method(anova,pm)
S3method(as.data.frame,dist)
S3method(as.list,procmod.frame)
S3method(as.procmod.frame,array)
S3method(as.procmod.frame,list)
S3method(as.procmod.frame,pm)
S3method(as.procmod.frame,procmod.frame)
S3method(deviance,pm)
S3method(dim,procmod.frame)
S3method(extractAIC,pm)
S3method(formula,procmod.frame)
S3method(na.exclude,procmod.frame)
S3method(na.omit,procmod.frame)
S3method(ortho,data.frame)
S3method(ortho,dist)
S3method(ortho,matrix)
S3method(ortho,procmod.frame)
S3method(plot,pm)
S3method(print,pm)
S3method(print,procmod.corls)
S3method(print,procmod.varls)
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(model.procmod.default)
export(nmds)
export(null_varls_sim)
export(ortho)
export(pca)
export(pcoa)
export(pm)
export(procmod.frame)
export(protate)
export(simulate_correlation)
export(simulate_matrix)
export(varls)
export(vars.procmod)
export(weighted.residuals)
import(MASS)
import(expm)
R/anova.pm.R deleted 100644 → 0
#' @include mprocuste.R
#' Anova function for `pm` objects.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
anova.pm = function(m, permutations = how(nperm = 999),...) {
mf = m$model
coef = m$coefficients
sdxy = sqrt(diag(m$cov))
mt = terms(mf)
reponse = attr(mt,"response")
explicatives = attr(mt,'term.labels')
nXs = length(explicatives)
coef.norm = coef * sdxy[explicatives] / sdxy[reponse]
cors = m$cov / outer(sdxy,sdxy)
#print(reponse)
#print(explicatives)
cors.Y = cors[reponse,explicatives]
SCT = m$SCT
SCR = sum(m$residuals^2)
SCE = SCT * abs(coef.norm * cors.Y)
names(SCE) = explicatives
SCET = sum(SCE)
SCIT = SCT - SCET - SCR
SumSq = c(SCE,Residuals=SCR)
Df = c(rep(1,nXs),m$df.residual)
MeanSq = SumSq/Df
Fvalue = c(MeanSq[-length(MeanSq)]/MeanSq[length(MeanSq)],NA)
Pvalue = 1-pf(Fvalue,
Df[-length(MeanSq)],
Df[length(MeanSq)])
result = data.frame(Df,
SumSq,
MeanSq,
Fvalue,
Pvalue)
colnames(result)=c("Df","Sum Sq","Mean Sq","F value","Pr(>F)")
attr(result, "heading")= paste("Analysis of Variance Table\n\nResponse:",
all.vars(m$terms)[1]
)
class(result)=c("anova",class(result))
return(result)
}
#' Generate permutation matrix according to a schema.
#'
#' The permutation schema is defined using the `how` function.
#' The implementation of this function is currectly extraced
#' from the VEGAN package and just copied pasted here to avoid
#' dependency on an hidden vegan function.
#'
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
}
R/covls.R
#' @include procmod.frame.R
#' @include procmod_frame.R
#' @include multivariate.R
#'
#' @import expm
......@@ -9,13 +9,12 @@ 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)
#' m <- matrix(1:16, nrow = 4)
#' ProcMod:::.Trace(m)
#'
#' @note Internal function do not use.
#'
......@@ -23,7 +22,7 @@ NULL
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#'
.Trace = function(X) sum(diag(X))
.Trace <- function(X) sum(diag(X))
#' Compute the variance, covariance matrix of K coordinate matrices.
#'
......@@ -36,111 +35,114 @@ NULL
#' @param rcovls
#'
#' @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)
#'
#' # 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
varls = function(...,nperm = 100,rcovls=FALSE) {
varls <- function(..., nperm = 100, rcovls = FALSE) {
if (is.numeric(nperm) &&
length(nperm)==1 &&
nperm==0)
nperm=NULL
length(nperm) == 1 &&
nperm == 0) {
nperm <- NULL
}
if (is.null(nperm))
rcovls=FALSE
if (is.null(nperm)) {
rcovls <- FALSE
}
Xs <- list(...)
if (length(Xs)==1) {
x = Xs[[1]]
if (is.procmod.frame(x))
Xs=x
else if (is.pm(x))
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 <- procmod.frame(x)
}
}
else {
Xs <- as.procmod.frame(Xs)
}
else
Xs=as.procmod.frame(Xs)
Xnames=names(Xs)
x_names <- names(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) - 1
CovXXs = matrix(0, nrow = nX, ncol = nX)
cov_xxs <- matrix(0, nrow = nx, ncol = nx)
# Computes the standard Covls covariance matrix
for (i in seq_len(nX))
for (j in i:nX) {
CovXXs[i,j] = Re(.Trace(sqrtm(XXs[[i]] %*% XXs[[j]])))
for (i in seq_len(nx))
for (j in i:nx) {
cov_xxs[i, j] <- Re(.Trace(sqrtm(XXs[[i]] %*% XXs[[j]])))
}
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))
CovXX2s = matrix(0, nrow = nX, ncol = nX)
if (! is.null(nperm)) {
rCovXXs = array(0,dim=c(nperm,nX,nX))
vXs=vector(mode="list",nX)
lXs=vapply(Xs,FUN = function(x) prod(dim(x)),0)
dXs=vapply(Xs,FUN = function(x) dim(ortho(x)),c(0,0))
for (i in seq_len(nX))
vXs[i]=list(t(sqrt(colMeans(Xs[i]^2) - colMeans(Xs[i])^2) %*% t(rep(1,nrow(Xs[i])))))
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])))))
for (k in seq_len(nperm)) {
# Generate the random matrix for the kth simulation
rXs=vector(mode="list",nX)
for (i in seq_len(nX))
rXs[i]= list(tcrossprod(scale(array(rnorm(lXs[i]),
dim = dXs[,i])
) * vXs[[i]]))
for (i in seq_len(nX))
for (j in i:nX)
rCovXXs[k,i,j] = .Trace(sqrtm(rXs[[i]] %*% rXs[[j]]))
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]]))
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]]))
}
}
rCovXXs = apply(rCovXXs,
MARGIN = c(2,3),
FUN=function(i) MASS::fitdistr(Re(i),
"normal")$estimate['mean']
)
for (i in seq_len(nX))
for (j in i:nX) {
CovXXs[i,j] = max(CovXXs[i,j] - rCovXXs[i,j],0)
CovXXs[j,i] = CovXXs[i,j]
rCovXXs[j,i] = rCovXXs[i,j]
r_cov_xxs <- apply(r_cov_xxs,
MARGIN = c(2, 3),
FUN = function(i) MASS::fitdistr(
Re(i),
"normal"
)$estimate["mean"]
)
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]
}
CovXXs = CovXXs / N
colnames(CovXXs)=Xnames
rownames(CovXXs)=Xnames
cov_xxs <- cov_xxs / N
colnames(cov_xxs) <- x_names
rownames(cov_xxs) <- x_names
CovXXs = Re(CovXXs)
cov_xxs <- Re(cov_xxs)
if (rcovls)
attr(CovXXs,"rcovls") = rCovXXs / N
if (rcovls) {
attr(cov_xxs, "rcovls") <- r_cov_xxs / N
}
attr(CovXXs,"nperm") = nperm
attr(cov_xxs, "nperm") <- nperm
return(make_subS3Class(Re(CovXXs), "procmod.varls"))
return(make_subS3Class(Re(cov_xxs), "procmod.varls"))
}
......@@ -149,16 +151,17 @@ 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)
s = sqrt(diag(cov))
vv= outer(s,s)
corls <- function(..., nperm = 100, rcovls = FALSE) {
cov <- varls(..., nperm = nperm, rcovls = rcovls)
s <- sqrt(diag(cov))
vv <- outer(s, s)
rls = make_subS3Class(cov/vv, "procmod.corls")
attr(rls,"nperm")=attr(cov,"nperm")
rls <- make_subS3Class(cov / vv, "procmod.corls")
attr(rls, "nperm") <- attr(cov, "nperm")
if (rcovls)
attr(rls,"rcovls") = attr(rls,"rcovls") / vv
if (rcovls) {
attr(rls, "rcovls") <- attr(rls, "rcovls") / vv
}
return(rls)
}
......@@ -168,14 +171,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)
C = solve(rls)
D = sqrt(diag(C) %o% diag(C))
rp = make_subS3Class(C/D,"procmod.corls")
attr(rp,"nperm")=attr(rls,"nperm")
if (rcovls)
attr(rp,"rcovls") = attr(rls,"rcovls")
corls.partial <- function(..., nperm = 100, rcovls = FALSE) {
rls <- corls(..., nperm = nperm, rcovls = rcovls)
C <- solve(rls)
D <- sqrt(diag(C) %o% diag(C))
rp <- make_subS3Class(C / D, "procmod.corls")
attr(rp, "nperm") <- attr(rls, "nperm")
if (rcovls) {
attr(rp, "rcovls") <- attr(rls, "rcovls")
}
return(rp)
}
......@@ -185,10 +191,10 @@ corls.partial = function(...,nperm = 100,rcovls=FALSE) {
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
print.procmod.varls = function (x, ...) {
class(x)="matrix"
attr(x,"permutations")=NULL
attr(x,"rcovls") =NULL
print.procmod.varls <- function(x, ...) {
class(x) <- "matrix"
attr(x, "permutations") <- NULL
attr(x, "rcovls") <- NULL
print(x)
}
......@@ -197,9 +203,9 @@ print.procmod.varls = function (x, ...) {
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
print.procmod.corls = function (x, ...) {
class(x)="matrix"
attr(x,"permutations")=NULL
attr(x,"rcovls") =NULL
print.procmod.corls <- function(x, ...) {
class(x) <- "matrix"
attr(x, "permutations") <- NULL
attr(x, "rcovls") <- NULL
print(x)
}
R/formula.procmod.frame.R deleted 100644 → 0
#' @include internals.R
NULL
#' Build a formula from a procmod.frame.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
formula.procmod.frame = function (x, ...)
{
nm <- sapply(names(x), as.name)
if (length(nm) > 1L) {
rhs <- nm[-1L]
lhs <- nm[1L]
}
else if (length(nm) == 1L) {
rhs <- nm[1L]
lhs <- NULL
}
else stop("cannot create a formula from an empty procmod.frame")
ff <- parse(text = paste(lhs,
paste(rhs,
collapse = "+"),
sep = "~"),
keep.source = FALSE)
ff <- eval(ff)
environment(ff) <- parent.frame()
ff
}
R/internals.R
make_subS3Class = function(obj,subclass) {
class(obj) = c(paste(subclass,
collapse = "_"),
class(obj))
make_subS3Class <- function(obj, subclass) {
class(obj) <- c(
paste(subclass,
collapse = "_"
),
class(obj)
)
return(obj)
}
dots.names=function(...) {
varnames = substitute(list(...))[-1L]
dots = list(...)
isname = sapply(varnames,is.name)
charname = as.character(varnames)
charname[!isname]=""
dots.names <- function(...) {
varnames <- substitute(list(...))[-1L]
dots <- list(...)
isname <- sapply(varnames, is.name)
charname <- as.character(varnames)
charname[!isname] <- ""
n=length(dots)
n <- length(dots)
explicit = names(dots)
explicit <- names(dots)
if (is.null(explicit))
explicit=character(n)
if (is.null(explicit)) {
explicit <- character(n)
}
ze = !nzchar(explicit)
ze <- !nzchar(explicit)
explicit[ze]=charname[ze]
ze = !nzchar(explicit)
explicit[ze] <- charname[ze]
ze <- !nzchar(explicit)
dnames <- paste('V',seq_len(n),sep='')
explicit[ze]=dnames[ze]
dnames <- paste("V", seq_len(n), sep = "")
explicit[ze] <- dnames[ze]
return(explicit)
}
make_procmod_subS3Class = function(obj,subclass) {
class(obj) = c(paste("procmod",subclass,
sep="_",collapse = "_"),
class(obj))
make_procmod_subS3Class <- function(obj, subclass) {
class(obj) <- c(
paste("procmod", subclass,
sep = "_", collapse = "_"
),
class(obj)
)
return(obj)
}
make_procmod_data = function(obj,subclass) {
make_procmod_data <- function(obj, subclass) {
eud <- inherits(obj, "procmod_data", which = TRUE)
eud = inherits(obj,'procmod_data',which = TRUE)
if (eud > 0) {
class(obj) <- class(obj)[-1:-(eud - 1)]
} else {
obj <- make_procmod_subS3Class(obj, "data")
}
if (eud > 0)
class(obj) = class(obj)[-1:-(eud-1)]
else
obj = make_procmod_subS3Class(obj,'data')
if (! missing(subclass))
obj = make_procmod_subS3Class(obj,subclass)
if (!missing(subclass)) {
obj <- make_procmod_subS3Class(obj, subclass)
}
return(obj)
}
R/mprocuste.R deleted 100644 → 0
#' @include internals.R
#' @include model.procmod.R
#'
#' @title ProcMod
#' @description blabla
#' @author Christelle Gonindard-Melodelima
#'
NULL
pm.fit = function(covmat,y,xs,
singular.ok = singular.ok,
tol) {
xy.cov = covmat[xs,y,drop=FALSE]
xx.cov = covmat[xs,xs,drop=FALSE]
#xx.qr = qr(xx.cov,tol=tol)
#if (!singular.ok && xx.qr$rank < ncol(xx.cov))
# stop("singular fit encountered")
coef = as.numeric(solve(xx.cov) %*% xy.cov)
names(coef)=rownames(xx.cov)
#coef = qr.coef(xx.qr,xy.cov)
#effects= qr.qty(xx.qr,xy.cov)
return(list(coefficients=coef
#rank=xx.qr$rank + 1
#qr=xx.qr,
)
)
}
.ctrace <- function(MAT) sum(MAT^2)
#' Tests that the object is as `pm` instance.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
is.pm = function(obj) {
inherits(obj, "pm")
}
#' Performs a procruste model on a set of coordinate matrices
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
pm = function (formula,data, subset, weights, na.action, method = "qr",
singular.ok = TRUE, tol = 1e-07,
model = TRUE, x = TRUE, y = TRUE, A = TRUE,
nperm = 100) {
ret.x <- x
ret.y <- y
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights",
"na.action"),
table = names(mf),
nomatch = 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- quote(ProcMod::model.procmod.default)
mf <- eval(mf, parent.frame())
if (method == "model.frame")
return(mf)
else if (method != "qr")
warning(gettextf("method = '%s' is not supported. Using 'qr'",
method), domain = NA)
mt <- attr(mf, "terms")
Y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
if (!is.null(w) && !is.numeric(w))
stop("'weights' must be a numeric vector")
if (is.empty.model(mt)) {
x <- NULL
z <- list(coefficients = NULL,
anova = NULL,
SCT = SCT,