Merge branch 'dev' to master

Release 1.1.0

DESCRIPTION

@@ -6,8 +6,12 @@ Maintainer: Mario J. Marques-Azevedo <mariojmaaz@gmail.com>
 Description: Simulate and fit, using Monte Carlo approach, rank abundance distribution to niche apportionment species abundance distributions models.
 License: GPL-2
 Encoding: UTF-8
-Depends: parallel
+Depends: parallel, methods
+Imports: graphics
+Suggests: knitr
+VignetteBuilder: knitr
+BuildVignettes: yes
 NeedsCompilation: no
-Version: 1.0.2
-Date: 2015-12-04
+Version: 1.1.0
+Date: 2015-12-07
 URL: http://mariojose.github.io/nicheApport/, https://github.com/mariojose/nicheApport

NAMESPACE

@@ -3,3 +3,18 @@ export(dominanceDecay, dominancePreemp, MacArthurFraction, powerFraction, randAs
 
 # Export fitting function
 export(fitmodel, fitmodelCl, findPFw)
+
+# Export Classes
+exportClasses(fitmodel)
+
+# Export methods
+exportMethods(plot, lines, qqplot, hist, show)
+
+# Import function from packages
+importFrom("graphics", plot, lines, hist)
+
+# Import packages
+import(
+    parallel,
+    methods
+)

NEWS

@@ -1,3 +1,15 @@
+1.1.0 (2015-12-14)
+* Improve runtime of functions
+* Joined documentation of randFraction and MacArthurFraction to powerFraction
+* Bug fix to dominancePreemp that return ranks with 0 when using continuous
+* Add news arguments to control plots colours and line type
+* Bux fix to plot and fit discrete data
+* Add vignette 'Quick Reference for nicheApport'
+* Update documentation
+* Improve fitmode and fitmodelCl functions
+* Add methods to plot rank abundance, simulated model line, qqplot of observed and simulated data and histogram of distribution of T value
+* Now fitmodel and fitmodelCl return object of class 'fitmodel'
+
 1.0.2 (2015-12-04)
 * fitmodel and fitmodelCl now return simulation range
 * Update interpretation of p-value in help files

R/MacArthurFraction.R

@@ -1,11 +1,4 @@
 # MacArthur Fraction
 MacArthurFraction <- function (N, S, count = FALSE){
-  k <- runif(n = S, min = 0, max = 1)
-  r <- N
-  for(i in 2:S){
-    j <- sample(x = 1:(i - 1), size = 1, prob = (r / sum(r)))
-    r[i] <- ifelse(count, floor(r[j] * k[i]), r[j] * k[i])
-    r[j] <- abs(r[j] - r[i])
-  }
-  return(sort(r, decreasing = TRUE))
+  return(powerFraction(N = N, S = S, w = 1, count = count))
 }

R/dominanceDecay.R

@@ -1,13 +1,23 @@
 # Dominance Decay
 dominanceDecay <- function (N, S, count = FALSE){
   k <- runif(n = S, min = 0, max = 1)
-  r <- N
-  r[1] <- ifelse(count, floor(r[1] * k[1]), r[1] * k[1])
-  r[2] <- N - r[1]
-  for(i in 3:S){
-    r <- sort(r, decreasing = FALSE)
-    r[i] <- ifelse(count, floor(r[i - 1] * k[i]), r[i - 1] * k[i])
-    r[i - 1] <- abs(r[i - 1] - r[i])
+  if(count){
+    r <- floor(N * k[1])
+    r[2] <- N - r[1]
+    for(i in 3:S){
+      r <- sort(r, decreasing = FALSE)
+      r[i] <- floor(r[i - 1] * k[i])
+      r[i - 1] <- abs(r[i - 1] - r[i])
+    }
+  }
+  else{
+    r <- N * k[1]
+    r[2] <- N - r[1]
+    for(i in 3:S){
+      r <- sort(r, decreasing = FALSE)
+      r[i] <- r[i - 1] * k[i]
+      r[i - 1] <- abs(r[i - 1] - r[i])
+    }
   }
   return(sort(r, decreasing = TRUE))
 }

R/dominancePreemp.R

@@ -2,9 +2,17 @@
 dominancePreemp <- function (N, S, count = FALSE){
   k <- runif(n = (S - 1), min = 0.5, max = 1)
   r <- N
-  for(i in 1:(S - 1)){
-    r[i] <- ifelse(count, floor(r[i] * k[i]), r[i] * k[i])
-    r[i + 1] <- N - sum(r[1:i])
+  if(count){
+    for(i in 1:(S-1)){
+      r[i + 1] <- floor(r[i] * (1 - k[i]))
+      r[i] <- floor(r[i] * k[i])
+    }
+  }
+  else{
+    for(i in 1:(S-1)){
+      r[i + 1] <- r[i] * (1 - k[i])
+      r[i] <- r[i] * k[i]
+    }
   }
   return(sort(r, decreasing = TRUE))
 }

R/findPFw.R

@@ -18,10 +18,10 @@ findPFw <- function(x, count, ws = 0, we = 1, f = 0.01, nRand = 99, cl = TRUE, n
     } else {
       m <- fitmodel(x, model = "powerFraction", count = count, nRand = nRand, w = out$w[i])
     }
-    out$TMobs[i] <- m$TMobs
-    out$TVobs[i] <- m$TVobs
-    out$pvalueM[i] <- m$stat[1, ]
-    out$pvalueV[i] <- m$stat[2, ]
+    out$TMobs[i] <- m@Tstats$TMobs
+    out$TVobs[i] <- m@Tstats$TVobs
+    out$pvalueM[i] <- m@Tstats$pvalue[1, ]
+    out$pvalueV[i] <- m@Tstats$pvalue[2, ]
   }
   return(out)
 }

R/fitmodel.R

@@ -16,24 +16,29 @@ fitmodel <- function(x, model, count, nRand = 999, ...){
     stop("\n'count' must be specified")
   }
 
-  # Number of replicates and total species.
+  # Arguments in dots
+  dots <- list(...)
+
+  # Number of replicates and ranks.
   n <- dim(x)[1]
-  S <- dim(x)[2]
-  # Total abundance of each replicates.
-  N <- apply(x, MARGIN = 1, FUN = sum)
+  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 = S)
-  V <- matrix(nrow = nRand, ncol = S)
+  M <- matrix(nrow = nRand, ncol = Rk)
+  V <- matrix(nrow = nRand, ncol = Rk)
   for(i in 1:nRand){
-    sim <- matrix(nrow = n, ncol = S)
+    sim <- matrix(0, nrow = n, ncol = Rk)
     for(j in 1:n){
-      sim[j, ] <- do.call(model, list(N = N[j], S = S, count = count, ... = ...))
+      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, ] <- sim[j, ] / sum(sim[j, ])
+      sim[j,1:S[j]] <- sim[j,1:S[j]] / sum(sim[j,1:S[j]])
     }
-    M[i, ] <- apply(sim, MARGIN = 2, FUN = mean)
-    V[i, ] <- apply(sim, MARGIN = 2, FUN = var)
+    M[i, ] <- apply(sim, 2, mean)
+    V[i, ] <- apply(sim, 2, var)
   }
 
   # Transform each replicate to relative abundance.
@@ -42,13 +47,13 @@ fitmodel <- function(x, model, count, nRand = 999, ...){
   }
 
   # Observed relative abundance mean and variance of replicates.
-  M0 <- apply(x, MARGIN = 2, FUN = mean)
-  V0 <- apply(x, MARGIN = 2, FUN = var)
+  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:S){
+  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))
@@ -66,7 +71,7 @@ fitmodel <- function(x, model, count, nRand = 999, ...){
   for(i in 1:nRand){
     pM <- c()
     pV <- c()
-    for(j in 1:S){
+    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)))
@@ -82,14 +87,24 @@ fitmodel <- function(x, model, count, nRand = 999, ...){
   pvalueV <- sum(dTV > TV0) / (nRand + 1)
 
   # Simulation range for mean and variance.
-  rM <- apply(M, MARGIN = 2, FUN = range)
-  rV <- apply(V, MARGIN = 2, FUN = range)
-  dimnames(rM) <- list(c("min", "max"), paste("rank", 1:S, sep = ""))
-  rownames(rV) <- list(c("min", "max"), paste("rank", 1:S, sep = ""))
+  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(list(dTmean = dTM, dTvar = dTV, TMobs = TM0, TVobs = TV0, rangeM = rM, rangeV = rV,
-              simulations = matrix(c(apply(M, 2, mean), apply(V, 2, mean)), nrow = 2, ncol = S, byrow = TRUE,
-                                   dimnames = list(c("mean", "variance"), paste("rank", 1:S, sep = ""))),
-              stat = matrix(c(pvalueM, pvalueV), nrow = 2, ncol = 1, 
-                          dimnames = list(c("mean", "variance"), "p-value"))))
+  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 = "")))))
 }