R/mom_function.R

In kathalexknuts/ntsims2:

#' Method of Moments Estimating of Migration Paramter
#'
#' This function estimates the migration parameter between a local and meta community for a single individual.
#' 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. 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 meta A meta community of OTU reads. 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.
#' @import vegan 
#' @return a list of 1 element. Element 1 is the estimate for the migration parameter. 
#' @export

mom2 <- function(local, meta)
{

  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)
  BC_dim_local <- mean(BC_local)
  BC_dim_meta <- mean(BC_meta)
  
  m_ests <- c()
  for(i in 2:ncol(local)){
  
    NT=sum(local[,i])
    d=1/NT
    oralMn=t(t(meta[,i])/sum(meta[,i]))
  
    lungFr=ifelse(local[,i] == 0, 0, 1)
  
    minFcn <- function(m){
      sum((lungFr-1+pbeta(d,NT*m*oralMn, NT*m*(1-oralMn)))^2) 
    }
  
    f1Oral=optimize(f=minFcn,lower=1.0e-7,upper=1-1.0e-7)
    m_ests <- c(m_ests, f1Oral$min)
  }
  
  return(list(mean(m_ests), BC_dim_local, BC_dim_meta))
}