MicEco：计算Sloan随机性的另一方法

1library(devtools)
2install_github("Russel88/MicEco")

 1#' Fit Sloan et al. (2006) Neutral Model
 2#'
 3#' Fit neutral model developed by Sloan et al. (2006, Environ Microbiol 8(4):732-740) and implemented by Burns et al. (2015, ISME J 10(3):655-664).
 4#' @param otu An OTU-table with taxa as columns and samples as rows.
 5#' @keywords neutral
 6#' @return A list of length two; first element contains fit statistics, the second element contains predictions.
 7#' @import bbmle
 8#' @export
 9
10neutral.fit <- function(otu){
11
12  # Frequency
13  otu.pa <- (otu>0)*1
14  freq <- apply(otu.pa, 2, mean)
15
16  # Individuals per community
17  N <- mean(apply(otu, 1, sum))
18
19  # Relative abundance
20  p <- apply(otu, 2, function(x) mean(x))/N
21
22  # Detection limit
23  d = 1/N
24
25  # Define likelihood function
26  neutral.ll <- function(m, sigma){
27    R = freq - pbeta(d, N*m*p, N*m*(1-p), lower.tail=FALSE)
28    -sum(dnorm(R, 0, sigma,log=TRUE))
29  }
30
31  # Fit neutral model
32  m.mle <- mle2(neutral.ll, start=list(m=0.01, sigma=0.1),method="Nelder-Mead")
33
34  # R-squared
35  gRsqr <- 1 - exp(-as.numeric(logLik(m.mle))/length(p))
36
37  # Predictions
38  freq.pred <- pbeta(d, N*m.mle@coef['m']*p, N*m.mle@coef['m']*(1-p), lower.tail=FALSE)
39  pred.ci <- binconf(freq.pred*nrow(otu), nrow(otu), return.df=TRUE)
40
41  # Bind results
42  stats <- cbind(m.mle@coef['m'], m.mle@details$value,gRsqr , N, nrow(otu), length(p), d)
43  pred <- cbind(p, freq, freq.pred, pred.ci[,2:3])
44
45  results <- list(stats,pred)
46  return(results)
47}