R/mom_nosource.R

In kathalexknuts/NeutralTheorySimulations: Simulations for Neutral Theory Analysis

#' Method of Moments Estimating of Migration Paramter
#'
#' This function estimates the migration parameter, given only data on a local community.
#' This method is an extension from the paper "Taxa-Abundance Distributions in Microbial Communities Using Environmental
#' Sequence Data" by Sloan et al.
#' @param local A local community of OTU reads for n samples. Rows should be OTUs, Columns should be Samples. A column of taxon indexes should be included at the beginning. Note that if this is not included, the first sample in your data will not be used.
#' @param pool_est to pool meta/source community data during estimation
#' @import vegan 
#' @return a list of 2 elements. Element 1 is the estimate for the migration parameter. Element 2 is the mean bray curtis dissimilarity in the local community sample.
#' @export

mom_nosource <- function(local, meta, pool_est = FALSE)
{
  
  BC_local <- vegdist(t(local[,-1]), method = "bray", diag = FALSE, upper = FALSE, na.rm = TRUE)
  BC_meta <- vegdist(t(meta[,-1]), method = "bray", diag = FALSE, upper = FALSE, na.rm = TRUE)
  #meanBC <- mean(BC)
  #BC_dim <- sum((c(BC) - meanBC)^2)
  BC_dim_local <- mean(BC_local)
  BC_dim_meta <- mean(BC_meta)
  
  lungTot=apply(local[,2:ncol(local)],2,sum) #sum of Nlocal for each sample
  NT=mean(lungTot) 
  d=1/NT
  
  if(pool_est == TRUE){
    oralMn <- rowSums(local[,2:ncol(local)])/sum(rowSums(local[,2:ncol(local)]))
  }
  else{
    oralpik=t(t(local[,2:ncol(local)])/apply(local[,2:ncol(local)],2,sum)) #relative abundances for each OTU within each sample, OTUs in rows
    oralMn=apply(oralpik,1,mean) #average relative abundance for each taxa taken across all samples
  }
  
  lungFr=apply(local[,2:ncol(local)],1,f1 <- function(x){sum(x>0)})/(ncol(local)-1) #proportion of samples with non-zero reads for each OTU
  
  minFcn <- function(m){
    sum((lungFr-1+pbeta(d,NT*m*oralMn, NT*m*(1-oralMn)))^2) #mimimizing diff between "observed" proportion of detection in local community & expected based on source community under Neutral Theory assumption 
  } #as the beta parameter gets large, peak shrinks and shifts left - aka when p is small
  
  f1Oral=optimize(f=minFcn,lower=1.0e-7,upper=1-1.0e-7)
  
  return(list(f1Oral$min, BC_dim_local, BC_dim_meta))
}