Statistical models and computational algorithms for discovering relationships in microbiome data

3.1.1 The Dirichlet-Multinomial distribution

The multinomial distribution is suitable for representing counts from a population where the probability of selecting a category, such as an OTU, remains (approximately) the same, independent of other selections, and the total number of selections is fixed. The multinomial probability mass function of the multinomial is given in Eq. (1)).

Here, x_j is the number of times OTU j is observed among p different types of OTUs, ϕ_j is the probability of observing OTU j, and Γ is the gamma function. Note that we will use the common notation y+=∑jyj, yi+=∑jyij, etc. throughout this paper when the indices are understood by the context, so in Eq. (1), x₊ is the total number of OTUs observed in the sample.

To account for overdispersion, proportions may be considered to be random variables rather than fixed values. In this case, the distribution to model proportions is often the conjugate prior of the multinomial, the Dirichlet distribution, shown in Eq. (2),

where α_j > 0 is a parameter dictating the dispersion of the proportion ϕ_j, and so it itself is a parameter associated with OTU j. A compound distribution results by multiplying the two and integrating out the random proportions over the unit p − 1 simplex, resulting in the DM of Eq. (3).

3.1.2 Finite mixture models

A finite mixture model (FMM) is a method of representing G distinct subpopulations, each of which can possess its own probability density (or mass) function, f_g, for g = 1, 2, …, G. The probability function of the entire population can then be expressed as f(x)=∑g=1Gπgfg(x;𝚯g) where x = (x₁, x₂, …, x_p) are the data, π_g is the proportion of the population in group g called the mixing proportion, and 𝚯_g is the set of distribution parameters that correspond to f_g. The flexibility of a finite mixture model can be thought of as an additional means of accommodating the heterogeneity of the data that cannot simply be attributed to covariates, such as through a linear model, and can even be superior to them (McLachlan and Peel, 2000). Finite mixtures of the DM were developed previously for microbiome analysis (Holmes et al., 2012) to identify subgroups based on abundance, similar to previously considered enterotypes (Arumugam et al., 2011). Parameter and group memberships of finite mixture models are typically solved using the Expectation-Maximization (EM) algorithm (Dempster et al., 1977) or one of its variants.

3.1.3 Evolutionary algorithms

An evolutionary algorithm (EA) is a means of finding optimal solutions to a problem that mimics natural selection whereby organisms (possible solutions) with better genes (properties) are favoured over successive generations (iterations). A didactic example of an EA is the string evolver. Imagine one person (evaluator) has a seven-digit phone number in mind and another person (guesser) is trying to guess it but only knows the phone number is seven digits. The guesser begins by giving a collection of 10 randomly-generated seven-digit numbers to the evaluator who ranks the numbers from most correct to least correct and tells the guesser the ranking of each phone number based on how correct it is. The guesser discards the five worst phone numbers, replacing them with a copy of each of the five best. The guesser then has two copies of each of the five phone numbers, so one random digit in one of the two copies is changed to a random number. Now the guesser has 10 numbers which could all be distinct. The process of the evaluator rating the phone numbers repeats with the five best phone numbers being kept, copied, etc.. Eventually, the same five phone numbers will be returned in several successive iterations. These five phone numbers will be potential solutions the guesser has to the phone number the evaluator has in mind. Despite being random, this process of guessing a phone number is far more efficient than having the guesser ask the evaluator to individually rate all ten million possible random phone numbers (the brute-force approach).

As mentioned, EAs are inspired by the process of propagation and selection based on fitness, found in nature. Although we do not provide a complete summary of these class of algorithms, we briefly outline EAs and highlight some of the terminology (Ashlock, 2006) which was also the source motiving our string evolver example above. The reader is referred to this book or expert overviews and further reference (Bäck, 1996; Bäck and Schwefel, 1996) on EAs.

There are several key components to an EA which we will define here (and parenthetically illustrate from the string evolver example). All of these terms are in reference to the EA, not the biological context. Although not common in EA literature, we will append the word “evolutionary” or “EA” in front of EA terms where necessary to distinguish them from their biological counterparts. An evolutionary gene is a trait that can differ between solutions (one digit of the phone-number). An evolutionary organism is one potential solution comprised of genes (a complete phone number). Variational operators are methods of constructing new phone numbers. The only variational operator we considered in the phone number-example was mutation: altering exactly one gene in an organism (altering one digit). Offspring are the organisms resulting from the application of variational operators (new phone numbers resulting from copying old ones and changing a digit). Fitness is some value assigned to an organism for the purpose of selection, pruning away organisms based on fitness. In our example, the evaluator performed the task of evaluating the fitness of individuals for the purpose of pruning away the lowest five scoring phone numbers.

In our microbiome context, we will propose an EA where compositional similarities, represented by constraints in the statistical model, are the EA organisms and a standard statistical model-selection criterion is our method of selection. This is akin to stepwise selection in multiple linear regression, but stepwise selection is a special case of our method. Rather than constraints, the goal of stepwise selection is to determine a subset of potential predictors in a model by successively adding a predictor at every iteration. In stepwise selection, the case is special because every possible mutation (adding/removing a potential predictor from the model) is considered at every iteration, and only the single best organism (potential model with predictors) survives selection (based on some criterion like Bayesian Information Criterion; BIC) to form a new generation. This reduces the random process to a deterministic one once the algorithm begins in stepwise regression. However, the more general undermentioned algorithm should be considered stochastic.

3.2.1 EA construction

We will represent constraints through partitions of parameters where each partition has only one degree of freedom. If two parameters exist in the same partition, then their value is necessarily forced to be equal to each other. In an unconstrained model, every parameter is in a partition by itself. Every partition can be represented by sets of equalities between every pair of parameters. Thus, the EA organism, a particular set of constraints on the DM, is faithfully represented as a collection of EA genes: where an EA gene is a pairwise constraint. A mutation in the EA will be imposing a new constraint to the EA organism and we avoid null mutations by ensuring we only allow mutations which can change the organism. Note that in some instances, some additional equalities will be imposed by transitivity. For instance if α₁ = α₂ already in one candidate model at one stage of the EA and the constraint α₂ = α₃ is imposed on the model, the constraint α₁ = α₃ must automatically be imposed. All three OTUs represented by α₁, α₂ and α₃ will be considered to exist in approximately the same proportions in the microbiome.

When a constraint is imposed, rather than replace the value of one parameter another and completely lose the information originally held by one parameter, we will replace both by the average of the two. In general, when any partition is formed, each parameter in the partition will be replaced by the arithmetic mean of original values of the parameters in the partition. This is intuitive but has a useful property specific to the DM in complete-data log-likelihood shown in Eq. (4).

Here, α_gj is the parameter corresponding to OTU j in group g where the proportion of group g comprised of OTU j is αgjαg+, z_ig represents the probability that observation i belongs in group g of the FMM, also called a responsibility, and k is a constant with respect to the DM parameters. We can see that the first pair of terms only involve α_g+. If we change the value of one α_gj but alter others so that the value of α_g+ remains the same, this aspect of the log-likelihood remains unchanged. Our goal will be to find ways to do this as much as possible without unduly detrimenting the log-likelihood by finding parameters that are very close to each other and replacing them with their average, thereby keeping the value of α_g+ the same. This concept of maintaining the sum of the α_g+ values is drawn from the theoretical property of complete neutrality (Mosimann, 1962). This property states that the removal of one variable in the model will not impact the values of the other parameters in the model, and so their relative ratios will remain stable. This will also be the case in our approach: altering the values of parameters as we describe will not influence other parameters nor our perception of the OTU’s composition in the microbiome. Hence, this computational property reflects a very important interpretative property.

Our goal in finding compositional similarities in the microbiome is reflected by sufficiently similar parameters in the model which we will constrain to be equal. We make two notes which aids the development of our algorithm. First, any final collection of constraints can be decomposed into a set of pairwise constraints. Second, there exists a sequence (not necessarily unique) of successive pairwise constraints that can be added incrementally to an unconstrained model to arrive at any constrained model. Therefore, we will focus only on pairwise constraints and represent the final set of constraints as sets of pairwise constraints. Model similarity is again like the case of the t distribution. Two t distributions with a large difference in degrees of freedom can be considered to be approximately the same distribution when the degrees of freedom parameter from each respective distribution is large. The same difference in degrees of freedom, however, may not be acceptable to consider approximate equality if one of the degrees of freedom is small. For instance, the difference between t distributions with 105 and 115 degrees of freedom is often negligible but the difference between t distributions with 5 and 15 degrees of freedom is far more considerable. By way of this analogy, it should then make sense to declare two distributions sufficiently similar when the reciprocal difference is sufficiently close to zero. This shall be our approach to compare different distributions that differ only by one of the α parameters.

Given a candidate set of constraints, the complete-data log-likelihood is re-evaluated and the effective number of parameters in light of the constraints re-counted to estimate the BIC (Schwartz, 1978). We use the BIC because it is well understood and enforces fewer free parameters, which we prefer, than than other well-known criteria like the AIC. It has also specifically been favoured in the context of FMMs (Leroux, 1992; Keribin, 2000) for many reasons including as a model selection criterion ((Fraley and Raftery, 1998, 2002).