Merge branch 'develop'

R/DAMOCLES_ML.R

@@ -94,6 +94,8 @@ DAMOCLES_all_loglik_choosepar <- function(trparsopt,
 #' package 'expm' or any of the numerical solvers, used in deSolve.
 #' @param model model used. Default is 0 (standard null model). Other options are 1 (binary traits)
 #' 2 (trinary environmental trait) or 3 (diversity-dependent colonization - beta version)
+#' @param num_cycles the number of cycles of opimization. If set at Inf, it will
+#' do as many cycles as needed to meet the tolerance set for the target function.
 #' @param verbose Whether intermediate output should be printed. Default is FALSE.
 #' @author Rampal S. Etienne
 #' @seealso \code{\link{DAMOCLES_loglik}} \code{\link{DAMOCLES_sim}}
@@ -115,6 +117,7 @@ DAMOCLES_ML <- DAMOCLES_all_ML <- function(
    edgeTList = NULL,
    methode = 'analytical',
    model = 0,
+   num_cycles = 1,
    verbose = FALSE)
 {
   locatenode <- NULL
@@ -211,7 +214,8 @@ DAMOCLES_ML <- DAMOCLES_all_ML <- function(
                            methode = methode,
                            pchoice = pchoice,
                            model = model,
-                           verbose = verbose)
+                           verbose = verbose,
+                           num_cycles = num_cycles)
       if(out$conv > 0)
       {
         cat("Optimization has not converged. Try again with different starting values.\n")

R/HERACLES_bootstrap.R

@@ -1,4 +1,4 @@
-HERACLES_find_trait_stait_combinations = function(pa,paTrait)
+HERACLES_find_trait_state_combinations = function(pa,paTrait)
 {
 	Np1 = length(which(pa  ==  1 & paTrait  ==  1))
 	Np2 = length(which(pa  ==  1 & paTrait  ==  2))
@@ -38,35 +38,98 @@ o.dectobin = function(y,ly = 0)
 
 combsUse = function(nRegional,nSample = 1000000)
 { 
-	nposs = 2^nRegional
-	nposs = nposs - 1
+	nposs = 2^nRegional - 1
 	combs = DDD::sample2(0:nposs,size = nSample,replace = TRUE)
 	return(combs)
 }	
 
 ######################################################################################################################################
 ######################################################################################################################################
 
-#this now includes the option to calculate trait based metrics of phylogenetic structure by setting traitMetric = TRUE and including traitdist, a trait distance matrix
-Heracles_ImportanceSampling = function(nSamples,n,regionalSpecies,S_regional,p,pa,phy,phydist,parsDAM,Mlist,model,pchoice,samptype,edgeObj,traitMetric,traitdist)
+#' Importance Sampling from the HERACLES model
+#' 
+#' Computes likelihood and metrics for randomly sampled presence-absence data of species in a local
+#' community for a given phylogeny of species in the region.
+#' 
+#' @param nsamples The number of samples used in importance sampling
+#' @param n
+#' @param regionalSpecies The list of species present in the regional community (SP)
+#' @param S_regional The number of species in the regional species pool
+#' @param p The probability used for the binomial distribution
+#' @param pa presence-absence table with the first column the species labels
+#' and the second column the presence (1) or absence (0) of the species
+#' @param phy phylogeny in phylo format
+#' @param phydist TBD
+#' @param parsDAM Vector of model parameters:\cr
+#' \code{pars[1]} corresponds to mu (extinction rate in local community)\cr
+#' \code{pars[2]} corresponds to gamma_0 in formula
+#' gamma(t) = gamma_0/(1 + gamma_1 * t) where gamma(t) is immigration rate
+#' into local community)\cr 
+#' \code{pars[3]} corresponds to
+#' gamma_1 in formula gamma(t) = gamma_0/(1 + gamma_1 * t) where gamma(t) is
+#' immigration rate into local community)
+#' @param Mlist list of M matrices that can be specified when methode = 'analytical'. If set
+#' at NULL (default) and methode = 'analytical', Mlist will be computed.
+#' @param model model used. Default is 0 (standard null model). Other options are 1 (binary traits)
+#' 2 (trinary environmental trait) or 3 (diversity-dependent colonization - beta version)
+#' @param pchoice sets the p-value to optimize:\cr
+#' pchoice == 0 corresponds to
+#' the sum of p_0f + p_1f\cr
+#' pchoice == 1 corresponds to p_0f\cr
+#' pchoice == 2 corresponds to p_1f\cr
+#' @param samptype Type of sampling distribution, can be either 'uniform' or 'binomial' in which case
+#' the local samples are uniformly or binomially generated, with the local diversity being a stochastic
+#' variable, or 'fixed' in which case the observed local diversity is used and configurations consistent
+#' with this diversity are sampled
+#' @param edgeObj list of edge lengths that need to be successively pruned; if
+#' not specified, it will computed using compute_edgeTList
+#' @param methode method used to solve the ODE. Either 'analytical' for the analytical
+#' solution, 'Matrix' for matrix exponentiation using package Matrix or 'expm' using
+#' package 'expm' or any of the numerical solvers, used in deSolve.
+#' @return A list containing attributes of the loglikelihood and importance sampling, and
+#' of the metrics (mntd and mpd, and TBD)
+#' @author Rampal S. Etienne & Alex L. Pigot
+#' @references Pigot, A.L. & R.S. Etienne (2015). A new dynamic null model for
+#' phylogenetic community structure. Ecology Letters 18: 153-163.
+#' @keywords models
+#' @examples TBD
+#' @export HERACLES_ImportanceSampling
+Heracles_ImportanceSampling <- function(nSamples,
+                                        n,
+                                        regionalSpecies,
+                                        S_regional,
+                                        p = n/S_regional,
+                                        pa,
+                                        phy,
+                                        phydist,
+                                        parsDAM,
+                                        Mlist = NULL,
+                                        model,
+                                        pchoice,
+                                        samptype,
+                                        edgeObj = NULL,
+                                        methode = 'analytical',
+                                        traitdist = NULL)
 {	
-	#create a matrix to store metrics of community structure
+  edgeObj <- DAMOCLES_check_edgeTList(phy,edgeObj)
+  Mlist <- DAMOCLES_check_Mlist(Mlist,parsDAM,model,methode)
+  #create a matrix to store metrics of community structure
 	#create a matrix to store the loglikelihood of each community and its sampling probability  
-	loglikMatrix = matrix(ncol = 1 + 2 * (samptype  ==  'binomial'),nrow = nSamples)
-	metricMatrix = matrix(ncol = 5,nrow = nSamples)
+	loglikMatrix <- matrix(ncol = 1 + 2 * (samptype == 'binomial'),nrow = nSamples)
+	metricMatrix <- matrix(ncol = 5,nrow = nSamples)
 
 	#I have now included the possibility of calculate two trait based metrics
 	if(samptype == 'binomial')
   {
-     colnames(loglikMatrix) = c("P_i","X_i","n")
+     colnames(loglikMatrix) <- c("P_i","X_i","n")
   } else
   {
-     colnames(loglikMatrix) = c("P_i")
+     colnames(loglikMatrix) <- c("P_i")
   }
-  colnames(metricMatrix) = c("mnpd","mpd","mntd","mtd","n")
+  colnames(metricMatrix) <- c("mnpd","mpd","mntd","mtd","n")
 
-	DAMOCLES.samp = matrix(rep(0,n),nrow = 1)
-	colnames(DAMOCLES.samp) = phy$tip.label
+	DAMOCLES.samp <- matrix(rep(0,n),nrow = 1)
+	colnames(DAMOCLES.samp) <- phy$tip.label
 
 	#I ADDED THESE IF STATEMENTS BECAUSE IF THE REGIONAL POOL IS  LARGE THEN WE CANNOT USE combsUse BECAUSE WE HAVE TOO MANY POSSIBLE CONFIGURATIONS
 	#FOR LARGE REGIONAL POOLS WE WILL HAVE TO USE BINOMIAL SAMPLING
@@ -78,60 +141,69 @@ Heracles_ImportanceSampling = function(nSamples,n,regionalSpecies,S_regional,p,p
   		 stop("The regional species pool is large. Uniform sampling will require exploring many community configurations and may exceed memory limitations. Use binomial sampling instead")
 		} else
     {
-			 combs = combsUse(S_regional,nSamples)
+			 combs <- combsUse(S_regional,nSamples)
 		}
 	}
+
+	if(samptype == 'fixed')
+	{
+	  S_loc <- sum(as.numeric(pa[,2]))
+	}
 
- 	pafoc = pa
-	pafoc[,2] = as.character(rep(0,n))
+ 	pafoc <- pa
+	pafoc[,2] <- as.character(rep(0,n))
 
 	for(i in 1:nSamples)
   {	
 	#if samptype is binomial, sample a local community richness S_loc from a binomial distribution with parameters S_regional and p
 	 	 if(samptype == 'binomial')
      {
-       	S_loc = stats::rbinom(n = 1,size = S_regional,prob = p)
-     		sam = sample(c(rep(1,S_loc),rep(0,S_regional - S_loc)))
+       	S_loc <- stats::rbinom(n = 1,size = S_regional,prob = p)
+     		sam <- sample(c(rep(1,S_loc),rep(0,S_regional - S_loc)))
   		 	#calculate the sampling probability of the sampled community under a binomial distribution with parameters S_regional and p
 
-	  	 	#the log probability of having local community richness S_loc is
-		   	loglikMatrix[i,2] = log(stats::dbinom(S_loc, size = S_regional, prob = p)) - (stats::dbinom(0, size = S_regional, prob = p) + stats::dbinom(1, size = S_regional, prob = p))
+	  	 	#the log probability of having local community richness S_loc conditional on this being larger than 1 is
+		   	loglikMatrix[i,2] <- stats::dbinom(S_loc, size = S_regional, prob = p, log = TRUE) - stats::pbinom(1, size = S_regional, prob = p, lower.tail = FALSE, log.p = TRUE)
 
 			  #the number of configurations with local richness S_loc is
-		  	loglikMatrix[i,3] = lgamma(S_regional + 1) - lgamma(S_loc + 1) - lgamma(S_regional - S_loc + 1)
+		  	loglikMatrix[i,3] <- lgamma(S_regional + 1) - lgamma(S_loc + 1) - lgamma(S_regional - S_loc + 1)
 		 } else
 	   if(samptype == 'uniform')
      {
-    		sam = o.dectobin(combs[i],S_regional)
-      	S_loc = sum(as.numeric(pa[,2]))
-     }  
+    		sam <- o.dectobin(combs[i],S_regional)
+    		S_loc <- sum(sam)
+	   }  
+	   if(samptype == 'fixed')
+	   {
+  	    sam <- sample(c(rep(1,S_loc),rep(0,S_regional - S_loc))) 
+	   }
 
 	   #enter the sampled local community into our pa dataframe 	
-		 pafoc[regionalSpecies,2] = sam
-	   DAMOCLES.samp[1,] = pafoc[,2]
+		 pafoc[regionalSpecies,2] <- sam
+	   DAMOCLES.samp[1,] <- pafoc[,2]
 
 		 #calculate phylogenetic metric for the sampled community	
-		 metricMatrix[i,1] = picante::mntd(DAMOCLES.samp,phydist)
-		 metricMatrix[i,2] = picante::mpd(DAMOCLES.samp,phydist)
+		 metricMatrix[i,1] <- picante::mntd(DAMOCLES.samp,phydist)
+		 metricMatrix[i,2] <- picante::mpd(DAMOCLES.samp,phydist)
 
-		 if(traitMetric == TRUE)
+		 if(!is.null(traitdist))
      {
         #calculate trait metrics for the sampled community	
-		    metricMatrix[i,3] = mntd(traitdist,pafoc)
-			  metricMatrix[i,4] = mtd(traitdist,pafoc)
+		    metricMatrix[i,3] <- DAMOCLES_mntd(traitdist,pafoc)
+			  metricMatrix[i,4] <- DAMOCLES_mtd(traitdist,pafoc)
 	   }
 
-		 metricMatrix[i,5] = S_loc
+		 metricMatrix[i,5] <- S_loc
 
 		 #calculate the loglikelihood of the sampled community under DAMOCLES
-		 loglikMatrix[i,1] = DAMOCLES_all_loglik(phy = phy,pa = pafoc,pars = parsDAM,pchoice = pchoice,edgeTList = edgeObj, methode = 'analytical', model = model,Mlist = Mlist)
+		 loglikMatrix[i,1] <- DAMOCLES_all_loglik(phy = phy,pa = pafoc,pars = parsDAM,pchoice = pchoice,edgeTList = edgeObj, methode = methode, model = model,Mlist = Mlist)
 
 		#print(i)
 	}	
 
-	resultsList = list()
-	resultsList[[1]] = metricMatrix
-	resultsList[[2]] = loglikMatrix
+	resultsList <- list()
+	resultsList[[1]] <- metricMatrix
+	resultsList[[2]] <- loglikMatrix
 	return(resultsList)
 }
 
@@ -183,13 +255,13 @@ HERACLES_extractCI = function(loglikMatrix,metricMatrix,ci_lower = 0.025,ci_uppe
 #functions to calculate metrics of trait structure. These are used in "Heracles_ImportanceSampling"
 
 #mean nearest trait distance
-mntd = function(traitdist,pa)
+DAMOCLES_mntd = function(traitdist,pa)
 {
 	present = which(pa[,2] == "1")
 	if(length(present)>1)
   {
 		traitdistfoc = traitdist[present,present]
-		mntd = mean(matrixStats::rowMins(traitdistfoc,na.rm=TRUE))
+		mntd = mean(matrixStats::rowMins(traitdistfoc,na.rm = TRUE))
 	} else {
 		mntd = NA
 	}	
@@ -198,13 +270,13 @@ mntd = function(traitdist,pa)
 
 #mean trait distance
 
-mtd = function(traitdist,pa)
+DAMOCLES_mtd = function(traitdist,pa)
 {
 	present = which(pa[,2] == "1")
 	if(length(present)>1)
   {
 		traitdistfoc = traitdist[present,present]
-		mtd = mean(rowMeans(traitdistfoc,na.rm=TRUE))
+		mtd = mean(rowMeans(traitdistfoc,na.rm = TRUE))
 	} else {
 		mtd = NA
 	}	

R/HERACLES_sim.R

@@ -71,7 +71,7 @@ root.state, root.trait.state, plotit = FALSE,keepExtinct = FALSE)
      #transition rates
      rates = c(mu_1,mu_2,gamma_1,gamma_2,q_0_to_2,q_2_to_0,q_1_to_2,q_2_to_1)
 
-     traitStateComb = HERACLES_find_trait_stait_combinations(pa,paTrait)
+     traitStateComb = HERACLES_find_trait_state_combinations(pa,paTrait)
 
      #gillespie algorithm    
      rates.real = (rates * traitStateComb)

tests/testthat/test_DAMOCLES.R

@@ -2,7 +2,7 @@
 #' 
 #' context("DAMOCLES")
 
-test_that("DAMOCLES_ML", {
+test_that("DAMOCLES_ML works", {
   data(NWPrimates_data)
   out <- DAMOCLES_ML(
     phy = NWPrimates_data[[1]],