R/subniche.R
In polak51/WitOMI: Within Outlying Mean Indexes: Refining the OMI Analysis
Defines functions
subniche
subparam.refor
rtestrefor
subparam.subor
rtestsubor
Documented in
rtestrefor
rtestsubor
subniche
subparam.refor
subparam.subor
#' @import ade4
#' @title The Within Outlying Mean Indexes calculation
#' @aliases subniche print.subniche summary.subniche refparam rtest.subniche subparam.refor rtestrefor subparam.subor rtestsubor
#' @description The indexes allows to divide the niche, estimated from the \link[ade4]{niche} function in the \link{ade4} package into subniches defined by a factor, which creates the
#' subsets. See details for more information.
#' @usage subniche(nic, factor)
#' @param nic an object of class \code{niche}.
#' @param factor a factor which will defined the subsets within which the subniches will be calculated (the same length of the number of sites)
#' @param x an object of class \code{subniche}.
#' @param xtest an object of class \code{subniche}.
#' @param object an object of class \code{subniche}.
#' @param ... further arguments passed to or from other methods
#' @param sim a numeric vector of simulated values
#' @param obs a numeric vector of an observed value
#' @param alter a character string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two-sided".he length must be equal to the length of the vector obs, values are recycled if shorter.
#' @param names a vector of names for tests
#' @param subpvalue the subset pvalue resulting from \code{subkrandtest} function
#' @param p.adjust.method a string indicating a method for multiple adjustment, see \link[stats]{p.adjust.methods} for possible choices.
#' @param call a call order
#' @param nrepet the number of permutations for the testing procedure
#' @return Adds items in the niche list and changing the class into \code{subniche} containing:
#' @return \code{factor} the factor use to divide the environmental and species matrix into submatrices.
#' @return \code{G_k} a dataframe with the sub-origins, \emph{G_k}.
#' @return \code{sub} a dataframe with the species subniche coordinates
#' @keywords subniche
#' @author Stephane Karasiewicz, \email{stephane.karasiewicz@wanadoo.fr}
#' @details The Within Outlying Mean Index analysis is a statistical exploratory niche analysis which provides observation of niche shift and/or conservatism, of an entire community,at different subcales
#' (temporal ,spatial and/or finer biological organisation level), and comparable under the same environmental gradients. This hindcasting multivariate analysis is based on the OMI analysis (Doledec \emph{et al.} 2000) which is used as reference.
#' The niches refinement is inspired by the K-select (Calenge \emph{et al.} 2005) which emphasizes the limiting factors in habitat use in design II and III (Thomas and Taylor, 1990).The different estimations should help understand:
#'
#' 1. the environmental factors defining a species' reference niche, under on the full scale, within a community.
#'
#' 2. the environmental factors defining a species' subniches, under each subsets, within a community.
#'
#' The subniches parameters can be calculated from both the reference origin,\emph{G}, which corresponds to the reference plan origin, and from \emph{G_k}, which corresponds to the
#' suborigins. \emph{G} is the graphical representation of the mean environmental conditions encountered over the full scale of the data. \emph{G_k} is the mean environmental conditions
#' encountered at a subset defined by the factor. They are complementary has you can compare:
#'
#' 1. a single species' subniches to \emph{G}.
#'
#' 2. the community' subniches to \emph{G_k} at a specific subset.
#'
#' The subniches of a single species can only be compared to \emph{G} as it is the common origin to all subsets. Whereas \emph{G_k} is only common to the species found within the
#' subset. So comparing different subniches of one species, found within different subsets, is only relevant to \emph{G}. The community's subniches can be compared to both \emph{G}
#' and \emph{G_k}, but \emph{G}, being the mean environmental conditions found within the full scale, will not express the specificity of the environmental conditions that the species
#' encountered at the subset. \emph{G_k}, being the mean environmental conditions of the subset, will reflect the atypical value of the environmental condition, making the
#' comparison of the community's subniches parameters more relevant. More information on the ecological concept can be found in Karasiewicz \emph{et al.} 2017.
#'
#' For more details description on the package use:\url{https://github.com/KarasiewiczStephane/WitOMI}.
#'
#' @references Karasiewicz S.,Doledec S.and Lefebvre S. (2017). Within outlying mean indexes: refining the OMI analysis for the realized niche decomposition. \emph{PeerJ} 5:e3364. \url{https://doi.org/10.7717/peerj.3364}.
#'
#' Calenge C., Dufour A.B. and Maillard D. (2005). K-select analysis: a new method to analyse habitat selection in radio-tracking studies. \emph{Ecological modelling}, \bold{186}, 143-153.
#'
#' Doledec S., Chessel D. and Gimaret C. (2000). Niche separation in community analysis: a new method. \emph{Ecology},\bold{81}, 2914-1927.
#'
#' Thomas, D.L., Taylor, E.J. (1990). Study Designs and Tests for Comparing Resource Use and Availability II. \emph{Natl. Widl.} \bold{54}(2), 322-330.
#'
#' @seealso \link[ade4]{niche} \link[ade4]{niche.param}
#' @examples
#'library(subniche)
#'data(doubs)
#'dudi1 <- dudi.pca(doubs$env, scale = TRUE, scan = FALSE, nf = 3)
#'nic1 <- niche(dudi1, doubs$fish, scann = FALSE)
#'# number of sites
#'N <- dim(nic1$ls)[1]
#'#Create a factor which defines the subsets
#'fact <- factor(c(rep(1,N/2),rep(2,N/2)))
#'# nic1 will be use as reference and fact will be use to define the subniches environment
#'subnic1 <- subniche(nic1, fact)
#'# the following two functions do the same display, plot.refniche is adapted to subniche objects
#'plot(nic1)
#'plot(subnic1)
#'#Display the marginality vector of the suborigins and the species subniche
#'margvect(subnic1)
#'#Display the subset's polygon, found within the overall environment's chull,
#'#and the corresponding species positions
#'subplot(subnic1)
#'# The following two functions do the same display, refparam is adapted to subniche objects
#'niche.param(nic1)
#'refparam(subnic1)
#'# The following two functions do the same display, rtest is adapted to subniche objects
#'rtest(nic1,10)
#'rtest(subnic1,10)
#'#Calculates the subniches' parameters from G with the corresponding rtest
#'subparam.refor(subnic1)
#'rtestrefor(subnic1,10)
#'#Calculates the subniches' parameters from G_k with the corresponding rtest
#'subparam.subor(subnic1)
#'rtestsubor(subnic1,10)
#' @export subniche
#' @rdname subniche
#' @importFrom stats weighted.mean
subniche <- function(nic, factor){
if (!inherits(nic, "niche"))
stop("Object of class niche expected")
appel <- as.list(nic$call)
Y <- eval.parent(appel[[3]])
liani <- split(as.data.frame(nic$ls), factor)
liwei <- sep.factor.row(Y,factor)
for(i in 1:length(liwei)){
w <- liwei[[i]]
w2 <- apply(w, 2, sum)
liwei[[i]] <- sweep(w, 2, w2, "/")
}
G_k <- as.data.frame(t(as.matrix(data.frame(lapply(liani,
function(x) apply(x, 2, mean))))))
names(G_k) <- names(nic$li)
mutemp <- list()
for (i in 1:length(liwei)){
mutemp[[i]] <- matrix(nrow=length(colnames(Y)), ncol=nic$nf)
rownames(mutemp[[i]]) <- paste(colnames(Y),levels(factor)[[i]],sep="")
for(j in 1:length(colnames(Y))){
z <- liwei[[i]][,j]
mutemp[[i]][j,] <- apply(liani[[i]],2, weighted.mean, z)
}
}
mutemp <- do.call("rbind", mutemp)
nic$sub <- mutemp
nic$G_k <- G_k
rownames(nic$G_k) <- paste(rep("G_k",dim(nic$G_k)[1]),levels(factor),sep="")
nic$factor <- factor
class(nic) <- c("subniche", "dudi")
return(nic)
}
#' @rdname subniche
#' @method print subkrandtest
#' @export
print.subkrandtest <- function (x, ...)
{
if (!inherits(x, "subkrandtest"))
stop("to be used with 'subkrandtest' object")
cat("class:", class(x), "\n")
cat("Monte-Carlo tests\n")
cat("Call: ")
print(x$call)
cat("\nNumber of tests: ", x$ntest, "\n")
cat("\nAdjustment method for multiple comparisons: ", x$adj.method,
"\n")
sumry <- list(Test = x$names, Obs = x$obs, Std.Obs = x$expvar[,
1], Alter = x$alter)
sumry <- as.data.frame(sumry)
row.names(sumry) <- 1:x$ntest
if (any(x$rep[1] != x$rep)) {
sumry <- cbind(sumry[, 1:4], N.perm = x$rep)
}
else {
cat("Permutation number: ", x$rep[1], "\n")
}
sumry <- cbind(sumry, Pvalue = x$pvalue)
if (x$adj.method != "none")
sumry <- cbind(sumry, Pvalue.adj = x$adj.pvalue)
print(sumry)
cat("Subsets Pvalue:", x$subpvalue, "\n")
if (length(names(x)) > 10) {
cat("other elements: ")
cat(names(x)[11:(length(x))], "\n")
}
}
#' @rdname subniche
#' @method print subnikrandtest
#' @export
print.subnikrandtest <- function (x, ...)
{
if (!inherits(x, "subnikrandtest"))
stop("to be used with 'subnikrandtest' object")
cat("class:", class(x), "\n")
cat("Monte-Carlo tests\n")
cat("Call: ")
print(x$call)
cat("\nNumber of tests: ", x$ntest, "\n")
cat("\nAdjustment method for multiple comparisons: ", x$adj.method,
"\n")
sumry <- list(Test = x$names, Obs = x$obs, Std.Obs = x$expvar[,
1], Alter = x$alter)
sumry <- as.data.frame(sumry)
row.names(sumry) <- 1:x$ntest
if (any(x$rep[1] != x$rep)) {
sumry <- cbind(sumry[, 1:4], N.perm = x$rep)
}
else {
cat("Permutation number: ", x$rep[1], "\n")
}
sumry <- cbind(sumry, Pvalue = x$pvalue, SubniPvalue = x$subni.pvalue)
if (x$adj.method != "none")
sumry <- cbind(sumry, Pvalue.adj = x$adj.pvalue, SubniPvalue = x$subni.pvalue)
print(sumry)
cat("\n")
if (length(names(x)) > 9) {
cat("other elements: ")
cat(names(x)[10:(length(x))], "\n")
}
}
#' @rdname subniche
#' @method print subniche
#' @export
print.subniche <- function (x, ...)
{
if (!inherits(x, "subniche"))
stop("to be used with 'niche' object")
cat("WitOMI calculation\n")
cat("call: ")
print(x$call)
cat("class: ")
cat(class(x), "\n")
cat("\n$rank (rank) :", x$rank)
cat("\n$nf (axis saved) :", x$nf)
cat("\n$RV (RV coeff) :", x$RV)
cat("\n\neigen values: ")
l0 <- length(x$eig)
cat(signif(x$eig, 4)[1:(min(5, l0))])
if (l0 > 5)
cat(" ...\n\n")
else cat("\n\n")
sumry <- array("", c(4, 4), list(1:4, c("vector", "length",
"mode", "content")))
sumry[1, ] <- c("$eig", length(x$eig), mode(x$eig), "eigen values")
sumry[2, ] <- c("$lw", length(x$lw), mode(x$lw), "row weigths (crossed array)")
sumry[3, ] <- c("$cw", length(x$cw), mode(x$cw), "col weigths (crossed array)")
sumry[4, ] <- c("$factor", length(x$factor), mode(x$factor), "factor used for creating subsets")
print(sumry, quote = FALSE)
cat("\n")
sumry <- array("", c(9, 4), list(1:9, c("data.frame", "nrow",
"ncol", "content")))
sumry[1, ] <- c("$tab", nrow(x$tab), ncol(x$tab), "crossed array (averaging species/sites)")
sumry[2, ] <- c("$li", nrow(x$li), ncol(x$li), "species coordinates")
sumry[3, ] <- c("$l1", nrow(x$l1), ncol(x$l1), "species normed scores")
sumry[4, ] <- c("$co", nrow(x$co), ncol(x$co), "variables coordinates")
sumry[5, ] <- c("$c1", nrow(x$c1), ncol(x$c1), "variables normed scores")
sumry[6, ] <- c("$ls", nrow(x$ls), ncol(x$ls), "sites coordinates")
sumry[7, ] <- c("$as", nrow(x$as), ncol(x$as), "axis upon niche axis")
sumry[8, ] <- c("$G_k", nrow(x$G_k), ncol(x$G_k), "G_k coordinates")
sumry[9, ] <- c("$sub", nrow(x$sub), ncol(x$sub), "species coordinates within each subset")
print(sumry, quote = FALSE)
cat("\n")
}
#' @rdname subniche
#' @method summary subniche
#' @export
summary.subniche <- function (object, ...)
{
cat("Class: ")
cat(class(object))
cat("\nCall: ")
print(object$call)
cat("\nTotal inertia: ")
cat(signif(sum(object$eig), 4))
cat("\n")
l0 <- length(object$eig)
cat("\nEigenvalues:\n")
vec <- object$eig[1:(min(5, l0))]
names(vec) <- paste("Ax", 1:length(vec), sep = "")
print(format(vec, digits = 4, trim = TRUE, width = 7), quote = FALSE)
cat("\nProjected inertia (%):\n")
vec <- (object$eig/sum(object$eig) * 100)[1:(min(5, l0))]
names(vec) <- paste("Ax", 1:length(vec), sep = "")
print(format(vec, digits = 4, trim = TRUE, width = 7), quote = FALSE)
cat("\nCumulative projected inertia (%):\n")
vec <- (cumsum(object$eig)/sum(object$eig) * 100)[1:(min(5,
l0))]
names(vec)[1] <- "Ax1"
if (l0 > 1)
names(vec)[2:length(vec)] <- paste("Ax1:", 2:length(vec),
sep = "")
print(format(vec, digits = 4, trim = TRUE, width = 7), quote = FALSE)
if (l0 > 5) {
cat("\n")
cat(paste("(Only 5 dimensions (out of ", l0, ") are shown)\n",
sep = "", collapse = ""))
}
cat("\n")
}
#' @rdname subniche
#' @export
refparam <- function (x)
{
if (!inherits(x, "subniche"))
stop("Object of class 'subniche' expected")
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab
Y <- eval.parent(appel[[3]])
w1 <- apply(Y, 2, sum)
if (any(w1 <= 0))
stop(paste("Column sum <=0 in Y"))
Y <- sweep(Y, 2, w1, "/")
calcul.param <- function(freq, mil) {
inertia <- sum(freq * mil * mil)
m <- apply(freq * mil, 2, sum)
margi <- sum(m^2)
mil <- t(t(mil) - m)
tolt <- sum(freq * mil * mil)
u <- m/sqrt(sum(m^2))
z <- mil %*% u
tolm <- sum(freq * z * z)
tolr <- tolt - tolm
w <- c(inertia, margi, tolm, tolr)
names(w) <- c("inertia", "OMI", "Tol", "Rtol")
w1 <- round(w[2:4]/w[1], digits = 3) * 100
names(w1) <- c("omi", "tol", "rtol")
return(c(w, w1))
}
res <- apply(Y, 2, calcul.param, mil = X)
t(res)
}
#' @rdname subniche
#' @method rtest subniche
#' @export
rtest.subniche <- function (xtest, nrepet = 99, ...)
{
if (!inherits(xtest, "dudi"))
stop("Object of class dudi expected")
if (!inherits(xtest, "subniche"))
stop("Type 'niche' expected")
appel <- as.list(xtest$call)
X <- eval.parent(appel$dudiX)$tab
Y <- eval.parent(appel$Y)
w1 <- apply(Y, 2, sum)
if (any(w1 <= 0))
stop(paste("Column sum <=0 in Y"))
Y <- sweep(Y, 2, w1, "/")
calcul.margi <- function(freq, mil) {
m <- apply(freq * mil, 2, sum)
return(sum(m^2))
}
obs <- apply(Y, 2, calcul.margi, mil = X)
obs <- c(obs, OMI.mean = mean(obs))
sim <- sapply(1:nrepet, function(x) apply(apply(Y, 2, sample),
2, calcul.margi, mil = X))
sim <- rbind(sim, OMI.mean = apply(sim, 2, mean))
res <- as.krandtest(obs = obs, sim = t(sim))
return(res)
}
#' @rdname subniche
#' @export
subparam.refor <- function(x){
res <- list()
appel <- as.list(x$call)
y <-eval.parent(appel[[3]])
factor <- x$factor
subniche.param <- function (x,y)
{
if (!inherits(x, "subniche"))
stop("Object of class 'subniche' expected")
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab[y,]
Y <- eval.parent(appel[[3]])[y,]
w1 <- apply(Y, 2, sum)
Y <- sweep(Y, 2, w1, "/")
calcul.param <- function(freq, mil) {
inertia <- sum(freq * mil * mil)
m <- apply(freq * mil, 2, sum)
margi <- sum(m^2)
mil <- t(t(mil) - m)
tolt <- sum(freq * mil * mil)
u <- m/sqrt(sum(m^2))
z <- mil %*% u
tolm <- sum(freq * z * z)
tolr <- tolt - tolm
w <- c(inertia, margi, tolm, tolr)
names(w) <- c("inertia", "WitOMIG", "Tol", "Rtol")
w1 <- round(w[2:4]/w[1], digits = 3) * 100
names(w1) <- c("witomig", "tol", "rtol")
return(c(w, w1))
}
res <- apply(Y, 2, calcul.param, mil = X)
t(res)
}
N <-length(levels(x$factor))
nam <- levels(factor)
for(i in 1:N){
res[[i]] <- subniche.param(x,array(factor)==nam[i])
rownames(res[[i]]) <- paste(rownames(x$li),nam[i],sep="")
}
res <- do.call("rbind",res)
return(res)
}
#' @rdname subniche
#' @export
rtestrefor <- function(x, nrepet){
factor <- x$factor
res <- list()
N <-length(levels(x$factor))
nam <- levels(factor)
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab
Y <- eval.parent(appel[[3]])
calcul.margi <- function(freq, mil) {
m <- apply(freq * mil, 2, sum)
return(sum(m^2))
}
for(i in 1:N){
X1 <- X[array(factor)==nam[i],]
Xwobs <- apply(X1, 2, mean)
Xsim <- sapply(1:nrepet, function(x) apply(t <- apply(X, 2, sample)[array(factor)==nam[i],], 2, mean))
Xtest <- subkrandtest(obs=Xwobs, sim=t(Xsim), "two-sided")
w1 <- apply(Y[array(factor)==nam[i],], 2, sum)
Y1 <- sweep(Y[array(factor)==nam[i],], 2, w1, "/")
obs <- apply(Y1, 2, calcul.margi, mil = X1)
obs <- c(obs, OMI.mean = mean(obs,na.rm=TRUE))
sim <- sapply(1:nrepet, function(x) apply(sweep(t <- apply(Y, 2, sample)[array(factor)==nam[i],], 2,apply(t, 2, sum),"/"), 2, calcul.margi, mil = X1))
sim <- rbind(sim, OMI.mean = apply(sim, 2, mean,na.rm=TRUE))
omitest <- subnikrandtest(obs = obs, sim = t(sim),subpvalue= Xtest$subpvalue)
res[[i]] <- list("Subsettest" = Xtest, "witomigtest" = omitest)
}
names(res) <- nam
return(res)
}
#' @rdname subniche
#' @export
subparam.subor <- function(x){
res <- list()
appel <- as.list(x$call)
y <-eval.parent(appel[[3]])
factor <- x$factor
subnichesub.param <- function (x,y)
{
if (!inherits(x, "subniche"))
stop("Object of class 'subniche' expected")
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab[y,]
X <- as.matrix(scale(X, center=T,scale=F))
Y <- eval.parent(appel[[3]])[y,]
w1 <- apply(Y, 2, sum)
Y <- sweep(Y, 2, w1, "/")
calcul.param <- function(freq, mil) {
inertia <- sum(freq * mil * mil)
m <- apply(freq * mil, 2, sum)
margi <- sum(m^2)
mil <- t(t(mil) - m)
tolt <- sum(freq * mil * mil)
u <- m/sqrt(sum(m^2))
z <- mil %*% u
tolm <- sum(freq * z * z)
tolr <- tolt - tolm
w <- c(inertia, margi, tolm, tolr)
names(w) <- c("inertia", "WitOMIG_k", "Tol", "Rtol")
w1 <- round(w[2:4]/w[1], digits = 3) * 100
names(w1) <- c("witomig_k", "tol", "rtol")
return(c(w, w1))
}
res <- apply(Y, 2, calcul.param, mil = X)
t(res)
}
N <-length(levels(x$factor))
nam <- levels(factor)
for(i in 1:N){
res[[i]] <- subnichesub.param(x,array(factor)==nam[i])
rownames(res[[i]]) <- paste(rownames(x$li),i,sep="")
}
res <- do.call("rbind",res)
return(res)
}
#' @rdname subniche
#' @export
rtestsubor <- function(x, nrepet){
factor <- x$factor
res <- list()
N <-length(levels(x$factor))
nam <- levels(factor)
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab
Y <- eval.parent(appel[[3]])
calcul.margi <- function(freq, mil) {
m <- apply(freq * mil, 2, sum)
return(sum(m^2))
}
for(i in 1:N){
X1 <- X[array(factor)==nam[i],]
Xwobs <- apply(X1, 2, mean)
Xobs <- sweep(X1, 2,Xwobs,"-")
Xsim <- sapply(1:nrepet, function(x) apply(t <- apply(X, 2, sample)[array(factor)==nam[i],], 2, mean))
Xtest <- subkrandtest(obs=Xwobs, sim=t(Xsim), "two-sided")
w1 <- apply(Y[array(factor)==nam[i],], 2, sum)
Y1 <- sweep(Y[array(factor)==nam[i],], 2, w1, "/")
obs <- apply(Y1, 2, calcul.margi, mil = Xobs)
obs <- c(obs, OMI.mean = mean(obs,na.rm=TRUE))
sim <- sapply(1:nrepet, function(x) apply(sweep(t <- apply(Y, 2, sample)[array(factor)==nam[i],], 2,apply(t, 2, sum),"/"), 2, calcul.margi, mil = Xobs))
sim <- rbind(sim, OMI.mean = apply(sim, 2, mean,na.rm=TRUE))
omitest <- subnikrandtest(obs = obs, sim = t(sim),subpvalue=Xtest$subpvalue)
res[[i]] <- list("Subsettest" = Xtest, "witomig_ktest" = omitest)
}
names(res) <- nam
return(res)
}
#' @rdname subniche
#' @export
#' @importFrom stats sd na.omit p.adjust p.adjust.methods
subkrandtest<-function (sim, obs, alter = "greater", call = match.call(), names = colnames(sim),
p.adjust.method = "none")
{
res <- list(sim = sim, obs = obs)
if (length(obs) != length(alter))
alter <- rep(alter, length = length(obs))
res$alter <- alter
res$rep <- apply(sim, 2, function(x) length(na.omit(x)))
res$ntest <- length(obs)
res$expvar <- data.frame(matrix(0, res$ntest, 3))
if (!is.null(names)) {
res$names <- names
}
else {
res$names <- paste("test", 1:res$ntest, sep = "")
}
names(res$expvar) <- c("Std.Obs", "Expectation", "Variance")
res$pvalue <- rep(0, length(obs))
for (i in 1:length(obs)) {
vec.sim <- na.omit(sim[, i])
res$alter[i] <- match.arg(res$alter[i], c("greater",
"less", "two-sided"))
res$expvar[i, 1] <- (res$obs[i] - mean(vec.sim))/sd(vec.sim)
res$expvar[i, 2] <- mean(vec.sim)
res$expvar[i, 3] <- sd(vec.sim)
if (res$alter[i] == "greater") {
res$pvalue[i] <- (sum(vec.sim >= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "less") {
res$pvalue[i] <- (sum(vec.sim <= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "two-sided") {
sim0 <- abs(vec.sim - mean(vec.sim))
obs0 <- abs(obs[i] - mean(vec.sim))
res$pvalue[i] <- (sum(sim0 >= obs0) + 1)/(res$rep[i] +
1)
}
}
p.adjust.method <- match.arg(p.adjust.method, p.adjust.methods)
res$adj.pvalue <- p.adjust(res$pvalue, method = p.adjust.method)
res$adj.method <- p.adjust.method
res$subpvalue <- prod(res$pvalue)
res$call <- call
class(res) <- "subkrandtest"
return(res)
}
#' @rdname subniche
#' @export
#' @importFrom stats sd na.omit p.adjust p.adjust.methods
subnikrandtest<-function (sim, obs, alter = "greater",subpvalue, call = match.call(), names = colnames(sim),
p.adjust.method = "none")
{
res <- list(sim = sim, obs = obs)
if (length(obs) != length(alter))
alter <- rep(alter, length = length(obs))
res$alter <- alter
res$rep <- apply(sim, 2, function(x) length(na.omit(x)))
res$ntest <- length(obs)
res$expvar <- data.frame(matrix(0, res$ntest, 3))
if (!is.null(names)) {
res$names <- names
}
else {
res$names <- paste("test", 1:res$ntest, sep = "")
}
names(res$expvar) <- c("Std.Obs", "Expectation", "Variance")
res$pvalue <- rep(0, length(obs))
for (i in 1:length(obs)) {
vec.sim <- na.omit(sim[, i])
res$alter[i] <- match.arg(res$alter[i], c("greater",
"less", "two-sided"))
res$expvar[i, 1] <- (res$obs[i] - mean(vec.sim))/sd(vec.sim)
res$expvar[i, 2] <- mean(vec.sim)
res$expvar[i, 3] <- sd(vec.sim)
if (res$alter[i] == "greater") {
res$pvalue[i] <- (sum(vec.sim >= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "less") {
res$pvalue[i] <- (sum(vec.sim <= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "two-sided") {
sim0 <- abs(vec.sim - mean(vec.sim))
obs0 <- abs(obs[i] - mean(vec.sim))
res$pvalue[i] <- (sum(sim0 >= obs0) + 1)/(res$rep[i] +
1)
}
}
p.adjust.method <- match.arg(p.adjust.method, p.adjust.methods)
res$adj.pvalue <- p.adjust(res$pvalue, method = p.adjust.method)
res$adj.method <- p.adjust.method
res$sub.pvalue <- subpvalue
res$subni.pvalue <- subpvalue*res$pvalue
res$call <- call
class(res) <- "subnikrandtest"
return(res)
}
polak51/WitOMI documentation built on April 4, 2020, 6:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions subniche subparam.refor rtestrefor subparam.subor rtestsubor
Documented in rtestrefor rtestsubor subniche subparam.refor subparam.subor
#' @import ade4
#' @title The Within Outlying Mean Indexes calculation
#' @aliases subniche print.subniche summary.subniche refparam rtest.subniche subparam.refor rtestrefor subparam.subor rtestsubor
#' @description The indexes allows to divide the niche, estimated from the \link[ade4]{niche} function in the \link{ade4} package into subniches defined by a factor, which creates the
#' subsets. See details for more information.
#' @usage subniche(nic, factor)
#' @param nic an object of class \code{niche}.
#' @param factor a factor which will defined the subsets within which the subniches will be calculated (the same length of the number of sites)
#' @param x an object of class \code{subniche}.
#' @param xtest an object of class \code{subniche}.
#' @param object an object of class \code{subniche}.
#' @param ... further arguments passed to or from other methods
#' @param sim a numeric vector of simulated values
#' @param obs a numeric vector of an observed value
#' @param alter a character string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two-sided".he length must be equal to the length of the vector obs, values are recycled if shorter.
#' @param names a vector of names for tests
#' @param subpvalue the subset pvalue resulting from \code{subkrandtest} function
#' @param p.adjust.method a string indicating a method for multiple adjustment, see \link[stats]{p.adjust.methods} for possible choices.
#' @param call a call order
#' @param nrepet the number of permutations for the testing procedure
#' @return Adds items in the niche list and changing the class into \code{subniche} containing:
#' @return \code{factor} the factor use to divide the environmental and species matrix into submatrices.
#' @return \code{G_k} a dataframe with the sub-origins, \emph{G_k}.
#' @return \code{sub} a dataframe with the species subniche coordinates
#' @keywords subniche
#' @author Stephane Karasiewicz, \email{stephane.karasiewicz@wanadoo.fr}
#' @details The Within Outlying Mean Index analysis is a statistical exploratory niche analysis which provides observation of niche shift and/or conservatism, of an entire community,at different subcales
#' (temporal ,spatial and/or finer biological organisation level), and comparable under the same environmental gradients. This hindcasting multivariate analysis is based on the OMI analysis (Doledec \emph{et al.} 2000) which is used as reference.
#' The niches refinement is inspired by the K-select (Calenge \emph{et al.} 2005) which emphasizes the limiting factors in habitat use in design II and III (Thomas and Taylor, 1990).The different estimations should help understand:
#'
#' 1. the environmental factors defining a species' reference niche, under on the full scale, within a community.
#'
#' 2. the environmental factors defining a species' subniches, under each subsets, within a community.
#'
#' The subniches parameters can be calculated from both the reference origin,\emph{G}, which corresponds to the reference plan origin, and from \emph{G_k}, which corresponds to the
#' suborigins. \emph{G} is the graphical representation of the mean environmental conditions encountered over the full scale of the data. \emph{G_k} is the mean environmental conditions
#' encountered at a subset defined by the factor. They are complementary has you can compare:
#'
#' 1. a single species' subniches to \emph{G}.
#'
#' 2. the community' subniches to \emph{G_k} at a specific subset.
#'
#' The subniches of a single species can only be compared to \emph{G} as it is the common origin to all subsets. Whereas \emph{G_k} is only common to the species found within the
#' subset. So comparing different subniches of one species, found within different subsets, is only relevant to \emph{G}. The community's subniches can be compared to both \emph{G}
#' and \emph{G_k}, but \emph{G}, being the mean environmental conditions found within the full scale, will not express the specificity of the environmental conditions that the species
#' encountered at the subset. \emph{G_k}, being the mean environmental conditions of the subset, will reflect the atypical value of the environmental condition, making the
#' comparison of the community's subniches parameters more relevant. More information on the ecological concept can be found in Karasiewicz \emph{et al.} 2017.
#'
#' For more details description on the package use:\url{https://github.com/KarasiewiczStephane/WitOMI}.
#'
#' @references Karasiewicz S.,Doledec S.and Lefebvre S. (2017). Within outlying mean indexes: refining the OMI analysis for the realized niche decomposition. \emph{PeerJ} 5:e3364. \url{https://doi.org/10.7717/peerj.3364}.
#'
#' Calenge C., Dufour A.B. and Maillard D. (2005). K-select analysis: a new method to analyse habitat selection in radio-tracking studies. \emph{Ecological modelling}, \bold{186}, 143-153.
#'
#' Doledec S., Chessel D. and Gimaret C. (2000). Niche separation in community analysis: a new method. \emph{Ecology},\bold{81}, 2914-1927.
#'
#' Thomas, D.L., Taylor, E.J. (1990). Study Designs and Tests for Comparing Resource Use and Availability II. \emph{Natl. Widl.} \bold{54}(2), 322-330.
#'
#' @seealso \link[ade4]{niche} \link[ade4]{niche.param}
#' @examples
#'library(subniche)
#'data(doubs)
#'dudi1 <- dudi.pca(doubs$env, scale = TRUE, scan = FALSE, nf = 3)
#'nic1 <- niche(dudi1, doubs$fish, scann = FALSE)
#'# number of sites
#'N <- dim(nic1$ls)[1]
#'#Create a factor which defines the subsets
#'fact <- factor(c(rep(1,N/2),rep(2,N/2)))
#'# nic1 will be use as reference and fact will be use to define the subniches environment
#'subnic1 <- subniche(nic1, fact)
#'# the following two functions do the same display, plot.refniche is adapted to subniche objects
#'plot(nic1)
#'plot(subnic1)
#'#Display the marginality vector of the suborigins and the species subniche
#'margvect(subnic1)
#'#Display the subset's polygon, found within the overall environment's chull,
#'#and the corresponding species positions
#'subplot(subnic1)
#'# The following two functions do the same display, refparam is adapted to subniche objects
#'niche.param(nic1)
#'refparam(subnic1)
#'# The following two functions do the same display, rtest is adapted to subniche objects
#'rtest(nic1,10)
#'rtest(subnic1,10)
#'#Calculates the subniches' parameters from G with the corresponding rtest
#'subparam.refor(subnic1)
#'rtestrefor(subnic1,10)
#'#Calculates the subniches' parameters from G_k with the corresponding rtest
#'subparam.subor(subnic1)
#'rtestsubor(subnic1,10)
#' @export subniche
#' @rdname subniche
#' @importFrom stats weighted.mean
subniche <- function(nic, factor){
if (!inherits(nic, "niche"))
stop("Object of class niche expected")
appel <- as.list(nic$call)
Y <- eval.parent(appel[[3]])
liani <- split(as.data.frame(nic$ls), factor)
liwei <- sep.factor.row(Y,factor)
for(i in 1:length(liwei)){
w <- liwei[[i]]
w2 <- apply(w, 2, sum)
liwei[[i]] <- sweep(w, 2, w2, "/")
}
G_k <- as.data.frame(t(as.matrix(data.frame(lapply(liani,
function(x) apply(x, 2, mean))))))
names(G_k) <- names(nic$li)
mutemp <- list()
for (i in 1:length(liwei)){
mutemp[[i]] <- matrix(nrow=length(colnames(Y)), ncol=nic$nf)
rownames(mutemp[[i]]) <- paste(colnames(Y),levels(factor)[[i]],sep="")
for(j in 1:length(colnames(Y))){
z <- liwei[[i]][,j]
mutemp[[i]][j,] <- apply(liani[[i]],2, weighted.mean, z)
}
}
mutemp <- do.call("rbind", mutemp)
nic$sub <- mutemp
nic$G_k <- G_k
rownames(nic$G_k) <- paste(rep("G_k",dim(nic$G_k)[1]),levels(factor),sep="")
nic$factor <- factor
class(nic) <- c("subniche", "dudi")
return(nic)
}
#' @rdname subniche
#' @method print subkrandtest
#' @export
print.subkrandtest <- function (x, ...)
{
if (!inherits(x, "subkrandtest"))
stop("to be used with 'subkrandtest' object")
cat("class:", class(x), "\n")
cat("Monte-Carlo tests\n")
cat("Call: ")
print(x$call)
cat("\nNumber of tests: ", x$ntest, "\n")
cat("\nAdjustment method for multiple comparisons: ", x$adj.method,
"\n")
sumry <- list(Test = x$names, Obs = x$obs, Std.Obs = x$expvar[,
1], Alter = x$alter)
sumry <- as.data.frame(sumry)
row.names(sumry) <- 1:x$ntest
if (any(x$rep[1] != x$rep)) {
sumry <- cbind(sumry[, 1:4], N.perm = x$rep)
}
else {
cat("Permutation number: ", x$rep[1], "\n")
}
sumry <- cbind(sumry, Pvalue = x$pvalue)
if (x$adj.method != "none")
sumry <- cbind(sumry, Pvalue.adj = x$adj.pvalue)
print(sumry)
cat("Subsets Pvalue:", x$subpvalue, "\n")
if (length(names(x)) > 10) {
cat("other elements: ")
cat(names(x)[11:(length(x))], "\n")
}
}
#' @rdname subniche
#' @method print subnikrandtest
#' @export
print.subnikrandtest <- function (x, ...)
{
if (!inherits(x, "subnikrandtest"))
stop("to be used with 'subnikrandtest' object")
cat("class:", class(x), "\n")
cat("Monte-Carlo tests\n")
cat("Call: ")
print(x$call)
cat("\nNumber of tests: ", x$ntest, "\n")
cat("\nAdjustment method for multiple comparisons: ", x$adj.method,
"\n")
sumry <- list(Test = x$names, Obs = x$obs, Std.Obs = x$expvar[,
1], Alter = x$alter)
sumry <- as.data.frame(sumry)
row.names(sumry) <- 1:x$ntest
if (any(x$rep[1] != x$rep)) {
sumry <- cbind(sumry[, 1:4], N.perm = x$rep)
}
else {
cat("Permutation number: ", x$rep[1], "\n")
}
sumry <- cbind(sumry, Pvalue = x$pvalue, SubniPvalue = x$subni.pvalue)
if (x$adj.method != "none")
sumry <- cbind(sumry, Pvalue.adj = x$adj.pvalue, SubniPvalue = x$subni.pvalue)
print(sumry)
cat("\n")
if (length(names(x)) > 9) {
cat("other elements: ")
cat(names(x)[10:(length(x))], "\n")
}
}
#' @rdname subniche
#' @method print subniche
#' @export
print.subniche <- function (x, ...)
{
if (!inherits(x, "subniche"))
stop("to be used with 'niche' object")
cat("WitOMI calculation\n")
cat("call: ")
print(x$call)
cat("class: ")
cat(class(x), "\n")
cat("\n$rank (rank) :", x$rank)
cat("\n$nf (axis saved) :", x$nf)
cat("\n$RV (RV coeff) :", x$RV)
cat("\n\neigen values: ")
l0 <- length(x$eig)
cat(signif(x$eig, 4)[1:(min(5, l0))])
if (l0 > 5)
cat(" ...\n\n")
else cat("\n\n")
sumry <- array("", c(4, 4), list(1:4, c("vector", "length",
"mode", "content")))
sumry[1, ] <- c("$eig", length(x$eig), mode(x$eig), "eigen values")
sumry[2, ] <- c("$lw", length(x$lw), mode(x$lw), "row weigths (crossed array)")
sumry[3, ] <- c("$cw", length(x$cw), mode(x$cw), "col weigths (crossed array)")
sumry[4, ] <- c("$factor", length(x$factor), mode(x$factor), "factor used for creating subsets")
print(sumry, quote = FALSE)
cat("\n")
sumry <- array("", c(9, 4), list(1:9, c("data.frame", "nrow",
"ncol", "content")))
sumry[1, ] <- c("$tab", nrow(x$tab), ncol(x$tab), "crossed array (averaging species/sites)")
sumry[2, ] <- c("$li", nrow(x$li), ncol(x$li), "species coordinates")
sumry[3, ] <- c("$l1", nrow(x$l1), ncol(x$l1), "species normed scores")
sumry[4, ] <- c("$co", nrow(x$co), ncol(x$co), "variables coordinates")
sumry[5, ] <- c("$c1", nrow(x$c1), ncol(x$c1), "variables normed scores")
sumry[6, ] <- c("$ls", nrow(x$ls), ncol(x$ls), "sites coordinates")
sumry[7, ] <- c("$as", nrow(x$as), ncol(x$as), "axis upon niche axis")
sumry[8, ] <- c("$G_k", nrow(x$G_k), ncol(x$G_k), "G_k coordinates")
sumry[9, ] <- c("$sub", nrow(x$sub), ncol(x$sub), "species coordinates within each subset")
print(sumry, quote = FALSE)
cat("\n")
}
#' @rdname subniche
#' @method summary subniche
#' @export
summary.subniche <- function (object, ...)
{
cat("Class: ")
cat(class(object))
cat("\nCall: ")
print(object$call)
cat("\nTotal inertia: ")
cat(signif(sum(object$eig), 4))
cat("\n")
l0 <- length(object$eig)
cat("\nEigenvalues:\n")
vec <- object$eig[1:(min(5, l0))]
names(vec) <- paste("Ax", 1:length(vec), sep = "")
print(format(vec, digits = 4, trim = TRUE, width = 7), quote = FALSE)
cat("\nProjected inertia (%):\n")
vec <- (object$eig/sum(object$eig) * 100)[1:(min(5, l0))]
names(vec) <- paste("Ax", 1:length(vec), sep = "")
print(format(vec, digits = 4, trim = TRUE, width = 7), quote = FALSE)
cat("\nCumulative projected inertia (%):\n")
vec <- (cumsum(object$eig)/sum(object$eig) * 100)[1:(min(5,
l0))]
names(vec)[1] <- "Ax1"
if (l0 > 1)
names(vec)[2:length(vec)] <- paste("Ax1:", 2:length(vec),
sep = "")
print(format(vec, digits = 4, trim = TRUE, width = 7), quote = FALSE)
if (l0 > 5) {
cat("\n")
cat(paste("(Only 5 dimensions (out of ", l0, ") are shown)\n",
sep = "", collapse = ""))
}
cat("\n")
}
#' @rdname subniche
#' @export
refparam <- function (x)
{
if (!inherits(x, "subniche"))
stop("Object of class 'subniche' expected")
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab
Y <- eval.parent(appel[[3]])
w1 <- apply(Y, 2, sum)
if (any(w1 <= 0))
stop(paste("Column sum <=0 in Y"))
Y <- sweep(Y, 2, w1, "/")
calcul.param <- function(freq, mil) {
inertia <- sum(freq * mil * mil)
m <- apply(freq * mil, 2, sum)
margi <- sum(m^2)
mil <- t(t(mil) - m)
tolt <- sum(freq * mil * mil)
u <- m/sqrt(sum(m^2))
z <- mil %*% u
tolm <- sum(freq * z * z)
tolr <- tolt - tolm
w <- c(inertia, margi, tolm, tolr)
names(w) <- c("inertia", "OMI", "Tol", "Rtol")
w1 <- round(w[2:4]/w[1], digits = 3) * 100
names(w1) <- c("omi", "tol", "rtol")
return(c(w, w1))
}
res <- apply(Y, 2, calcul.param, mil = X)
t(res)
}
#' @rdname subniche
#' @method rtest subniche
#' @export
rtest.subniche <- function (xtest, nrepet = 99, ...)
{
if (!inherits(xtest, "dudi"))
stop("Object of class dudi expected")
if (!inherits(xtest, "subniche"))
stop("Type 'niche' expected")
appel <- as.list(xtest$call)
X <- eval.parent(appel$dudiX)$tab
Y <- eval.parent(appel$Y)
w1 <- apply(Y, 2, sum)
if (any(w1 <= 0))
stop(paste("Column sum <=0 in Y"))
Y <- sweep(Y, 2, w1, "/")
calcul.margi <- function(freq, mil) {
m <- apply(freq * mil, 2, sum)
return(sum(m^2))
}
obs <- apply(Y, 2, calcul.margi, mil = X)
obs <- c(obs, OMI.mean = mean(obs))
sim <- sapply(1:nrepet, function(x) apply(apply(Y, 2, sample),
2, calcul.margi, mil = X))
sim <- rbind(sim, OMI.mean = apply(sim, 2, mean))
res <- as.krandtest(obs = obs, sim = t(sim))
return(res)
}
#' @rdname subniche
#' @export
subparam.refor <- function(x){
res <- list()
appel <- as.list(x$call)
y <-eval.parent(appel[[3]])
factor <- x$factor
subniche.param <- function (x,y)
{
if (!inherits(x, "subniche"))
stop("Object of class 'subniche' expected")
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab[y,]
Y <- eval.parent(appel[[3]])[y,]
w1 <- apply(Y, 2, sum)
Y <- sweep(Y, 2, w1, "/")
calcul.param <- function(freq, mil) {
inertia <- sum(freq * mil * mil)
m <- apply(freq * mil, 2, sum)
margi <- sum(m^2)
mil <- t(t(mil) - m)
tolt <- sum(freq * mil * mil)
u <- m/sqrt(sum(m^2))
z <- mil %*% u
tolm <- sum(freq * z * z)
tolr <- tolt - tolm
w <- c(inertia, margi, tolm, tolr)
names(w) <- c("inertia", "WitOMIG", "Tol", "Rtol")
w1 <- round(w[2:4]/w[1], digits = 3) * 100
names(w1) <- c("witomig", "tol", "rtol")
return(c(w, w1))
}
res <- apply(Y, 2, calcul.param, mil = X)
t(res)
}
N <-length(levels(x$factor))
nam <- levels(factor)
for(i in 1:N){
res[[i]] <- subniche.param(x,array(factor)==nam[i])
rownames(res[[i]]) <- paste(rownames(x$li),nam[i],sep="")
}
res <- do.call("rbind",res)
return(res)
}
#' @rdname subniche
#' @export
rtestrefor <- function(x, nrepet){
factor <- x$factor
res <- list()
N <-length(levels(x$factor))
nam <- levels(factor)
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab
Y <- eval.parent(appel[[3]])
calcul.margi <- function(freq, mil) {
m <- apply(freq * mil, 2, sum)
return(sum(m^2))
}
for(i in 1:N){
X1 <- X[array(factor)==nam[i],]
Xwobs <- apply(X1, 2, mean)
Xsim <- sapply(1:nrepet, function(x) apply(t <- apply(X, 2, sample)[array(factor)==nam[i],], 2, mean))
Xtest <- subkrandtest(obs=Xwobs, sim=t(Xsim), "two-sided")
w1 <- apply(Y[array(factor)==nam[i],], 2, sum)
Y1 <- sweep(Y[array(factor)==nam[i],], 2, w1, "/")
obs <- apply(Y1, 2, calcul.margi, mil = X1)
obs <- c(obs, OMI.mean = mean(obs,na.rm=TRUE))
sim <- sapply(1:nrepet, function(x) apply(sweep(t <- apply(Y, 2, sample)[array(factor)==nam[i],], 2,apply(t, 2, sum),"/"), 2, calcul.margi, mil = X1))
sim <- rbind(sim, OMI.mean = apply(sim, 2, mean,na.rm=TRUE))
omitest <- subnikrandtest(obs = obs, sim = t(sim),subpvalue= Xtest$subpvalue)
res[[i]] <- list("Subsettest" = Xtest, "witomigtest" = omitest)
}
names(res) <- nam
return(res)
}
#' @rdname subniche
#' @export
subparam.subor <- function(x){
res <- list()
appel <- as.list(x$call)
y <-eval.parent(appel[[3]])
factor <- x$factor
subnichesub.param <- function (x,y)
{
if (!inherits(x, "subniche"))
stop("Object of class 'subniche' expected")
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab[y,]
X <- as.matrix(scale(X, center=T,scale=F))
Y <- eval.parent(appel[[3]])[y,]
w1 <- apply(Y, 2, sum)
Y <- sweep(Y, 2, w1, "/")
calcul.param <- function(freq, mil) {
inertia <- sum(freq * mil * mil)
m <- apply(freq * mil, 2, sum)
margi <- sum(m^2)
mil <- t(t(mil) - m)
tolt <- sum(freq * mil * mil)
u <- m/sqrt(sum(m^2))
z <- mil %*% u
tolm <- sum(freq * z * z)
tolr <- tolt - tolm
w <- c(inertia, margi, tolm, tolr)
names(w) <- c("inertia", "WitOMIG_k", "Tol", "Rtol")
w1 <- round(w[2:4]/w[1], digits = 3) * 100
names(w1) <- c("witomig_k", "tol", "rtol")
return(c(w, w1))
}
res <- apply(Y, 2, calcul.param, mil = X)
t(res)
}
N <-length(levels(x$factor))
nam <- levels(factor)
for(i in 1:N){
res[[i]] <- subnichesub.param(x,array(factor)==nam[i])
rownames(res[[i]]) <- paste(rownames(x$li),i,sep="")
}
res <- do.call("rbind",res)
return(res)
}
#' @rdname subniche
#' @export
rtestsubor <- function(x, nrepet){
factor <- x$factor
res <- list()
N <-length(levels(x$factor))
nam <- levels(factor)
appel <- as.list(x$call)
X <- eval.parent(appel[[2]])$tab
Y <- eval.parent(appel[[3]])
calcul.margi <- function(freq, mil) {
m <- apply(freq * mil, 2, sum)
return(sum(m^2))
}
for(i in 1:N){
X1 <- X[array(factor)==nam[i],]
Xwobs <- apply(X1, 2, mean)
Xobs <- sweep(X1, 2,Xwobs,"-")
Xsim <- sapply(1:nrepet, function(x) apply(t <- apply(X, 2, sample)[array(factor)==nam[i],], 2, mean))
Xtest <- subkrandtest(obs=Xwobs, sim=t(Xsim), "two-sided")
w1 <- apply(Y[array(factor)==nam[i],], 2, sum)
Y1 <- sweep(Y[array(factor)==nam[i],], 2, w1, "/")
obs <- apply(Y1, 2, calcul.margi, mil = Xobs)
obs <- c(obs, OMI.mean = mean(obs,na.rm=TRUE))
sim <- sapply(1:nrepet, function(x) apply(sweep(t <- apply(Y, 2, sample)[array(factor)==nam[i],], 2,apply(t, 2, sum),"/"), 2, calcul.margi, mil = Xobs))
sim <- rbind(sim, OMI.mean = apply(sim, 2, mean,na.rm=TRUE))
omitest <- subnikrandtest(obs = obs, sim = t(sim),subpvalue=Xtest$subpvalue)
res[[i]] <- list("Subsettest" = Xtest, "witomig_ktest" = omitest)
}
names(res) <- nam
return(res)
}
#' @rdname subniche
#' @export
#' @importFrom stats sd na.omit p.adjust p.adjust.methods
subkrandtest<-function (sim, obs, alter = "greater", call = match.call(), names = colnames(sim),
p.adjust.method = "none")
{
res <- list(sim = sim, obs = obs)
if (length(obs) != length(alter))
alter <- rep(alter, length = length(obs))
res$alter <- alter
res$rep <- apply(sim, 2, function(x) length(na.omit(x)))
res$ntest <- length(obs)
res$expvar <- data.frame(matrix(0, res$ntest, 3))
if (!is.null(names)) {
res$names <- names
}
else {
res$names <- paste("test", 1:res$ntest, sep = "")
}
names(res$expvar) <- c("Std.Obs", "Expectation", "Variance")
res$pvalue <- rep(0, length(obs))
for (i in 1:length(obs)) {
vec.sim <- na.omit(sim[, i])
res$alter[i] <- match.arg(res$alter[i], c("greater",
"less", "two-sided"))
res$expvar[i, 1] <- (res$obs[i] - mean(vec.sim))/sd(vec.sim)
res$expvar[i, 2] <- mean(vec.sim)
res$expvar[i, 3] <- sd(vec.sim)
if (res$alter[i] == "greater") {
res$pvalue[i] <- (sum(vec.sim >= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "less") {
res$pvalue[i] <- (sum(vec.sim <= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "two-sided") {
sim0 <- abs(vec.sim - mean(vec.sim))
obs0 <- abs(obs[i] - mean(vec.sim))
res$pvalue[i] <- (sum(sim0 >= obs0) + 1)/(res$rep[i] +
1)
}
}
p.adjust.method <- match.arg(p.adjust.method, p.adjust.methods)
res$adj.pvalue <- p.adjust(res$pvalue, method = p.adjust.method)
res$adj.method <- p.adjust.method
res$subpvalue <- prod(res$pvalue)
res$call <- call
class(res) <- "subkrandtest"
return(res)
}
#' @rdname subniche
#' @export
#' @importFrom stats sd na.omit p.adjust p.adjust.methods
subnikrandtest<-function (sim, obs, alter = "greater",subpvalue, call = match.call(), names = colnames(sim),
p.adjust.method = "none")
{
res <- list(sim = sim, obs = obs)
if (length(obs) != length(alter))
alter <- rep(alter, length = length(obs))
res$alter <- alter
res$rep <- apply(sim, 2, function(x) length(na.omit(x)))
res$ntest <- length(obs)
res$expvar <- data.frame(matrix(0, res$ntest, 3))
if (!is.null(names)) {
res$names <- names
}
else {
res$names <- paste("test", 1:res$ntest, sep = "")
}
names(res$expvar) <- c("Std.Obs", "Expectation", "Variance")
res$pvalue <- rep(0, length(obs))
for (i in 1:length(obs)) {
vec.sim <- na.omit(sim[, i])
res$alter[i] <- match.arg(res$alter[i], c("greater",
"less", "two-sided"))
res$expvar[i, 1] <- (res$obs[i] - mean(vec.sim))/sd(vec.sim)
res$expvar[i, 2] <- mean(vec.sim)
res$expvar[i, 3] <- sd(vec.sim)
if (res$alter[i] == "greater") {
res$pvalue[i] <- (sum(vec.sim >= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "less") {
res$pvalue[i] <- (sum(vec.sim <= obs[i]) + 1)/(res$rep[i] +
1)
}
else if (res$alter[i] == "two-sided") {
sim0 <- abs(vec.sim - mean(vec.sim))
obs0 <- abs(obs[i] - mean(vec.sim))
res$pvalue[i] <- (sum(sim0 >= obs0) + 1)/(res$rep[i] +
1)
}
}
p.adjust.method <- match.arg(p.adjust.method, p.adjust.methods)
res$adj.pvalue <- p.adjust(res$pvalue, method = p.adjust.method)
res$adj.method <- p.adjust.method
res$sub.pvalue <- subpvalue
res$subni.pvalue <- subpvalue*res$pvalue
res$call <- call
class(res) <- "subnikrandtest"
return(res)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.