Because demographic rates contain probabilities, demographic outcomes are random variables. Because a heterogeneous population is a mixture, the distribution of the outcome is a mixture of the distributions within each group. This mixture is characterised by the moments of each group and the distribution of individuals among the groups. The latter distribution is called the mixing distribution.

The variance in outcome in a heterogeneous population reflects both stochasticity within each group and heterogeneity among groups (remember that the groups are homogeneous, and thus differences in outcome within a group are due to are strictly due to chance). In any outcome calculated from a life table, a Markov chain, or some equivalent machinery, individual stochasticity is the only source of within-group variance, because the calculation explicitly applies the same probabilities to every individual at every age. See Caswell (2023) and references therein for a discussion of the roles of heterogeneity and stochasticity and the issues of interpretation that arise from them.

Variances in outcomes are often called “inequalities,” although that identification is more subtle than is usually appreciated (Caswell, 2023). Not all differences are inequalities.¹ In economic terms (e.g., Atkinson, 2015), variance due to heterogeneity corresponds to inequality of opportunity. Individuals in, e.g., different income groups have different opportunities to live a long life. Variance due to stochasticity corresponds to inequality of outcome. Individuals in the same income group may differ in longevity because of the random outcomes of mortality risks, even though they all experience the same group-specific risks.

The relative magnitude of the contributions of stochasticity and heterogeneity depend on many factors, including the type and magnitude of the heterogeneity. Thus it is essential to partition the variance into within-group and between-group components.

Our intent is not to review the concepts of heterogeneity, stochasticity, and inequality, as a previous paper has explored these issues in some detail (Caswell, 2023). That paper surveyed the contributions of heterogeneity (among groups) and individual stochasticity (with groups) to: (1) variance in human longevity among latent frailty classes, income groups, education groups, and neighbourhood deprivation groups; (2) variance in adult longevity among insects differing in early life nutrition; (3) variance in longevity and lifetime reproduction in a wild seabird population; (4) variance in longevity and lifetime reproduction in laboratory populations of a rotifer, with groups defined by maternal age; (5) variance in healthy longevity with groups defined by European countries; and (6) variance in longevity and lifetime reproduction with groups defined by species and/or populations of plants and animals.

Instead, our goal is to provide the calculations necessary to extend variance partitioning from single factors to multiple interacting factors. Studies of the variance in longevity and in lifetime reproduction typically define heterogeneity in terms of a single factor. When studies examine multiple factors, they are usually treated one at a time. However, individuals are heterogeneous in multiple factors operating simultaneously, and the resulting variance is affected by all of those factors and their interactions. The inability to evaluate the contributions of interactions is a major limitation to the study of heterogeneity, and our goal here is to show how to analyse multi-factor studies, in which individuals are heterogeneous in two or more factors. We will present some two-factor and three-factor examples of the calculation. There is no limitation to the number of factors.

The variance in outcome among a heterogeneous set of individuals can be partitioned into between-group and within-group components using classical results from conditional probability. Consider a heterogeneous population that is a mixture of groups, the relative abundances of which are given by a mixing distribution π. Let ξ denote some outcome. The variance in ξ is V(ξ)=Eπ[V(ξ|group)]+Vπ[E(ξ|group)](1) =Vwithin(ξ)+Vbetween(ξ).(2)The within-group variance Vwithin is the expectation of the variance within each group, weighted by the distribution π. The between-group variance Vbetween is the variance among the group means, again weighted by the distribution π. This variance decomposition is a standard result in conditional probability (e.g., Rényi, 1970; Frühwirth-Schnatter, 2006) and in analysis of variance (ANOVA) (e.g., Kempthorne, 1957). Indeed, the decomposition was introduced, along with the term “variance” itself, by Fisher (1918). This analysis has now been used repeatedly to examine demographic inequality among groups defined by a single factor (a “one-way” design).

We will show results in this paper for variance decomposition for longevity and lifetime reproductive output (LRO). The means and variances (and other moments) of longevity are readily calculated by expressing the demographic rates in terms of absorbing Markov chains (pioneered by Feichtinger, 1971; see Caswell, 2001, 2006, 2009). It is now known that, as a general rule, even extreme differences among groups, including those created by important socioeconomic variables, account for only a small fraction of the variance in longevity (see the overview in Caswell, 2023).

The means and variances (and other moments) of LRO can be calculated using Markov chains with rewards (Caswell, 2011; van Daalen and Caswell, 2017). Stochasticity in LRO reflects both survival (which determines how long a woman has to reproduce) and fertility (which determines whether or not she reproduces at each age). These two components can be separated, and in developed, low-mortality countries the variance in LRO is increasingly accounted for by stochasticity in fertility (van Daalen and Caswell, 2015). However, the contributions of heterogeneity and stochasticity to the variance in LRO are not yet well understood. In one case, in which the groups are defined by maternal age in a laboratory population, the variance components depend strongly on the environmental conditions (van Daalen et al., 2022). We will provide some examples below, in Section Lifetime reproductive output.

The following notation is used throughout this paper. Matrices are denoted by upper-case bold characters (e.g., U) and vectors are denoted by lower-case bold characters (e.g., a). Vectors are column vectors by default; xT is the transpose of x. The vector 1 is a vector of ones, and the matrix I is the identity matrix. When necessary, subscripts are used to denote the size of a vector or matrix; e.g., Iω is an identity matrix of size ω×ω. Matrices and vectors with a tilde (e.g., U˜ or a˜) are block-structured; in this paper, blocks correspond to different factors. The notation ∥x∥ denotes the 1-norm of x (i.e., the sum of the absolute values of the entries). The symbol ⊗ denotes the Kronecker product. The vec operator stacks the columns of a m×n matrix into a mn×1 column vector. When applied to an array with more than two dimensions, it stacks columns from all dimensions. We will make use of a reshape operator that is the inverse of the vec operator, changing the vector back into the array; see equation (17). On occasion, Matlab notation will be used to refer to the entries, rows, and columns of matrices. For example, F(i,j) is the (i,j) entry of F, while F(i,:) and F(:,j) refer to the ith row and jth column of the matrix.

Element-by-element operations apply to arrays of any dimension. The symbol ○ denotes the Hadamard, or element-by-element product (implemented by .* in Matlab and by * in R). The symbol ⊘ is used to denote the Hadamard, or element-by-element quotient. Thus, for two matrices A and B, A○B=(aijbij)A⊘B=(aijbij)(3)with the obvious restrictions that the objects must be of the same size and, that for the quotient, none of the entries of B can be zero.

Let x be a vector of numbers and π a probability vector of the same length. The mean and variance of the entries of x, over the mixing distribution π are, in matrix notation, Eπ(x)=πTx(4) Vπ(x)=πT(x○x)−(πTx)2.(5)The second of these terms, the variance of X over a mixing distribution, will appear so commonly that we define it as a function V(x, π) that returns the variance, over the distribution π, of the vector x: V(x,π)≡πT(x○x)−(πTx)2.(6)

Let us review the simplest case: a one-factor design in which a population (or some other set of individuals) is divided into groups based on a single factor that we call A (e.g., income), with levels a=a1,a2,…,anA. The population is a mixture with a mixing distribution π. Let ξ be the demographic outcome of interest. The vectors containing the means and variances of ξ are m=(E(ξ|a1)⋮E(ξ|anA))=(m1⋮mnA)(7) v=(V(ξ|a1)⋮V(ξ|anA))=(v1⋮vnA)(8)In terms of these vectors, the variance decomposition in equation (2) becomes Vwithin=πTv(9) Vbetween=πT(m○m)−(πTm)2(10) =V(m,π).(11)The relative contribution of the within- and between-group components to the variance V(ξ) is measured by K=VbetweenVwithin+Vbetween(12)which measures the fraction of the total variance due to heterogeneity among groups. The ratio is referred to as the intraclass correlation coefficient in quantitative genetics (Falconer, 1960). Its square root is called the correlation ratio in probability theory (Rényi, 1970).

Treating each of the nAnB combination of factors A and B as a group, the within- and between-group variance components are calculated from the entries of M, V, and Π. The within-group variance is the mean, over the mixing distribution, of the variances in each group, Vwithin=(vec Π)Tvec V.(18)The between-group variance is given by applying equation (6) to the means of all the factor combinations: Vbetween=V(vec M,vec Π).(19)This between-group variance measures the overall contribution of heterogeneity, over all the groups, of all the factors, with individuals distributed according to Π. It makes no distinction between the contributions of A, B, and the AB interaction, but if the goal is a measure of how heterogeneity in factors contributes to inequality of outcomes, Vbetween is the answer.

In the two-factor case, Vbetween is due to contributions from each factor (VA,VB) and their interaction (VAB). These components are calculated from the marginal means and marginal mixing distributions corresponding to each factor and the interaction. Use subscripts to identify the marginal means (e.g., mA as the vector of means for each level of A, marginalising over levels of B), and let a,b denote the levels of A and B, respectively. The array of marginal means for factor A is obtained by calculating the average, weighted by the mixing distribution, over the levels of factor B. Recalling the definitions of the Hadamard product and the Hadamard quotient in Section Notation, we have mA=[∑b(M○Π)]⊘[∑b Π]nA×1.(20)The same pattern holds for the marginal mean for B: mB=[∑a(M○Π)]⊘[∑a Π]nB×1(21)and for the marginal mean for AB mAB=[M○Π]⊘ΠnA×nB.(22)The marginal mean for factor A is obtained by summing over levels of B; the marginal mean for factor B is obtained by summing over levels of A. The marginal mean for the combination AB is not really marginal, because there are no other factors over which to sum. The marginal means for factors A and B are one-dimensional vectors. The marginal mean for the AB interaction is a two-dimensional array (it is, in fact, just the array M).

The marginal mixing distributions are obtained from the array Π, πA=∑b ΠnA×1(23) πB=∑a ΠnB×1(24) πAB=ΠnA×nB.(25)

Matlab makes these calculations easy to implement; the command corresponding to equation (20), for example, is mA=sum(M .*Pi,[2]) ./sum(Pi,[2]).(26)

Finally, the variance components due to factors A and B and the interaction AB are calculated by applying the function V(·,·) defined in equation (6) to the marginal means and marginal mixing distributions: VA=V(mA,πA)(27) VB=V(mB,πB)(28) VAB=V(vec mAB,vec πAB)−VA−VB.(29)The marginal mean array mAB contains the effects of A and B as well as the interaction; thus the interaction variance VAB is obtained by subtracting VA and VB from the variance among all factor combinations.

The calculation of Vbetween treats vec M and vec V as vectors containing values for each of the groups. There is no restriction on the mixing distribution Π except, of course, that its entries sum to one.

Here, however, we want to partition Vbetween into its components, and to do this the mixing distribution Π requires some careful attention. Statistics texts are unanimous in stressing the importance of equal sample sizes in all treatment combinations in the analysis of variance (ANOVA) in a factorial experiment. The sample sizes play the role of the mixing distribution in our probability calculations. It has long been known (Yates, 1934) that unequal sample sizes make it impossible to calculate variance components. A recent text explains

“When the sample sizes for each cell are unequal, the two-way analysis of variance for factor effects becomes complex. The component sums of squares in the analysis of variance are no longer orthogonal; that is, they do not sum to the total sum of squares. The least squares method for obtaining the best estimates of the parameters is rather complicated in the fixed effects model and the best analysis has not been and probably will not be found for the random effects models” (Sahai and Ageel, 2012).

There exist two classes of mixing distributions that permit the calculation of variance components. One is the flat, or balanced distribution, which corresponds to equal sample sizes and assigns equal weight to all factor combinations. The other is the class of rank-one, or proportional mixing distributions (Yates, 1934; Kirk, 1982, and many others). Such a mixing distribution, when written as an array as in equation (15), has proportional rows and proportional columns. The flat distribution is a special case of the rank-one distribution.

A rank-one mixing distribution array can be assembled from its marginals. Let πA and πB be the marginal distributions among the levels of factors A, B. The mixing distribution array is Π=reshape(πB⊗πA,nA,nB).(30)Every column is proportional to πA and every row is proportional to πB. Note the order of the subscripts.

As in Section Factorial designs: Two factors, the within-group and between-group variances are calculated by treating all nAnBnC factor combinations as groups. Then, just as in the two-factor case, Vwithin=(vec Π)Tvec V(33) Vbetween=V(vec M,vec Π).(34)As in the two-factor case, Vbetween gives the contribution to variance of all the heterogeneity among groups, but no information on the contributions of the factors and their interactions.

The between-group variance Vbetween is partitioned into components due to the main effects of each factor (VA,VB,VC), the two-way interactions (VAB,VAC,VBC) between pairs of factors, and the three-way interaction (VABC). These components are calculated from the marginal means and marginal mixing distributions corresponding to each factor and each interaction.

The array of marginal means for factor A is obtained by averaging over the dimensions other than those for factor A: mA=[∑b∑c(M○Π)]⊘[∑b∑c Π]nA×1.(35)The same pattern holds for the other marginal means: mB=[∑a∑c(M○Π)]⊘[∑a∑c Π]nB×1(36) mC=[∑a∑b(M○Π)]⊘[∑a∑b Π]nC×1(37) mAB=[∑c(M○Π)]⊘[∑c Π]nA×nB(38) mAC=[∑b(M○Π)]⊘[∑b Π]nA×nC(39) mBC=[∑a(M○Π)]⊘[∑a Π]nB×nC(40) mABC=[(M○Π)]⊘ΠnA×nB×nC.(41)

The marginal mixing distributions are obtained from the array Π, πA=∑b∑c ΠnA×1(42) πB=∑a∑c ΠnB×1(43) πC=∑a∑b ΠnC×1(44) πAB=∑c ΠnA×nB(45) πAC=∑b ΠnA×nC(46) πBC=∑a ΠnB×nC(47) πABC=ΠnA×nB×nC.(48)

The calculations are readily expressed as Matlab commands; for example, the command corresponding to equation (35) is mA=sum(M .*Pi,[2,3]) ./sum(Pi,[2,3]).(49)The second argument ([2, 3]) in the sum command indicates the dimensions over which summation takes place.

The variance components due to each of the factors and interactions are calculated by applying the function V(·,·) in equation (6) to the vectors obtained by applying the vec operator to the marginal mean arrays VA=V(mA,πA)(50) VB=V(mB,πB)(51) VC=V(mA,πA)(52) VAB=V(vec mAB,vec πAB)−VA−VB(53) VAC=V(vec mAC,vec πAC)−VA−VC(54) VBC=V(vec mBC,vec πBC)−VB−VC(55) VABC=V(vec mABC,vec πABC)−VA−VB−VC−VAB−VAC−VBC.(56)The arrays for the two-way interactions also include the one-factor effects, so the one-factor variances are subtracted to obtain the two-factor variances. The variance due to the three-factor interaction has the one-factor and two-factor variances subtracted.

The interpretation of interactions in factorial experiments has always been a challenge (e.g. Steel and Torrie, 1960; Sahai and Ageel, 2012). A large component of the variance due to an AB interaction makes it difficult to say what the effects of A and B are because the effect of A depends on the level of B, and vice versa. The situation becomes even more difficult, of course, for three-way or higher interactions. If the contributions of interactions to variance are small, they can be ignored. Note that we have no operational definition of “small” such as is provided in ANOVA by tests of the statistical significance of the interactions. However, the difficulty of interpreting interactions does not change the fact that they are substantively interesting. Knowing that two factors interact is an important substantive finding and invites further study to understand how that interaction works.

Note that if the factors have additive effects on ξ, the interaction variances are all zero (Caswell, unpublished results).

Variance in longevity, often referred to as inequality or disparity in lifespan, has been analysed across a variety of social, economic, and biological variables. Partitioning of the variance into components has revealed that, even when group differences are very large, they contribute only a small fraction of the total variance (see an overview in Caswell, 2023). Because longevity is the outcome of a lifetime of repeated probabilistic survival events, it is subject to a large degree of individual stochasticity.

However, these analyses have been limited to single factors. Here, we report two examples of multi-factor studies: one examining variance due to sex and race, and the other examining sex, race, and state of residence. In these examples, we calculated the mean and variance of longevity using Markov chain methods (Feichtinger, 1971; Caswell, 2001, 2006), but they could equally have been calculated from a set of life tables.

Lifetime reproductive output (LRO), sometimes called lifetime reproductive success (LRS), is the number of offspring produced by a female over her lifetime. When offspring are measured as female children and are subject to mortality, the mean of LRO is the net reproductive rate R0 (or NRR). When offspring are measured as children of either sex and mortality is excluded, the mean of LRO is the total fertility rate (TFR). See Keilman et al. (2014) for a discussion of measures of lifetime reproduction including different combinations of sexes, and Caswell (2009, p. 1771) for the criteria that must be met for the net reproductive rate.

Lifetime reproductive “output” is a purposely flexible term. Output could be defined as children of either sex, or of both sexes, or of any other categories. Thus, it is important to specify how output is defined. Other categories can be analysed the same way. If, for example, we were interested in left-handed blue-eyed children and had age-specific numbers of such births per female, the lifetime output of such children could be calculated. Lifetime production of left-handed blue-eyed children is purposely frivolous, but the same logic applies to, for example, children born with trisomy-21 (the chromosome duplication that causes Down syndrome). The statistics of the lifetime production of children with trisomy-21 could be very interesting indeed.

Because survival and reproduction include stochastic processes, LRO is a random variable, with a variance, skewness and higher moments. This variance could be considered a “reproductive inequality” just as variance in longevity is considered “lifespan inequality.” However, this type of inequality of outcome seems to have drawn little interest (but see early work on the concentration of reproduction by Vaupel and Goodwin, 1987 and Shkolnikov et al., 2007). Variation in outcomes such as childlessness and completed parity has attracted attention, and integrating these outcomes into a variance partitioning framework is an interesting research challenge.

The variance among individuals in lifetime reproduction has been approached by population biologists in several ways (e.g., Caswell, 2011; Steiner and Tuljapurkar, 2012; van Daalen and Caswell, 2017; Snyder and Ellner, 2018; Tuljapurkar et al., 2020). Here we use Markov chains with rewards (MCWR), which provide the mean and variance (indeed, all the moments) of lifetime reproduction. See Caswell (2011); van Daalen and Caswell (2017, 2020) for details of the calculation. Briefly, individual development, survival and transitions among the possible life history states are described by an absorbing Markov chain. Rewards, in the form of offspring produced along the life history trajectory, are accumulated until death. The model provides the mean and variance of lifetime fertility as a function of the demographic rates, depending on the fate of the individual and on the outcomes of its chances of reproduction at each point in its lifetime.

While some comparisons of the variance, skewness and other statistics of LRO for humans and other species, have been reported (Caswell, 2011; van Daalen and Caswell, 2015; van Daalen and Caswell, 2017; van Daalen and Caswell, 2024; Varas Enríquez et al., 2022), only a few studies have partitioned variance in LRO into contributions from heterogeneity and stochasticity. One example is the study by van Daalen et al. (2022) of the contribution of heterogeneity in maternal age to variance in LRO in a rotifer.

As with longevity, our goal here is to extend one-factor analyses to multi-factor studies in which the variance in lifetime reproduction can be partitioned into contributions from heterogeneity and stochasticity. We show two examples. One is a comparison of lifetime fertility among human females in 31 developed countries over a 40-year time interval, treating country and time as factors. The other is a laboratory study of the effects of food limitation and pesticide exposure on a rotifer.

1.

Define the factors that characterise the heterogeneity among individuals (e.g., sex, race, income, education, resource levels) and the levels of each factor.

There are many ways to measure the variation (“inequality” in a broad sense) of some quantity. Economists have developed many indices to address specifically economic issues related to income, transfers, and so on (e.g., Sen, 1997; Jenkins and van Kerm, 2009; Atkinson, 2015). In demographic contexts, these measures are highly correlated (Van Raalte and Caswell, 2013). Thus if the aim is to focus only on the value of the index, there is little basis on which to choose among them. But if the aim is to go beyond the values, there are properties of the variance that make it an attractive choice.

As is well known, the variance is decomposable into within- and between-group components (as are indices based on entropy). Additive decomposability is more than a mathematical nicety; it is fundamental to our attempts to understand how heterogeneity contributes to inequality of outcomes. It was for this purpose that Fisher introduced it as the basis for the analysis of experimental data (Fisher, 1936). Statistical ANOVA has developed into an enormous variety of study designs, referred to as “experimental designs” in the experimental sciences. In our context, these designs correspond to arrangements of factors within and across populations. The factorial study design we explore here only scratches the surface of the possibilities.

The most popular alternative to the variance is the Gini coefficient. The Gini coefficient is based on the mean of the absolute values of the deviations between pairs of variates (e.g., individuals) selected from the distribution. Because of this, Permanyer et al. (2023) call the Gini coefficient an individual measure. They contrast it with the variance, which they call a group measure because the variance is usually written as the mean of the squared deviations from the mean, rather than as a difference among individual values. However, the variance and the Gini coefficient are both individual measures in this sense. The mean difference between two individuals from the distribution f(·) is Δ1=∫−∞∞∫−∞∞|x−y| f(x)f(y)dxdy.(65)The Gini coefficient is a standardised, dimensionless version of the mean difference G=Δ12μ1(66)where μ1 is the mean.

The variance can also be written as a difference among individuals, “without reference to deviations from a central value, the mean” (Kendall and Stuart, 1969, p. 47) as V=12∫−∞∞∫−∞∞(x−y)2 f(x)f(y)dxdy.(67)

The similarity between the mean difference and the variance is clear; the former uses absolute values of the differences between individuals, while the latter uses the squares of the differences. Just as the Gini coefficient is a standardised mean difference, the variance can be standardised as the familiar coefficient of variation CV=Vμ1(68)or the standardised variance Vs=Vμ12.(69)

The standardised variance is known in evolution and anthropology as Crow’s index of the opportunity for selection (Crow, 1958; Courtiol et al., 2012; van Daalen and Caswell, 2024). Like the variance, it is additively decomposable into between-group and within-group components as a measure of inequality (Rosenbluth, 1951).

The variance also has the advantage of being a central moment of the distribution, inviting connections to other moments. These moments highlight other aspects of the distribution (e.g., skewness), and at least some of them can also be decomposed into within- and between-group components. There seems to be no comparable linkage of absolute deviations to such other properties of the distribution.

The variance partitioning presented here is related to, but is not the same as, statistical analysis of variance (ANOVA). Our results do not provide hypothesis tests, because they lack an underlying sampling theory to obtain distributions of the variance components under a null hypothesis. If they did, most of the between-group components would be statistically insignificant. However, the factors investigated in demographic studies are usually known a priori to be of social or biological importance. The observation that the contribution of heterogeneity to variance is small does not imply that the factors are not worthy of attention; rather, it merely implies that stochasticity is itself an important factor and is worthy of study (see Caswell 2023, Section 8.5).

The mixing distribution plays a central role in variance partitioning. It is tempting to wonder if the results are somehow an artefact of that distribution. Insight into the potential effect of the mixing distribution can be obtained by generating marginal distributions at random, combining them to form a rank-one distribution, and then computing the variance ratio. Figure 2 shows the distribution of the variance ratio K obtained from a sample of 10,000 mixing distributions distributed uniformly over the simplex, for the case of the U.S. sex, race, and states (Table 2). This is a sample from the set of all possible rank-one mixing distributions (including many that are demographically absurd). Even so, the largest variance ratio in this sample is only 0.075. Since the variance ratio with a flat mixing distribution, is K=0.059, it seems unlikely from Figure 2 that this ratio could be increased very much by any other choice of mixing distribution.

A rank-one mixing distribution is required in order to partition the between-group variance into contributions from main effects and interactions. However, the calculation of the between-group variance itself, as in equation (19), places no restriction on the mixing distribution. It simply treats all the treatment combinations as levels of a single factor, reducing that portion of the calculation to a one-way design.

A recent survey of single-factor studies found that heterogeneity usually explained only 5%–10% of the variance in longevity (Caswell, 2023). The results of the present paper suggest that including multiple factors and their interactions may not increase K very much. Further comparative research is needed. The variance ratio for lifetime reproduction appears to be larger than that for longevity, and the variance components for LRO appear to be strongly influenced by conditions. van Daalen et al. (2022) found that heterogeneity in maternal age in a rotifer explained about 26% of the variance in LRO under laboratory conditions, but as little as 2% under most conditions that would lead to a stationary population. The sensitivity analysis of variance components (van Daalen and Caswell, 2020) may be helpful for exploring these relationships.

We encourage the use of the methods we present here to explore the contributions of multiple factors and their interactions to variance in demographic outcomes. Patterns await discovery.