R/fitmodel.R

Defines functions fitmodel

Documented in fitmodel

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 = "")))))
}
MarioJose/nicheApport documentation built on May 7, 2019, 2:52 p.m.