Estimating Probabilities via Maximum Entropy: Improved Iterative Scaling
Description
maxent
returns the probabilities that maximize the entropy conditional on a series of constraints that are linear in the features. It relies on the Improved Iterative Scaling algorithm of Della Pietra et al. (1997). It has been used to predict the relative abundances of a set of species given the trait values of each species and the communityaggregated trait values at a site (Shipley et al. 2006; Shipley 2009; Sonnier et al. 2009).
Usage
1 
Arguments
constr 
vector of macroscopical constraints (e.g. communityaggregated trait values). Can also be a matrix or data frame, with constraints as columns and data sets (e.g. sites) as rows. 
states 
vector, matrix or data frame of states (columns) and their attributes (rows). 
prior 
vector, matrix or data frame of prior probabilities of states (columns). Can be missing, in which case a maximally uninformative prior is assumed (i.e. uniform distribution). 
tol 
tolerance threshold to determine convergence. See ‘details’ section. 
lambda 
Logical. Should lambdavalues be returned? 
Details
The biological model of community assembly through traitbased habitat filtering (Keddy 1992) has been translated mathematically via a maximum entropy (maxent) model by Shipley et al. (2006) and Shipley (2009). A maxent model contains three components: (i) a set of possible states and their attributes, (ii) a set of macroscopic empirical constraints, and (iii) a prior probability distribution q = [qj].
In the context of community assembly, states are species, macroscopic empirical constraints are communityaggregated traits, and prior probabilities q are the relative abundances of species of the regional pool (Shipley et al. 2006, Shipley 2009). By default, these prior probabilities q are maximally uninformative (i.e. a uniform distribution), but can be specificied otherwise (Shipley 2009, Sonnier et al. 2009).
To facilitate the link between the biological model and the mathematical model, in the following description of the algorithm states are species and constraints are traits.
Note that if constr
is a matrix or data frame containing several sets (rows),
a maxent model is run on each individual set. In this case if prior
is a vector,
the same prior is used for each set. A different prior can also be specified for each set.
In this case, the number of rows in prior
must be equal to the number of rows in constr
.
If q is not specified, set pj = 1 / S for each of the S species (i.e. a uniform distribution), where pj is the probability of species j, otherwise pj = qj.
Calulate a vector c = [ci] = {c1, c2, ..., cT}, where ci = sum(tij); i.e. each ci is the sum of the values of trait i over all species, and T is the number of traits.
Repeat for each iteration k until convergence:
1. For each trait ti (i.e. row of the constraint matrix) calculate:
gamma_i(k) = ln(t.mean_i / [sum(pj(k) tij)]) (1 / ci)
This is simply the natural log of the known communityaggregated trait value to the calculated communityaggregated trait value at this step in the iteration, given the current values of the probabilities. The whole thing is divided by the sum of the known values of the trait over all species.
2. Calculate the normalization term Z:
Z(k) = sum(pj(k) e^(gamma_i(k) tij) )
3. Calculate the new probabilities pj of each species at iteration k+1:
pj(k+1) = [pj(k) e^(gamma_i(k) tij)/ Z]
4. If max(pj(k+1)  pj(k)) <= tolerance threshold (i.e. argument tol
) then stop, else repeat steps 1 to 3.
When convergence is achieved then the resulting probabilities (pj_hat) are those that are as close as possible to qj while simultaneously maximize the entropy conditional on the communityaggregated traits. The solution to this problem is the Gibbs distribution:
Note: equation not shown in HTML help file: please refer to PDF manual.
This means that one can solve for the Langrange multipliers (i.e. weights on the traits, lamda_i) by solving the linear system of equations:
Note: equation not shown in HTML help file: please refer to PDF manual.
This system of linear equations has T+1 unknowns (the T values of lambda plus ln(Z)) and S equations. So long as the number of traits is less than S1, this system is soluble. In fact, the solution is the wellknown least squares regression: simply regress the values ln(pj_hat of each species on the trait values of each species in a multiple regression.
The intercept is the value of ln(Z) and the slopes are the values of lambda_i and these slopes (Lagrange multipliers) measure by how much the ln(pj_hat), i.e. the ln(relative abundances), changes as the value of the trait changes.
maxent.test
provides permutation tests for maxent models (Shipley 2010).
Value
prob 
vector of predicted probabilities 
moments 
vector of final moments 
entropy 
Shannon entropy of 
iter 
number of iterations required to reach convergence 
lambda 
lambdavalues, only returned if 
constr 
macroscopical constraints 
states 
states and their attributes 
prior 
prior probabilities 
Author(s)
Bill Shipley bill.shipley@usherbrooke.ca
http://pages.usherbrooke.ca/jshipley/recherche/
Ported to FD by Etienne Laliberté.
References
Della Pietra, S., V. Della Pietra, and J. Lafferty (1997) Inducing features of random fields. IEEE Transactions Pattern Analysis and Machine Intelligence 19:113.
Keddy, P. A. (1992) Assembly and response rules: two goals for predictive community ecology. Journal of Vegetation Science 3:157164.
Shipley, B., D. Vile, and É. Garnier (2006) From plant traits to plant communities: a statistical mechanistic approach to biodiversity. Science 314: 812–814.
Shipley, B. (2009) From Plant Traits to Vegetation Structure: Chance and Selection in the Assembly of Ecological Communities. Cambridge University Press, Cambridge, UK. 290 pages.
Shipley, B. (2010) Inferential permutation tests for maximum entropy models in ecology. Ecology in press.
Sonnier, G., Shipley, B., and M. L. Navas. 2009. Plant traits, species pools and the prediction of relative abundance in plant communities: a maximum entropy approach. Journal of Vegetation Science in press.
See Also
functcomp
to compute communityaggregated traits,
and maxent.test
for the permutation tests proposed by Shipley (2010).
Another faster version of maxent
for multicore processors called maxentMC
is available from Etienne Laliberté (etiennelaliberte@gmail.com). It's exactly the same as maxent
but makes use of the multicore, doMC, and foreach packages. Because of this, maxentMC
only works on POSIXcompliant OS's (essentially anything but Windows).
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  # an unbiased 6sided dice, with mean = 3.5
# what is the probability associated with each side,
# given this constraint?
maxent(3.5, 1:6)
# a biased 6sided dice, with mean = 4
maxent(4, 1:6)
# example with tussock dataset
traits < tussock$trait[, c(2:7, 11)] # use only continuous traits
traits < na.omit(traits) # remove 2 species with NA's
abun < tussock$abun[, rownames(traits)] # abundance matrix
abun < t(apply(abun, 1, function(x) x / sum(x) )) # relative abundances
agg < functcomp(traits, abun) # communityaggregated traits
traits < t(traits) # transpose matrix
# run maxent on site 1 (first row of abun), all species
pred.abun < maxent(agg[1,], traits)
## Not run:
# do the constraints give predictive ability?
maxent.test(pred.abun, obs = abun[1,], nperm = 49)
## End(Not run)
