fitmodel <- function(x, model, count, nRand = 999, ...){
if(!(model %in% c("dominanceDecay", "dominancePreemp", "MacArthurFraction",
"powerFraction", "randAssort", "randFraction"))){
stop(paste("\n'model' must be one of this:", '"dominanceDecay", "dominancePreemp", "MacArthurFraction", "powerFraction", "randAssort", "randFraction"'))
}
if(!is.matrix(x) & !is.data.frame(x)){
stop("\n'x' must be a data frame or a matrix")
}
if(length(x[ ,1]) < 2){
stop("\nYou must have more than one replicates in lines")
}
if(is.null(count)){
stop("\n'count' must be specified")
}
# Arguments in dots
dots <- list(...)
# Number of replicates and ranks.
n <- dim(x)[1]
Rk <- dim(x)[2]
# Total abundance and species of each replicates.
N <- apply(x, 1, sum)
S <- apply(x > 0, 1, sum)
# 'nRand' means and variances of 'n' simulations to the model.
M <- matrix(nrow = nRand, ncol = Rk)
V <- matrix(nrow = nRand, ncol = Rk)
for(i in 1:nRand){
sim <- matrix(0, nrow = n, ncol = Rk)
for(j in 1:n){
sim[j,1:S[j]] <- do.call(model, c(list(N = N[j], S = S[j], count = count), dots))
# Transform to relative abundance
sim[j,1:S[j]] <- sim[j,1:S[j]] / sum(sim[j,1:S[j]])
}
M[i, ] <- apply(sim, 2, mean)
V[i, ] <- apply(sim, 2, var)
}
# Transform each replicate to relative abundance.
for(i in 1:n){
x[i, ] <- sort(x[i, ] / sum(x[i, ]), decreasing = TRUE)
}
# Observed relative abundance mean and variance of replicates.
M0 <- apply(x, 2, mean)
V0 <- apply(x, 2, var)
# Probability that the observed mean and variance are predicted by the model.
pM0 <- c()
pV0 <- c()
for(i in 1:Rk){
# p = (b+1)/(m+1)
pM0[i] <- 2 * min((sum(M[ ,i] < M0[i]) + 1) / (nRand + 1),
(sum(M[ ,i] > M0[i]) + 1) / (nRand + 1))
pV0[i] <- 2 * min((sum(V[ ,i] < V0[i]) + 1) / (nRand + 1),
(sum(V[ ,i] > V0[i]) + 1) / (nRand + 1))
}
# Global statistic (T observed) to combined p of all ranks.
TM0 <- -2 * sum(log(pM0))
TV0 <- -2 * sum(log(pV0))
# Distribution of T values.
dTM <- c()
dTV <- c()
for(i in 1:nRand){
pM <- c()
pV <- c()
for(j in 1:Rk){
# p = (b+1)/(m+1)
pM[j] <- 2 * min(((sum(c(M[-i,j], M0[j]) < M[i,j]) + 1) / (nRand + 1)),
((sum(c(M[-i,j], M0[j]) > M[i,j]) + 1) / (nRand + 1)))
pV[j] <- 2 * min(((sum(c(V[-i,j], V0[j]) < V[i,j]) + 1) / (nRand + 1)),
((sum(c(V[-i,j], V0[j]) > V[i,j]) + 1) / (nRand + 1)))
}
dTM[i] <- -2 * sum(log(pM))
dTV[i] <- -2 * sum(log(pV))
}
# Probability that T observed is drawn from T values generated by the model.
pvalueM <- sum(dTM > TM0) / (nRand + 1)
pvalueV <- sum(dTV > TV0) / (nRand + 1)
# Simulation range for mean and variance.
rM <- matrix(c(apply(M, 2, min), apply(M, 2, max)), nrow = 2, ncol = Rk, byrow = TRUE,
dimnames = list(c("min", "max"), paste("rank", 1:Rk, sep = "")))
rV <- matrix(c(apply(V, 2, min), apply(V, 2, max)), nrow = 2, ncol = Rk, byrow = TRUE,
dimnames = list(c("min", "max"), paste("rank", 1:Rk, sep = "")))
return(new("fitmodel",
call = list(model = model, nRepl = n, nRank = Rk, nRand = nRand,
count = count),
Tstats = list(dTmean = dTM, dTvar = dTV, TMobs = TM0, TVobs = TV0,
pvalue = matrix(c(pvalueM, pvalueV), nrow = 2, ncol = 1,
dimnames = list(c("mean", "variance"),
"p-value"))),
sim.stats = matrix(c(apply(M, 2, mean), apply(V, 2, mean)), nrow = 2,
ncol = Rk, byrow = TRUE,
dimnames = list(c("mean", "variance"),
paste("rank", 1:Rk, sep = ""))),
sim.range = list(mean = rM, variance = rV),
obs.stats = matrix(c(M0, V0), nrow = 2, ncol = Rk, byrow = TRUE,
dimnames = list(c("mean", "variance"),
paste("rank", 1:Rk, sep = "")))))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.