Exact Statistical Distribution and Correlation of Human Height and Weight: Analysis and Experimental Confirmation

Mark P. Silverman

Department of Physics, Trinity College, Hartford, USA.

DOI: 10.4236/ojs.2022.125044                

Abstract

The statistical relationship between human height and weight is of especial importance to clinical medicine, epidemiology, and the biology of human development. Yet, after more than a century of anthropometric measurements and analyses, there has been no consensus on this relationship. The purpose of this article is to provide a definitive statistical distribution function from which all desired statistics (probabilities, moments, and correlation functions) can be determined. The statistical analysis reported in this article provides strong evidence that height and weight in a diverse population of healthy adults constitute correlated bivariate lognormal random variables. This conclusion is supported by a battery of independent tests comparing empirical values of 1) probability density patterns, 2) linear and higher order correlation coefficients, 3) statistical and hyperstatistics moments up to 6th order, and 4) distance correlation (dCor) values to corresponding theoretical quantities: 1) predicted by the lognormal distribution and 2) simulated by use of appropriate random number generators. Furthermore, calculation of the conditional expectation of weight, given height, yields a theoretical power law that specifies conditions under which body mass index (BMI) can be a valid proxy of obesity. The consistency of the empirical data from a large, diverse anthropometric survey partitioned by gender with the predictions of a correlated bivariate lognormal distribution was found to be so extensive and close as to suggest that this outcome is not coincidental or approximate, but may be a consequence of some underlying biophysical mechanism.

Keywords

Correlation of Height and Weight, Distribution of Height and Weight, Body Mass Index, Lognormal Distribution, Distance Correlation (dCor), Hyperstatistics

1. Introduction

Scientific interest in the values and correlations of anthropometric data trace back to the beginnings of modern statistics in the late 19th and early 20th Centuries with the researches of Quetelet, Galton, Pearson, and others [1] [2] [3]. These many studies established the Gaussian function as the mathematical expression best approximating the distribution of such human features as height, weight, and other biometric attributes. So pervasive has been the Gaussian distribution that it is ubiquitously referred to as the “normal distribution”, a reference probably dating back to Quetelet’s influential study of “the average man” (L’homme Moyen) in 1835 [4].

Although more refined studies have revealed that anthropometric data can show deviations from normality, attempts to find relationships between human height and weight remained uncertain, controversial and based on approximate or indirect methods such as data fitting [5], mechanical modeling [6], and gene identification [7]. Height and weight are of particular importance since they directly relate to the body mass index (BMI) [8], which is a measure of obesity and a risk factor for metabolic disease [9] and Alzheimer’s Disease [10]. A previous paper by Silverman and Lipscombe [11], to be referred to as Part I, determined the mathematically exact statistical distribution of BMI.

The present paper, to be regarded as Part II, provides evidence for the proposition that, in a healthy adult human population with access to adequate nutrition, height and weight are distributed as correlated bivariate lognormal random variables. This conclusion is supported by a comprehensive investigation comprising four independent components:

1) Tests of the correlation functions of the height and weight of a large anthropometric data set of individuals, partitioned by gender, against predictions of the bivariate lognormal distribution;

2) Search for nonlinear correlations, not attributable to the bivariate lognormal distribution, by means of a sensitive nonparametric algorithm known as distance correlation [12] [13];

3) Comparison of statistical tests of the empirical anthropometric data set with identical tests performed on comparably sized populations artificially created with correlated lognormal random number generators (RNGs);

4) Tests of the marginal distributions of height and weight against predictions of associated univariate lognormal distributions, and of the natural logarithms of height and weight against predictions of associated univariate normal distributions.

The outcome of this four-part analysis shows that linear and higher-order correlations of human height and weight are predictable in terms of the single Pearson correlation coefficient for height and weight employed in the lognormal probability density function (PDF). Moreover, agreement between the empirical data and the predictions of lognormal theory is so extensive as to suggest that the lognormal distribution of adult human height and weight is not approximate, but an exact distribution possibly characteristic of a more fundamental underlying biophysical mechanism.

1.1. Marginal Distributions of Height and Weight

Part I [11] reported the exact probability density function of BMI that follows mathematically from the defining relation

$B=W/H^2$ (1)

in which height H is expressed in meters (m), the corresponding weight W is expressed in the mass unit kilograms (kg), and B is the BMI expressed in kg/m². It is to be stressed that H, W, and therefore B, are random variables, which means that information and interpretations extrapolated from the BMI PDF refer to populations, and not to individuals, an essential point not always understood by the lay news media [11] [14].

The specific form that the general BMI density function takes depends on the statistical distributions of H and W. Such empirical distributions are often represented visually as histograms. However, if two random variables are not independent, then the histogram of each is a graphical representation of the marginal distribution of that variable, and provides no information regarding the correlation of the two variables. In Part I evidence was provided to show that height and weight of individuals measured in the Anthropometric Survey of U.S. Army Personnel (ANSUR)—a large data base comprising 4082 males and 1986 females [15] —were highly correlated. Figure 1 shows scatter plots of W against H for the separate male and female cohorts. The two patterns suggest a significant linear correlation. The superposed curves, to be discussed later, are the lines of regression (dashed red) and the conditional expectation functions (solid blue). Descriptive statistics are given in Table 1 for the two cohorts, together with theoretically predicted values, where appropriate. Details of Table 1 will be discussed at relevant points throughout the paper.

Figure 1. Correlation of weight and height for males (left) and females (right) of the ANSUR population. Lines of regression (dashed red) are obtained by the method of least squares. The conditional expectation functions of weight given height (solid blue) are calculated from the lognormal PDF (17), using ANSUR parameters in Table 1.

Table 1. Descriptive statistics of height and weight of ANSUR population.

It is to be recalled that a random variable X is lognormal if its natural logarithm, symbolized by $Y=\ln X$, is normal. As a matter of standard notation used in this paper, random variables are represented by upper case letters (e.g. X), and realizations of that variable (referred to as variates) are represented by lower case letters (e.g. x). Histograms of the natural logarithms of the ANSUR heights and weights, partitioned by gender, were shown in Part I to be satisfactorily described by PDFs of Gaussian form

$p_H^{\left(N\right)}\left(y\right)=\frac{1}{\sqrt{2\pi s_H^2}}{\text{e}}^{-\frac{{\left(y-m_H\right)}^2}{2s_H^2}}$ (2)

$p_W^{\left(N\right)}\left(y\right)=\frac{1}{\sqrt{2\pi s_W^2}}{\text{e}}^{-\frac{{\left(y-m_W\right)}^2}{2s_W^2}}$. (3)

A more detailed demonstration of the normality of lnH and lnW will be given in Section 6. From relations (2) and (3) follow the parent lognormal PDFs

$p_H^{\left(\Lambda \right)}\left(x\right)=\frac{1}{\sqrt{2\pi s_H^2}}\frac{{\text{e}}^{-\frac{{\left(\ln \left(x\right)-m_H\right)}^2}{2s_H^2}}}{x}$ (4)

$p_W^{\left(\Lambda \right)}\left(x\right)=\frac{1}{\sqrt{2\pi s_W^2}}\frac{{\text{e}}^{-\frac{{\left(\ln \left(x\right)-m_W\right)}^2}{2s_W^2}}}{x}$ (5)

with location parameters $\left(m_H,m_W\right)$ and scale parameters and $\left(s_H,s_W\right)$ for height and weight, respectively. Numerical values of these parameters are given in Table 1.

Superscripts N and Λ in the above PDFs signify normal and lognormal distributions, as well as symbolize the associated random variables (RVs)

$\begin{array}{l}X=\Lambda \left(m,s^2\right)\\ Y=N\left(m,s^2\right)\end{array}\}\Rightarrow \{\begin{array}{l}Y=\ln X\\ X={\text{e}}^Y\end{array}$ (6)

Note that parameters $\left(m,s^2\right)$ defining the random variables Y and X are the mean and variance of the normal variable Y. All statistics of the marginal distribution of X are predictable in terms of the parameters $\left(m,s^2\right)$ of Y [11] [16]. For example, the mean ${\mu }_X$, variance ${\sigma }_X^2$, skewness $S{k}_X$, and kurtosis $K_X$ of X take the forms [17]

${\mu }_X={\text{e}}^{m+\frac{1}{2}{s}^2}$. (7)

${\sigma }_X^2={\text{e}}^{2m}\left({\text{e}}^{2s^2}-{\text{e}}^{s^2}\right)$ (8)

$S{k}_X=\left({\text{e}}^{s^2}+2\right)\sqrt{{\text{e}}^{s^2}-1}$ (9)

$K_X={\text{e}}^{4s^2}+2{\text{e}}^{3s^2}+3{\text{e}}^{2s^2}-3$. (10)

1.2. Correlation of Height and Weight

When analyzing lognormal variates, it is often strategically easier—indeed necessary—to work with the logarithms of the variates, since these are distributed normally. Figure 2 shows scatter plots (black) of the scaled variates of $Y_W\equiv \ln W$ against scaled variates of $Y_H\equiv \ln H$ for males and females respectively in the ANSUR data set. The scaled variables $\left(U,V\right)$

$\begin{array}{l}U\equiv \left(\ln H-m_H\right)/s_H\\ V\equiv \left(\ln W-m_W\right)/s_W\end{array}$ (11)

with variates $\left(u,v\right)$ are measured with respect to their means and divided by their standard deviations, and are therefore dimensionless quantities distributed as standard normal variables of mean 0 and variance 1 if the variables H, W are lognormal. Figure 2 likewise clearly shows a strong linear correlation of U and V. The slope of the line of regression (dashed red) in each black scatter plot directly yields the corresponding Pearson correlation coefficient r [11] defined by

Figure 2. Correlation of scaled log weight and scaled log height for males (top left) and females (top right) of the ANSUR population. Lower panels show corresponding scatter plots created with correlated lognormal random number generators using the same empirical parameters. The patterns display a strong linear correlation. The Pearson correlation coefficient is equal to the slope of the associated lines of regression (dashed).

$r\equiv \langle UV\rangle =\frac{\langle \left(Y_H-m_H\right)\left(Y_W-m_W\right)\rangle }{s_H{s}_W}\equiv \frac{\mathrm{cov}\left(Y_H,Y_W\right)}{s_H{s}_W}$ (12)

where cov signifies covariance, as defined in Equation (12). Angular brackets are used in this paper to indicate expectation values. The scatter plots in red in Figure 2 were obtained by computer simulation using correlated lognormal RNGs, the details of which will be discussed in a later section. Suffice it to say at this point that the computer simulations employed the same distribution parameters that were extracted from the ANSUR height and weight data and are labeled $\left(m_H,m_W,s_H,s_W,r\right)$ for both male and female cohorts in Table 1. Lines of regression (dashed blue) to the simulated scatter plots are nearly identical to those of the empirical plots.

Once the correlation coefficient r has been determined empirically from the sample of normal variates $\left(u,v\right)$, the correlation coefficient $\rho$ of the parent lognormal variables $X_H\equiv {\text{e}}^{Y_H}=H$ and $X_W\equiv {\text{e}}^{Y_W}=W$ can be calculated theoretically from the relation [11]

${\rho }_{thy}=\frac{{\text{e}}^{r{s}_H{s}_W}-1}{\sqrt{\left({\text{e}}^{s_H^2}-1\right)\left({\text{e}}^{s_W^2}-1\right)}}$ (13)

and compared with the empirical correlation coefficient obtained directly from the data according to

${\rho }_{emp}\equiv \frac{\mathrm{cov}\left(X_H,X_W\right)}{{\sigma }_H{\sigma }_W}=\frac{\langle \left(H-{\mu }_H\right)\left(W-{\mu }_W\right)\rangle }{{\sigma }_H{\sigma }_W}$ (14)

in analogy to Equation (12). The implementation of Equation (14) can be achieved algebraically and geometrically:

1) Algebraic method: If $h_i$ and $w_i$ are variates of H and W, such as plotted in Figure 1, where $i=1,\cdots ,n$, then

${\rho }_{emp}=\frac{1}{n{\sigma }_H{\sigma }_W}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(h_i-{\mu }_H\right)\left(w_i-{\mu }_W\right)}$ (15)

in the limit of large n.

2) Geometric method: ${\rho }_{emp}$ is equal to the slope of the line of regression in a scatter plot of the scaled variables $\left(W-{\mu }_W\right)/{\sigma }_W$ against $\left(H-{\mu }_H\right)/{\sigma }_H$. To deduce ${\rho }_{emp}$ from the line of regression in the plot of the unscaled variables (H, W) in Figure 1, one multiplies the slope by the ratio of standard deviations ${\sigma }_H/{\sigma }_W$.

Five parameters $\left(m_1,m_2,s_1,s_2,r\right)$ are required to specify the PDF of two correlated bivariate normal RVs $\left(Y_1,Y_2\right)$

$\begin{array}{l}{p}_{Y_1,Y_2}^{\left(N,N\right)}\left(y_1,y_2\right)=\frac{1}{2\pi s_1{s}_2\sqrt{1-r^2}}{\text{e}}^{-q_Y/2}\\ q_Y=\frac{1}{1-r^2}\left[{\left(\frac{y_1-m_1}{s_1}\right)}^2-2r\left(\frac{y_1-m_1}{s_1}\right)\left(\frac{y_2-m_2}{s_2}\right)+{\left(\frac{y_2-m_2}{s_2}\right)}^2\right]\end{array}$ (16)

from which is derived the PDF of the parent bivariate lognormal RVs $\left(X_1,X_2\right)$ [11]

$\begin{array}{l}{p}_{X_1,X_2}^{\left(\Lambda ,\Lambda \right)}\left(x_1,x_2\right)=\frac{1}{2\pi s_1\text{}\text{}\text{}{s}_2\sqrt{1-r^2}}\frac{{\text{e}}^{-q_X/2}}{x_1{x}_2}\\ q_X=\frac{1}{1-r^2}[{\left(\frac{\ln \left(x_1\right)-m_1}{s_1}\right)}^2-2r\left(\frac{\ln \left(x_1\right)-m_1}{s_1}\right)\left(\frac{\ln \left(x_2\right)-m_2}{s_2}\right)\\ +{\left(\frac{\ln \left(x_2\right)-m_2}{s_2}\right)}^2]\end{array}$ (17)

Double superscripts N and Λ signify that both variables are normal in PDF (16) and lognormal in PDF (17). The expectation operations in Equation (12) and Equation (14) are performed respectively with PDF (16) and PDF (17).

If, as proposed in this paper, H and W are correlated bivariate lognormal variables, then all measurable statistical information concerning adult human height, weight, and their correlations, should be predictable from their joint distribution Equation (17) in terms of the five parameters (2 means, 2 variances, and 1 Pearson correlation) that define a given population. This statement has important implications for the study of obesity and its associated illnesses.

The BMI (1) was introduced by Quetelet in 1835 [18] and has been widely used up to present times by clinicians and epidemiologists as a proxy for obesity under the assumption that it correlates strongly with weight, but is independent of height. This assumption is itself predicated on a by-no-means obvious assumption that human weight in a healthy adult population varies as the square of an individual’s height. These assumptions will be examined later in this paper both empirically and theoretically. It is to be noted at this point, however, that both assumptions have elicited criticism, e.g. [19] [20] [21], leading to proposals of alternative power-law measures such as the Benn Index [22] and Rohrer’s Index [23], non-power law correlations such as [24] [25], and empirical parametric models such as [26], all purporting to determine more satisfactorily than BMI a single optimal relationship between human weight and height.

With regard to the goal of capturing the relationship between height and weight, the following general statistical principles must be emphasized. First, an exact PDF of the bivariate distribution of two correlated random variables provides all the statistical information that can be learned about these two correlated variables. And second, there is no single optimal mathematical expression—apart from the PDF and its equivalent transformations¹—that completely captures the statistical relation between two correlated random variables. Rather, the PDF provides a potentially infinite number of mathematical expressions that, together, characterize the complete relation between the two variables. From a practical standpoint, however, the number of testable expressions that can meaningfully characterize the correlation of two variables is limited by the size of the sample, since the intrinsic uncertainty increases with the order (i.e. power) of the variables, and can eventually exceed the mean value for a fixed sample size. These points will be elaborated on in the following sections.

1.3. Organization

The remainder of this paper is organized as follows.

In Section 2 the relation between weight (W) and height (H) is examined by means of the conditional expectation functions of W, given H.

In Section 3 the proposition that human height and weight are correlated lognormal variables is tested by examining generalized correlation functions of data sets $\left(H,W\right)$ and $\left(\ln H,\ln W\right)$.

In Section 4 the preceding data sets are each examined for nonlinear correlations beyond those attributable to the bivariate lognormal distribution by a procedure known as distance correlation.

In Section 5 the marginal distributions of $\left(H,W\right)$ and $\left(\ln H,\ln W\right)$ are tested against predictions of the univariate lognormal and normal distributions, respectively.

Section 6 examines the implications of the distribution of $\left(H,W\right)$ for the body mass index.

Section 7 discusses the computer simulation of correlated lognormal variables.

And last, the results of this comprehensive investigation are summarized and interpreted in Section 8.

2. Conditional Expectation of Weight, Given Height

The conditional expectation $\langle W^p|H\rangle$ of $W^p$ (for $p=1,2,\cdots$ ) given H, defined by the ratio

$W_p\left(h\right)\equiv \frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}{w}^p{p}_{H,W}^{\left(\Lambda ,\Lambda \right)}\left(h,w\right)\text{d}w}}{{\displaystyle \underset{0}{\overset{\infty }{\int }}{p}_{H,W}^{\left(\Lambda ,\Lambda \right)}\left(h,w\right)\text{d}w}}$, (18)

is a function $W_p\left(h\right)$ of the continuous variate of H. Since $W_p\left(h\right)$ derives from the joint PDF of W and H, it is more informative than an empirical line of regression such as obtained by the method of least squares or, more generally, the method of maximum likelihood [27].

Calculation of $W_p\left(h\right)$ in Equation (18) requires evaluation of two integrals whose kernel is the PDF $p_{H,W}^{\left(\Lambda ,\Lambda \right)}\left(h,w\right)$ given by Equation (17). The integrals can be greatly simplified by the transformation (11) to variables² $u=\left(\ln \left(h\right)-m_H\right)/s_H$ and $v=\left(\ln \left(w\right)-m_W\right)/s_W$, which re-expresses the bivariate normal PDF (16) more simply in the form

$f_{U,V}\left(u,v\right)=\frac{1}{2\pi \sqrt{1-r^2}}\exp \left(-\frac{1}{2\left(1-r^2\right)}\left(u^2-2ruv+v^2\right)\right)$ (19)

and, through the inverse transformation

$\begin{array}{l}h=\exp \left(m_H+s_{H}u\right)\\ w=\exp \left(m_W+s_{W}v\right)\end{array}$ (20)

leads to the conditional expectation

$W_p\left(u\right)\equiv \frac{{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left(\exp \left(s_{W}v+m_W\right)\right)}^{p}f\left(u,v\right)\text{d}v}}{{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f\left(u,v\right)\text{d}v}}$ (21)

as a function of u. Both integrals in (21) are readily evaluated in closed form., after which replacement of the normal variable u in terms of the lognormal variable h yields the general relation

$W_p\left(h\right)=h^{pr{s}_W/s_H}\exp \left(p{m}_W+\frac{1}{2}{p}^2{s}_W^2\left(1-r^2\right)-pr{m}_H{s}_W/s_H\right)$. (22)

The conditional expectation of the $p^{th}$ power of W is thus seen to have a power-law dependence on H with exponent $pr{s}_W/s_H$. The lowest two orders p = 1, 2 are of primary interest

$W_1\left(h\right)=h^{\frac{r{s}_W}{s_H}}\exp \left(m_W+\frac{1}{2}{s}_W^2\left(1-r^2\right)-\frac{r{m}_H{s}_W}{s_H}\right)$ (23)

$W_2\left(h\right)=h^{\frac{2r{s}_W}{s_H}}\exp \left(2m_W+2s_W^2\left(1-r^2\right)-\frac{2r{m}_H{s}_W}{s_H}\right)$ (24)

and yield the conditional variance and standard deviation

$\begin{array}{c}\mathrm{var}\left(W|H\right)=W_2\left(h\right)-{\left(W_1\left(h\right)\right)}^2\\ =h^{\frac{2r{s}_W}{s_H}}{\text{e}}^{\left(2m_W-\frac{2r{m}_H{s}_W}{s_H}\right)}\left({\text{e}}^{2s_W^2\left(1-r^2\right)}-{\text{e}}^{s_W^2\left(1-r^2\right)}\right)\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\end{array}$ (25)

$\begin{array}{c}{\sigma }_W\left(h\right)\equiv \sigma \left(W|H\right)=\sqrt{\mathrm{var}\left(W|H\right)}\\ =h^{\frac{r{s}_W}{s_H}}{\text{e}}^{\left(m_W-\frac{r{m}_H{s}_W}{s_H}\right)}{\left({\text{e}}^{2s_W^2\left(1-r^2\right)}-{\text{e}}^{s_W^2\left(1-r^2\right)}\right)}^{\frac{1}{2}}\end{array}$ (26)

Substituting in Equations (23) and (26) the bivariant lognormal parameters for each gender cohort of the ANSUR sample listed in Table 1 leads to the expressions

$\begin{array}{l}{\left(W_1\left(h\right)\pm {\sigma }_W\left(h\right)\right)}_{\text{M}}=\left(27.7071\pm 4.0621\right)h^{1.9987}\\ {\left(W_1\left(h\right)\pm {\sigma }_W\left(h\right)\right)}_{\text{F}}=\left(23.2156\pm 3.1523\right)h^{2.1923}\end{array}$ (27)

where subscripts M and F signify male and female, respectively.

Plots of the conditional expectations $W_1\left(h\right)$ in Equation (27) comprise the solid blue curves superposed on the scatter plots in Figure 1. Although $W_1\left(h\right)$ is a power law and the line of regression (dashed red) is linear, the two curves are virtually indistinguishable over the densest part of the plots. However, it is important to bear in mind the conceptual difference between the two curves: the line of regression is merely a fit to data, whereas the mathematical relation (22), the specific exponent $pr{s}_W/s_H$, and the numerical values in expressions (27) are predictions drawn from the lognormal PDF. Figure 3 graphically displays the full information content of relation (27) by displaying the regions of ±1 standard deviation about the means for the two cohorts.

The numerical values of the exponents in relations (27) bear out the fundamental assumption underlying the use of BMI that weight is a quadratic function of height in a healthy adult population. Nevertheless, for a different set of values of the parameters $\left(s_H,s_W,r\right)$, such as may characterize a demographic different from the one represented by the ANSUR population, the lognormal predicted exponent could be different.

3. Tests of Correlation Functions of Height (H) and Weight (W)

Correlation functions of order (p, q), defined as follows

$C_{p,q}\left(s_H,s_W,r\right)\equiv \langle {\left(\frac{H-{\mu }_H}{{\sigma }_H}\right)}^p{\left(\frac{W-{\mu }_W}{{\sigma }_W}\right)}^q\rangle$ (28)

$R_{p,q}\left(r\right)\equiv \langle {\left(\frac{\ln \left(H\right)-m_H}{s_H}\right)}^p{\left(\frac{\ln \left(W\right)-m_W}{s_W}\right)}^q\rangle$, (29)

Figure 3. Plots of the conditional expectation $W_1\left(h\right)$ for male (solid blue line) and female (solid red line) of the ANSUR populations centered on regions (blue for males, red for females) of ±1 standard deviation ${\sigma }_W\left(h\right)$. The region of overlap appears purple.

generalize the standard covariance ( $p=q=1$ ) expressed in relations (12) and (14). Expectations (28) and (29) are to be implemented with the bivariate lognormal PDF (17) where indices (p, q) independently take on integer values (1, 2, …).

As a matter of notation and terminology, evaluation of functions $C_{p,q}$ and $R_{p,q}$ by substitution of empirical parameters for the arguments yield the numerical correlation coefficients $c_{p,q}$, $r_{p,q}$, where $r_{1,1}\equiv r$ and $c_{1,1}\equiv \rho$ as conventionally defined. If the proposition that $\left(H,W\right)$ are bivariate lognormal variables is valid, then $c_{p,q}$ and $r_{p,q}$ should be predictable from the arguments shown in Equations (28) and (29) and the parameters $\left(s_H,s_W,r\right)$ listed by gender in Table 1. It is to be recalled that the first two parameters (standard deviations) were obtained empirically from the variates of the marginal distributions $\ln H$ and $\ln W$ of the ANSUR populations, whereas the third parameter (Pearson correlation coefficient) was obtained empirically from the joint distribution of $\ln H$ and $\ln W$, such as exhibited in Figure 2. The associated sets of means $\left(m_H,m_W\right)$, which are also listed by gender in Table 1, drop out of relations (28) and (29) by virtue of their expressions as ratios. Further details are given in Part I [11]. To facilitate reading the tables to follow, indices of the correlation functions and coefficients will be expressed as arguments, e.g. $R\left(p,q\right)$ and $r\left(p,q\right)$ in the tables.

Calculation of the correlation functions $C_{p,q}$ and $R_{p,q}$ requires evaluation of the expectation values

$C_{p,q}={\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\left(\frac{h-{\mu }_H}{{\sigma }_H}\right)}^p}}{\left(\frac{w-{\mu }_W}{{\sigma }_W}\right)}^q{p}_{\left(H,W\right)}^{\left(\Lambda ,\Lambda \right)}\left(h,w\right)\text{d}h\text{d}w$ (30)

$R_{p,q}={\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\left(\frac{\ln \left(h\right)-m_H}{s_H}\right)}^p}}{\left(\frac{\ln \left(w\right)-m_W}{s_W}\right)}^q{p}_{\left(H,W\right)}^{\left(\Lambda ,\Lambda \right)}\left(h,w\right)\text{d}h\text{d}w$. (31)

As in the previous section, the two integrals can be greatly simplified by the transformation (11) to variables $u=\left(\ln \left(h\right)-m_H\right)/s_H$ and $v=\left(\ln \left(w\right)-m_W\right)/s_W$, which results in the bivariate normal PDF (19). Substitution for the means $\left({\mu }_H,{\mu }_W\right)$ and standard deviations $\left({\sigma }_H,{\sigma }_W\right)$ by use of Equations (7) and (8) then leads to the operational expressions

$C_{p,q}={\left({\text{e}}^{s_H^2}-1\right)}^{-p/2}{\left({\text{e}}^{s_W^2}-1\right)}^{-q/2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left({\text{e}}^{s_{H}u-\frac{1}{2}{s}_H^2}-1\right)}^p}}{\left({\text{e}}^{s_{W}v-\frac{1}{2}{s}_W^2}-1\right)}^q{f}_{U,V}\left(u,v\right)\text{d}u\text{d}v$ (32)

and

$R_{p,q}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{u}^p{v}^q{f}_{U,V}\left(u,v\right)\text{d}u\text{d}v}}$. (33)

In the symmetric case ( $p=q$ ), the two correlation functions will simply be designated $C_p$ and $R_p$.

The integral in relation (33) displays a number of symmetries: 1) $R_{p,q}=R_{q,p}$, 2) $R_{p,q}=0$ if $p+q$ is odd, otherwise 3) if $p+q$ is even, then $R_{p,q}$ is of order $r^d$ where d is the smaller of p and q. Symmetries (2) and (3) do not necessarily hold for $C_{p,q}$. As a matter of terminology in the sections to follow, $R_{p,q}$ and $C_{p,q}$ will be referred to as “asymmetric odd” if $p+q$ is odd and “asymmetric even” if $p+q$ is even.

3.1. Calculation and Measurement of Correlation Functions $R_{p,q}$

Because the correlation of U and V is determined entirely by the parameter r in the probability density $f\left(u,v\right)$, it is useful to examine the structure of $f\left(u,v\right)$ graphically. The left panels of Figure 4 display a sequence of theoretical 3-dimensional density plots for correlation coefficients r = 0.0, 0.5, and 0.95, which span nearly the entire positive range of r. The maximum r = 0.95 was chosen, rather than r = 1, because the latter value is simply a straight line coincident with the diagonal axis. The right panels view the underside of the density plots, or, equivalently, the projection of the plots onto the $\left(u,v\right)$ plane. In the top panels, vectors U and V are independent (r = 0), and $f_{U,V}\left(u,v\right)$ factors into a product $f_U\left(u\right)f_V\left(v\right)$. Since the exponential in $f_{U,V}\left(u,v\right)$ for r = 0 is the sum $u^2+v^2$, a transformation of variables converts that sum into the square of a radial variable, which accounts for the circular symmetry of the top right panel. As the correlation coefficient increases toward +1, the density function profile becomes increasingly linear along the diagonal for which $uv>0$, as clearly shown in the bottom panels. For $r\to -1$, the profile (not shown) would approach linearity along the diagonal for which $uv<0$.

Whereas the left side of Figure 4 shows the actual probability density profiles, the images on the right side show smooth shapes to which scatter plots of a sample of discrete paired variates $\left(u_i,v_i\right)$, $i=1,\cdots ,n$, approach as $n\to \infty$. This is borne out by Figure 5, which shows simulated plots of n = 2000 pairs of correlated standard normal variates of increasing correlation coefficient r. In each panel, the quantity “rho” designates the correlation parameter supplied to the RNGs; the quantity “slope” is the slope of the line of regression (dashed blue), which equals the actual correlation coefficient of the simulated sample. The two numbers in each scatter plot are close, but not identical because the scatter plot comprises a finite random sample. The connection between the PDF $f_{U,V}\left(u,v\right)$ and the scatterplot of $\left(U,V\right)$ for each value of r is particularly clear when one compares the right side panels of Figure 4 to the corresponding plots of Figure 5. The orientations of the two sets of figures may be different, but the distributions are invariant to orientation.

Correlation function $R_{p,q}\left(r\right)$ expressed in Equation (33) can be evaluated in closed form, and depends only on the Pearson coefficient r. Table 2 lists the first six orders of the symmetric correlation functions and the most pertinent of the asymmetric correlation functions. It is seen that beyond the basic covariance (12), the higher-order symmetric correlations are increasingly nonlinear in r. Asymmetric correlation functions of the form $R_{p,1}$ in Table 2, where p is the exponent of the scaled variable for $\ln H$, address the issue (referred to in the

Figure 4. Left Panels: Probability density profiles $f\left(u,v\right)$ of correlated standard normal variates $\left(u,v\right)$ for correlation coefficients r = 0 (top), 0.5 (center), 0.95 (bottom). Right Panels: Profiles as viewed from the underside of the associated density plots. The patterns are the smooth shapes of scatter plots of discrete samples in the limit of infinite sample size.

Figure 5. Simulated scatter plots of 2000 correlated pairs of standard normal variates with correlation coefficients increasing from 0 to nearly 1. In each plot, “rho” is the numerical value of the correlation parameter supplied to the RNGs, and “slope” (of the line of regression) is the actual correlation coefficient produced by the simulation.

Table 2. Correlation functions R (p, q) of powers of standard normal variables: U^p and V^q.

Introduction) concerning how weight correlates with powers of height. Also listed are asymmetric correlation functions characterizing how powers of weight correlate with powers of height.

Although the Pearson correlation coefficient r is the expectation value of a composite random variable UV, it is regarded here as a nondistributed quantity since the values of r for male and female cohorts used throughout this paper are fixed parameters extracted from the ANSUR data. The same is true of the means $\left(m_H,m_W\right)$ and variances $\left(s_H^2,s_W^2\right)$ for the male and female cohorts. The variance of the correlation coefficient $r_{p,q}$ therefore characterizes only the variation of the product $U^p{V}^q$ in the defining integral and can be expressed by [28]

$\mathrm{var}\left(r_{p,q}\right)=r_{2p,2q}-r_{p,q}^2$. (34)

It then follows that the standard error (se) of $r_{p,q}$ takes the form

$se\left(r_{p,q}\right)=\sqrt{\frac{\mathrm{var}\left(r_{p,q}\right)}{n}}=\sqrt{\frac{r_{2p,2q}-r_{p,q}^2}{n}}$ (35)

where n is the sample size. By the same reasoning, the standard error of the correlation coefficient $c_{p,q}$ is

$se\left(c_{p,q}\right)=\sqrt{\frac{\mathrm{var}\left(c_{p,q}\right)}{n}}=\sqrt{\frac{c_{2p,2q}-c_{p,q}^2}{n}}$ (36)

since $c_{p,q}$ likewise depends on the fixed ANSUR parameters. In general, however, the distribution function and moments of even the lowest correlation coefficient $r_{1,1}$ are difficult to obtain in closed form [29], and, at the time of writing, the author knows of no calculation in the literature of the exact closed-form distribution functions and moments of higher-order correlation coefficients of two normal or two lognormal variables.

Theoretical and empirical ANSUR correlation coefficients $r_{p,q}$ are displayed in Table 3 for symmetric indices up to p = 6 and for a selection of asymmetric indices ranging from (2, 1) to (6, 4). For other correlation orders, the corresponding standard errors were too large relative to the means for a comparison of theory and experiment to be meaningful. Examination of Table 3 shows striking agreement between lognormal theory and the ANSUR data.

For all but two correlation coefficients of the female cohort listed in the table, the magnitude of the difference between theoretical (thy) and empirical (emp) coefficients did not exceed 1 standard error. In other words, assuming, as justified by the Central Limit Theorem [30], that the relative error

$z=\left({r_{p,q}|}_{emp}-{r_{p,q}|}_{thy}\right)/se\left(r\right)$ follows a normal distribution, then rejection of the hypothesis ${r_{p,q}|}_{emp}={r_{p,q}|}_{thy}$ at the conventional 5% threshold would require

$\left|z\right|\ge 1.65$ [31], which was not the case for any of the correlation coefficients of the female cohort. The relative errors of the two exceptional coefficients $r_{3,2}$ and $r_{4,3}$ were approximately 1.074 and 1.243, respectively.

Table 3. Correlation coefficients r (p, q) of (Scaled lnH)^pand (Scaled lnW)^q.

Relative errors of the male cohort were overall larger than those of the female cohort, although most did not exceed the 5% threshhold value of 1.65. Exceptions occurred primarily among the higher orders of the asymmetric odd coefficients whose theoretical means were zero and standard errors large. Under such circumstances, large deviations are to be expected and would require a larger sample size for resolution. Moreover, expressions (35) and (36) give lower limits for the standard errors since, in accordance with stated assumptions, they do not take account of the variation in lognormal parameters. It is therefore likely that a more exact estimate of the relative errors would be lower than those listed in Table 3.

In summary, to appreciate how extensive and close is the agreement of the theoretical and empirical correlations displayed in Table 3, one must bear in mind the following context. Theoretical predictions of $r_{p,q}$ were obtained from the bivariate normal distribution of $\ln H$ and $\ln W$ ; empirical values of $r_{p,q}$ were obtained from the natural logarithms of the raw ANSUR sample. If adult heights and weights were not distributed lognormally, then the comprehensive correspondence by pure chance, especially in the female cohort, of these two sets of numbers would be extremely improbable. For example, if either or both of the attributes of height and weight were themselves normally distributed, as had long been assumed, the logarithm of the variates would depart significantly from a normal distribution, as was demonstrated in Part I [11].

Nevertheless, the results in Table 3 raise a curious question. Why does the female cohort appear to bear out the predictions of lognormal theory more closely than the male cohort despite the fact that the number of men sampled is about twice that of women? As discussed in Part I [11], the Anthropometric Survey of U.S. Army Personnel was undertaken to obtain data representative of the “Total Army” [15] with regard to making accurate decisions concerning clothing, protective equipment, workspaces, and other size-dependent, work-related matters. The survey measured more than 90 human attributes directly and compiled data demographically in terms of race, ethnicity, gender, age, and geographic location. For the analyses in this paper and in Part I, the data were partitioned by gender only. Therefore, both the male and female cohorts can be regarded as diverse samples of fundamentally healthy adults. However, since there are considerably fewer women in the U.S. Army than men, it is conceivable that, irrespective of other demographic characteristics, the women who joined the U. S. Army and took part in the ANSUR sample formed a more homogeneous group in regard to body type and physical fitness than the men. Such an explanation would seem likely, since there is no biophysical basis to believe that male height and weight would be statistically distributed by a probability function of different mathematical form than female height and weight.

3.2. Density Plots Associated with Correlation Functions $R_{p,q}$

Whereas the correlation coefficients $r_{p,q}$ are single numbers quantifying the correlation of $U^p$ and $V^q$ (i.e. powers of the scaled variables of $\ln H$ and $\ln W$ ) for a specified sample, the actual scatter plots of the two sets of variates yield a more comprehensive visual perspective of their correlation. Figure 6 shows such plots for symmetric correlation coefficients of orders p = 2, 3, 4, 5 of the ANSUR male cohort (sample size 4082). The patterns for the female cohort are similar, although less dense (sample size 1986) and not shown.

The correlations expressed in Figure 6 are highly nonlinear in two ways. Geometrically, the overall patterns of order $p>1$ do not show a well-defined linear variation such as seen in Figure 2 for $p=1$. And algebraically, the associated theoretical correlation function is of order $r^p$ in the Pearson correlation parameter, as summarized in Table 2. The four plots in Figure 6 illustrate the characteristic property that, with increasing order p, the density of points clusters

Figure 6. Empirical scatter plots of V^p against U^p for p = 2 (top left), 3 (top right), 4 (bottom left), 5 (bottom right), compiled from the ANSUR data. Patterns show a highly nonlinear correlation of weight and height. U is the scaled variable for lnH; V is the scaled variable for lnW.

more tightly about the coordinate axes. Fluctuations for even p extend primarily into the first quadrant (since both variates are positive). Fluctuations for odd p extend primarily into the first and third quadrants for $r>0$ (and into the second and fourth quadrants for $r<0$ ; not shown).

The empirical patterns and properties in Figure 6 are reproduced nearly identically (apart from random fluctuations) in the computer simulated patterns shown in Figure 7. The simulations were created by means of correlated lognormal RNGs using the ANSUR lognormal parameters in Table 1. The reproduction of empirical correlation scatter plots (and correlation coefficients) by computer simulation extends as well to asymmetric even and odd correlation functions, as displayed in Figure 8 for the pairs of variables $\left(U^2,V\right)$, $\left(U^3,V\right)$, and $\left(U^3,V^2\right)$. Plots in black are empirical; those in red are simulated for a population of corresponding size. Altogether, these results support the conclusion that the nonlinear correlations of height and weight stem exclusively from the properties of the lognormal distribution function and depend on no correlation parameters other than the Pearson coefficient r.

Probability density functions for the correlated powers $\left(U^p,V^q\right)$, which yield

Figure 7. Simulated scatter plots of V^p against U^p for p = 2 (top left), 3 (top right), 4 (bottom left), 5 (bottom right), obtained from correlated lognormal RNGs using the lognormal parameters in Table 1. The shapes of the patterns and extent of fluctuations closely resemble the corresponding empirical plots in Figure 6.

Figure 8. Empirical (black) and simulated (red) scatter plots of V^q against U^p for the asymmetric orders (p, q) = (2, 1) (top), (3, 1) (middle), (3, 2) (bottom).

the patterns approached by scatter plots in Figures 6-8 in the limit of infinite sample size, can be constructed by means of the Dirac delta function as follows

$f_{X,Y}^{\left(p,q\right)}\left(x,y\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{f}_{U,V}\left(u,v\right)\delta \left(u^p-x\right)\delta \left(v^q-y\right)\text{d}u\text{d}v}}$. (37)

The subscripts $\left(X,Y\right)$ represent the random variables whose lower-case variates are respectively $x=u^p$, $y=v^q$. Powers of standard normal variables like U, V do not, in general, follow known, named distributions to which variables X and Y can be assigned [32]. A brief recapitulation of the properties and identities of the Dirac delta function is given in Part I [11] and in mathematical physics books [33]. The integral (37) can be evaluated in closed form, but gives rise to long, cumbersome expressions for $p>2$, which will not be reproduced here.

Plots of probability density (37) for the symmetric cases p = 2 and p = 3 are shown respectively in Figure 9 and Figure 10. Left-side panels show the density

Figure 9. Left panel: Plot of the probability density $f_{X,Y}^{\left(2\right)}\left(u^2,v^2\right)$ for r = 0.5. Right panel: View from the underside highlights the profile to which the scatter plot of V² against U² in the top left of Figure 6 approches in the limit of infinite sample size.

Figure 10. Left panel: Plot of the probability density $f_{X,Y}^{\left(3\right)}\left(u^3,v^3\right)$ for r = 0.5. Right panel: View from the underside shows the profile approached by the scatter plot of V³ against U³ in the top right of Figure 6 in the limit of infinite sample size.

patterns from above the $\left(u^p,v^p\right)$ plane; right-side panels give complementary images from below. The function $f_{X,Y}^{\left(2\right)}\left(u^2,v^2\right)$ in Figure 9 shows a concentration of probability along the vertical axis with density decreasing with distance into the first quadrant, as shown empirically in Figure 7 (top left). Function $f_{X,Y}^{\left(3\right)}\left(u^3,v^3\right)$ in Figure 10 (left side) sharply delineates the “cross” of probability density along the coordinate axes, whereas the projection of the pattern onto the $\left(u^3,v^3\right)$ plane (right side) captures the point distribution in the corresponding plot of Figure 7 (top right).

3.3. Calculation and Measurement of Correlation Functions $C_{p,q}$

$C_{p,q}\left(s_H,s_W,r\right)$ in Equation (32) expresses the correlation of height and weight directly, rather than of their logarithms. The integral in Equation (32) can be evaluated in closed form of which the lowest orders pertinent to this paper are given in Table 4 for two arbitrary, but correlated, lognormal variables $X_1$ and $X_2$ with bivariate parameter set $\left(m_1,m_2,s_1,s_2,r\right)$. Beyond the symmetric order $p=4$ and asymmetric order $\left(p,q\right)=\left(4,2\right)$ expressions for the functions become overly long and not especially informative. Some points to note: 1) symmetric functions $C_p\left(s_H,s_W,r\right)$ are invariant under the interchange of parameters $s_H$ and $s_W$ ; 2) asymmetric odd functions $C_{p,q}$ do not identically vanish, as do asymmetric odd $R_{p,q}$ ; 3) none of the functions $C_{p,q}$ vanishes for r = 0, as do the functions $R_{p,q}$ (except those of the form $R\left(2m,2n\right)$ where m, n are integers > 0).

Substitution of the ANSUR parameters of Table 1 for male and female cohorts into the functions of Table 4 and variance of Equation (36) yield the corresponding empirical correlation coefficients summarized in Table 5. Agreement of the empirical values with lognormal theory is again excellent, apart from the highest orders where the standard errors are large relative to the means and signify that a larger sample size is required.

4. Test for Nonlinear Correlations by the Method of Distance Correlation

If adult human height and weight are bivariate lognormal variables, then all measures of their correlation must be calculable from the PDF (17). Evidence for the proposition of bivariate lognormality has been supported up to this point by the agreement of measured and lognormally predicted correlation coefficients and comparison of empirical and lognormally simulated probability density plots. The question remains, however, as to whether there may be nonlinear correlations beyond those intrinsic to the bivariate lognormal distribution. Correlation of distances [12], provides a sensitive method for testing the independence of random vectors.

As initially presented by its developers, the term distance covariance (dCov) of two random vectors X and Y, defined by

$V\left(X,Y\right)\equiv \Vert g_{X,Y}\left(t,s\right)-g_X\left(t\right)g_Y\left(s\right)\Vert$, (38)

Table 4. Correlation functions C (p, q) of Lognormal variables $X_1^p$ and $X_2^q$ (Notation: $E\left(x\right)\equiv \exp \left(x\right)$ ).

Table 5. Correlation coefficients c (p, q) of (Scaled H)^p and (Scaled W)^q.

is a measure of the difference between the joint characteristic function (CF) [28]

$g_{X,Y}\left(t,s\right)={\displaystyle \iint f_{X,Y}\left(x,y\right)\exp \left(i\left(tx+sy\right)\right)\text{d}x\text{d}y}$ (39)

and the product of the marginal CFs

$\begin{array}{l}{g}_X\left(t\right)={\displaystyle \int f_X\left(x\right)\exp \left(itx\right)\text{d}x}\\ g_Y\left(s\right)={\displaystyle \int f_Y\left(x\right)\exp \left(isy\right)\text{d}y}\end{array}$ (40)

As indicated in Equations (39) and (40), the CF $g_Z$ is the Fourier transform of the corresponding probability density $f_Z$ of some specific random variable or set of random variables Z .

The actual evaluation of $V\left(X,Y\right)$ in Equation (38), together with its properties and associated theorems, is given in Ref. [12]. For the purposes of this paper, sufice it to say that $V\left(X,Y\right)$ involves a weighted integral of ${\left|g_{X,Y}\left(t,s\right)-g_X\left(t\right)g_Y\left(s\right)\right|}^2$ over the Fourier coordinates t, s. The associated quantity of distance variance (dVar), given by $V\left(X,X\right)$, would then be the same weighted integral over ${\left|g_{X,X}\left(t,s\right)-g_X\left(t\right)g_X\left(s\right)\right|}^2$ with analogous expressions for $V\left(Y,Y\right)$. The distance correlation (dCor), expressed by $R\left(X,Y\right)$, is then defined in terms of dCov and dVar by the relation

$R\left(X,Y\right)\equiv \frac{\left|V\left(X,Y\right)\right|}{\sqrt{\left|V\left(X,X\right)\right|\left|V\left(Y,Y\right)\right|}}$. (41)

Equation (41) resembles in form Equation (12) or Equation (14) for the Pearson correlation coefficient, but its properties are significantly different, as well as its empirical evaluation. The author is unaware of any theoretical derivations of the probability density function or statistical moments of $R\left(X,Y\right)$. However, if X and Y are standard normal variables, then $R\left(X,Y\right)$ has been evaluated in the following closed form [12]

$R^{\left(N,N\right)}\left(X,Y\right)=\frac{\rho \arcsin \left(\rho \right)+\sqrt{1-{\rho }^2}-\rho \arcsin \left(\rho /2\right)-\sqrt{4-{\rho }^2}+1}{1+\left(\pi /3\right)-\sqrt{3}}$ (42)

where the superscript (N, N) signifies the special case of bivariate standard normality, and $\rho$ is the associated Pearson correlation coefficient.

4.1. Statistical Application of Distance Correlation (dCor) to Height and Weight

The procedure for applying dCor to a statistical system is as follows: Given samples $\left(x_i,y_i\right)$ for $i=1,\cdots ,n$ of two random variables X and Y, construct the statistic

$A_{k,l}=a_{k,l}-a_{k·}-a_{·l}+a$, (43)

where $\left(k,l=1,\cdots ,n\right)$ and

$\{\begin{array}{l}{a}_{k,l}=\left|x_k-x_l\right|,a_{k·}=\frac{1}{n}{\displaystyle \underset{l=1}{\overset{n}{\sum }}{a}_{k,l}}\\ a_{·l}=\frac{1}{n}{\displaystyle \underset{k=1}{\overset{n}{\sum }}{a}_{k,l}},\text{}\text{}a=\frac{1}{n^2}{\displaystyle \underset{k,l=1}{\overset{n}{\sum }}{a}_{k,l}}\end{array}$ (44)

and the associated statistic

$B_{k,l}=b_{k,l}-b_{k·}-b_{·l}+b$ (45)

where

$\{\begin{array}{l}{b}_{k,l}=\left|y_k-y_l\right|,b_{k·}=\frac{1}{n}{\displaystyle \underset{l=1}{\overset{n}{\sum }}{b}_{k,l}}\\ b_{·l}=\frac{1}{n}{\displaystyle \underset{k=1}{\overset{n}{\sum }}{b}_{k,l}},\text{}\text{}\text{}b=\frac{1}{n^2}{\displaystyle \underset{k,l=1}{\overset{n}{\sum }}{b}_{k,l}}\end{array}$ (46)

The squares of the empirical dCov and dVar are then given by

$V_n^2\left(X,Y\right)=\frac{1}{n^2}{\displaystyle \underset{k,l=1}{\overset{n}{\sum }}{A}_{k,l}{B}_{k,l}}$ (47)

$\begin{array}{l}{V}_n^2\left(X,X\right)=\frac{1}{n^2}{\displaystyle \underset{k,l=1}{\overset{n}{\sum }}{A}_{k,l}^2}\\ V_n^2\left(Y,Y\right)=\frac{1}{n^2}{\displaystyle \underset{k,l=1}{\overset{n}{\sum }}{B}_{k,l}^2}\end{array}$ (48)

from which follows the square of the empirical dCor

$R_n^2\left(X,Y\right)=\frac{V{\left(X,Y\right)}^2}{\sqrt{V_n^2\left(X,X\right)V_n^2\left(Y,Y\right)}}$ (49)

corresponding to Equation (41). The statistic $R_n\left(X,Y\right)$ 1) ranges between 0 and 1, 2) is 0 only if X and Y are independent, and 3) approaches the theoretical dCor $R\left(X,Y\right)$ in the limit of infinite sample size n [12].

4.2. Distance Correlation Test of the ANSUR Data

In the sequence of analyses of the ANSUR data to follow, the pair of variables $\left(X,Y\right)$ was taken to be the scaled sets 1) $\left(H,W\right)$, 2) $\left(\ln H,\ln W\right)$, and 3) $\left(\ln H,{\left(\ln W\right)}_{\text{red}}\right)$. It is to be recalled that a “scaled” variable is in dimensionless form of zero mean and unit variance. The subscript “red” in set (3) indicates that the scaled variates of $\ln W$ were reduced by subtraction of the regression of $\ln W$ on $\ln H$. The reason for this reduction and the way it was implemented will be clarified shortly.

The algorithm (Equations (43) to (48)) leading to Equation (49) is straightforward to implement by computer. However, for sample sizes on the order of thousands, the computation time is impractically long. To circumvent this difficulty, a sampling procedure analogous to bootstrapping [34] [35] was employed.

Each bivariate pair $\left(h_i,w_i\right)$ of height and weight in the ANSUR data set is labeled by an index i , referred to here as the “participant number”, that ranges from 1 to the full sample size n, which is $n_m=4082$ for the male cohort and $n_f=1986$ for the female cohort. Participants in the survey were apparently measured and recorded in random order, as shown in Figure 11, which displays scatter plots by gender of the scaled height and weight vs participant number. Although histograms of the variates, analyzed in part I [11], are precisely matched by lognormal distributions, the point density plots in Figure 11 are well represented by uniform distributions across the entire range of participants. In other words, each vertical slice of points of sufficient width contains a statistical spread of variates equivalent to any other vertical slice of the same width. Given this uniform density, the ranges $n_m$ and $n_f$ were respectively partitioned into 20 and 10 subgroups of 200 participants each, as shown in Figure 12. Distance correlation of height and weight was then evaluated for 50 consecutive participants in each subgroup, starting at the participant numbers marked by diamond plotting symbols in the figure.

For example, $dCo{r}_1$ was calculated from participants [200 - 249], $dCo{r}_2$ from partipants [400 - 449], and so on up to $dCo{r}_{20}$ from participants [4000 to 4049] for males and up to $dCo{r}_9$ from participants [1800 to 1849] for females. As with standard bootstrapping, this method of calculating dCor by repeated sampling not only circumvented what otherwise would have been an excessively long computer calculation, but it also provided a vector of dCor values from which to estimate dCor uncertainty in the absence of a known statistical distribution.

The results of the analyses of distance correlation are summarized in Table 6.

Section I of the table records the distance correlation of the scaled variables H and W. Particularly striking is the close agreement between the empirical values

Figure 11. Plots of height (left panels) and weight (right panels) of individual participants in the ANSUR sample: male cohort (top panels); female cohort (bottom panels).

(a) (b)

Figure 12. Illustration of the resampling strategy for calculation of distance correlation of participants’ height and weight in the ANSUR population. Diamond plotting symbols mark participant numbers which begin each subgroup of 50 participants to be sampled over the range (a) 200 to 4000 (male cohort), (b) 200 to 1850 (female cohort).

obtained from the ANSUR sample (columns 2 and 3 for male and female cohorts, respectively) and values created by computer simulation using a bivariate lognormal RNG (columns 4 and 5 for the closely corresponding sample sizes 4000 and 2000, respectively). Numerical values designated by “rho” at the top of columns 4 and 5 are the correlation parameters supplied to the RNG. These values

Table 6. Test of nonlinear correlations by means of distance correlation: $\text{dCor}\left(X,Y\right)$.

(from Table 1) correspond to the empirical Pearson correlation coefficients of the variables $\left(\ln H,\ln W\right)$. The empirical Pearson correlation coefficients $\rho$ (columns 2 and 3) of the variables $\left(H,W\right)$ are closely matched by the corresponding values (columns 4 and 5) obtained by computer simulation. In short, the dCor values of Section I are consistent with attributing the entire correlation of height and weight, including any nonlinear contributions, to the bivariate lognormal distribution.

Section II of the table records distance correlation of the scaled variables $\left(\ln H,\ln W\right)$. Computer simulated populations of sizes approximating the male and female ANSUR cohorts were generated by use of bivariate normal RNGs. Agreement between empirical and corresponding computer simulated values is again very close, especially for the female cohort. The outcome indicates that the full correlation of $\ln H$ and $\ln W$ is attributable to the bivariate normal distribution which, in turn, derives from the parent lognormal distribution function.

The distance correlation procedure tests random vectors for independence irrespective of their specific distributions, provided the first moments are finite [13]. It is a nonparametric test that can reveal correlations even when the Pearson correlation coefficient is null. However, if the Pearson correlation of two normal vectors is null, then those vectors are fully independent—i.e. there is no latent nonlinear correlation. Section III of the table exploits this point to ascertain whether the correlation between height and weight exhibited in Figure 1 and Figure 2 have a nonlinear contribution not attributable to the parent lognormal or derived normal distributions.

The basic idea is to subtract from the ordinate of each point in the scatter plots in Figure 2 the corresponding ordinate of the line of regression. For correlations of scaled standard normal variables $\left(X,Y\right)$, the line of regression takes the simple form $y=rx$, where r is the slope of the line and is equal to the Pearson correlation coefficient. A scatter plot is then made of the reduced scaled log weight

${\left(\ln W\right)}_{red}\equiv \left(\ln W-m_W\right)/s_W-r\left(\ln H-m_H\right)/s_H$ (50)

against the scaled log height $\left(\ln H-m_H\right)/s_H$, as shown in Figure 13. The scatter plots (black points) for male (left) and female (right) ANSUR cohorts exhibit the same isotropic patterns (apart from fluctuations) as the simulated scatter plot for null correlation (rho = 0) in Figure 5. Quantitatively, the slopes of the lines of regression (dashed red) of the reduced plots are respectively $-4.1118\times {10}^{-11}$ (male cohort) and $-4.1011\times {10}^{-11}$ (female cohort). In other words, removal of the lines of regression from the empirical scatter plots of

Figure 13. Scatter pattern (black points) and associated line of regression (dashed red) of zero slope signifying a null Pearson correlation of log weight and log height when the log weight variates were reduced by corresponding values of the line of regression (dashed blue) of the original scatter plot (cyan points).

scaled $\left(\ln H,\ln W\right)$ has resulted in two normal random vectors of null Pearson correlation coefficient—and therefore presumably statistically independent. For comparison, the reduced scatterplots (black points) in Figure 13 are superposed on the original scatterplots (cyan points) of Figure 2 with their lines of regression (dashed blue).

The dCor values in Section III provide quantitative confirmation of the statistical independence of height and weight upon removal of the Pearson linear correlation. In the subsection A based on repeated sampling of the entire population in samples of size 50, the empirical dCor values (approximately 0.25) agree closely with dCor values (approximately 0.23 to 0.24) produced by two independent standard normal RNGs. The small deviations are attributable in part to the fact that the resulting Pearson correlation coefficients of the simulated populations were not precisely 0, but in the range 0.01 to 0.02, a consequence of the fluctuations intrinsic to finite sampling.

Potentially problematic is the discrepancy between the empirical (as well as simulated) dCor values obtained by repeated sampling and the values predicted by Equation (42), which should be close to 0 for two independent normal random variables. However, Equation (42) is strictly valid only in the limit of an infinite population. Subsection B, based on single sampling of a much larger subpopulation of 1000 participants, shows that empirical dCor values dropped to approximately 0.06, in much closer agreement with Equation (42). Evaluation of the distance correlation of a sample of 1000 required computation times longer than 8 hours. Thus, to test rigorously whether dCor approaches 0 asymptotically as a function of sample size would require impractically long computation times.

Altogether, the three sections of Table 6 consistently support the conclusion that the observed correlation between height and weight can be accounted for entirely by a bivariate lognormal distribution. In other words, the five parameters defining the bivariate lognormal distribution of height and weight suffice to predict any measureable function or test of the correlation of height and weight of a healthy adult human population.

5. Marginal Statistics of ANSUR Height and Weight Data

Previous sections concentrated on the bivariate lognormal correlation of height and weight. This section examines the marginal statistics of H and W, which are predicted to follow the respective univariate lognormal distributions $\Lambda \left(m_H,s_H^2\right)$ and $\Lambda \left(m_W,s_W^2\right)$, and on lnH and lnW, which are predicted to follow the respective univariate normal distributions $N\left(m_H,s_H^2\right)$ and $N\left(m_W,s_W^2\right)$ as discussed in Section 1.1.

Table 7 summarizes the outcomes of chi-square tests of fitness of the four variables H, W, lnH, lnW to their respective distributions, identified explicitly in column 3. As a reminder, the chi-square statistic ${\chi }_{\nu }^2$, in column 7 of the table, is determined empirically from the relation

${\left({\chi }_{\nu }^2\right)}_{emp}=\frac{1}{\kappa }{\displaystyle \underset{i=1}{\overset{\kappa }{\sum }}\frac{{\left(O_i-E_i\right)}^2}{E_i}}$ (51)

Table 7. Chi-square tests of goodness of fit of ANSUR height and weight.

where $\kappa$ is the number of test categories (bins), $O_i$ is the observed value in the i^th bin, and $E_i$ is the expected value in the i^th bin. The subscript $\nu$ (Greek nu) is the number of degrees of freedom (d.o.f.) in column 5 equal to $\kappa -1$. The chi-square tests were implemented with the Maple Statistics Package, which determined the number of bins as the integer closest to the square root of the sample size. The critical values in column 4 of the table are the values of ${\chi }_{\nu }^2$ resulting in P-values of 5%, which is the conventional threshhold of statistical significance; i.e. a tested hypothesis is deemed unsupported if the P-value is below threshhold. A P-value is the probability of obtaining a test result at least as extreme as the observed result ${\chi }_{\text{obs}}^2$, and is calculated from the expression

$P={\displaystyle \underset{{\chi }_{\text{obs}}^2}{\overset{\infty }{\int }}{p}_{{\chi }_{\nu }^2}\left(z\right)\text{d}z}$ (52)

with chi-square PDF

$p_{{\chi }_{\nu }^2}\left(z\right)=\frac{z^{\frac{1}{2}\nu -1}{\text{e}}^{-z/2}}{2^{\nu /2}\Gamma \left(\nu /2\right)}$. (53)

The outcomes summarized in the table show that all 8 propositions (the distributions of 4 variables of 2 genders) passed their respective chi-square tests with P-values far above threshhold. This means that the propositions cannot be rejected on the basis of these tests. It does not necessarily mean, however, that the propositions are true.

For further confirmation, consider again the information in Table 1. In the first section of the table, empirical values of the mean, standard deviation, skewness, and kurtosis of ANSUR heights H and weights W are compared with corresponding values predicted by lognormal expressions (7) to (10), based on the parameters in the second section of the table, derived from the variates of lnH and lnW. Agreement of experiment and theory is seen to be within 1 standard error (se) in most cases. Table 1 employed the following published estimators of the standard errors for skewness and kurtosis [36]

$se{\left(n\right)}_{Sk}=\sqrt{\frac{6n\left(n-1\right)}{\left(n-2\right)\left(n+1\right)\left(n+3\right)}}$ (54)

$se{\left(n\right)}_K=2\sqrt{\frac{6n{\left(n-1\right)}^2}{\left(n-2\right)\left(n+5\right)\left(n^2-9\right)}}$. (55)

In the second section of Table 1 the empirical skewness of both lnH and lnW for both cohorts are close to zero, as expected for normal variables. Likewise, the empirical kurtosis is very close to 3, as expected for normal variables. Skewness is a measure of the asymmetry of a distribution about the mean. Kurtosis (from the Greek root for “bulging”) is a measure of the curving or arching of the tails of a distribution; in other words, kurtosis is an indicator of the extent of outliers, relative to the normal distribution.

Ordinarily, standardized statistical moments employed in physical science and medicine include at most only the first four orders (mean, variance, skewness, kurtosis). Beyond these, higher standardized moments, such as “hyperskewness” and “hyperkurtosis” [37], are rarely used in the author’s experience, presumably because they are less readily interpetable as well as have greater measures of uncertainty for a given sample size. Nevertheless, in testing a proposed statistical distribution, it is useful to examine these higher moments, particularly if the validity of all lower moments has been confirmed.

In the terminology and notation of this paper, the hyperstatistic $S_p\left(X\right)$ of the random variable X is the p^th standardized central moment defined by the relation

$S_p\left(X\right)\equiv \frac{\langle {\left(X-{\mu }_X\right)}^p\rangle }{{\langle {\left(X-{\mu }_X\right)}^2\rangle }^{p/2}}\equiv \frac{{\mu }_p\left(X\right)}{{\mu }_2^{p/2}\left(X\right)}$ (56)

where ${\mu }_X$ is the mean of X, and ${\mu }_p\left(X\right)$ is the mean of the random variable

$m_p\left(X\right)\equiv {\left(X-{\mu }_X\right)}^p$ (57)

referred to as the p^th central moment. To simplify symbolic notation in the ensuing text, the argument (X) will be omitted whenever the context is clear.

Substitution in relation (56) of the univariate lognormal PDF, mean, and variance leads to the operational expression

$S_p\left(X\right)=\frac{{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left({\text{e}}^{su}-{\text{e}}^{\frac{1}{2}{s}^2}\right)}^p{p}_X^{\Lambda }\left(u\right)\text{d}u}}{{\left({\text{e}}^{2s^2}-{\text{e}}^{s^2}\right)}^{\frac{p}{2}}}$ (58)

where s is the standard deviation of the variable lnX. Skewness and kurtosis correspond respectively to p = 3, 4. Table 8 lists the theoretical expressions for hyperstatistics of order 1 through 6 derived from relation (58).

Table 8. Theoretical hyperstatistics of univariate lognormal distribution. (Notation: $E\left(x\right)\equiv \exp \left(x\right)$ ).

Table 9 summarizes the empirical results for hyperstatistics of orders p = 5, 6 for the variables H and W of the ANSUR population. The empirical entries (column 3) are the sample statistics

$S_p{\left(X\right)}_{emp}=\frac{\overline{m_p}}{{\left(\overline{m_2}\right)}^{\frac{p}{2}}}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(x_i-\bar{x}\right)}^p}/{\left(\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(x_i-\bar{x}\right)}^2}\right)}^{\frac{p}{2}}$ (59)

where $\bar{x}$ is the sample mean

$\bar{x}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}$ (60)

and the numerator

$\overline{m_p}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(x_i-\bar{x}\right)}^p}$ (61)

is the expectation of the sample p^th central moment (57). Theoretical entries in column 4 were calculated from Equation (58). Overall, agreement between theory and experiment appears reasonably close, but several statistics show what may be significant deviations. To ascertain whether any deviation between theory and experiment is statistically significant requires knowing the standard error of the mean statistic, but the author is unaware of any published expressions for the distributions or standard errors of hyperskewness and hyperkurtosis. To estimate the pertinent standard errors, three independent approaches were taken.

The first approach was to use the approximations of error propagation theory [38] together with expressions for variance and covariance of central moments in Chapter 10 of Ref. [28] to derive an estimate of the variance of hyperstatistic $S_p(\; X\; )$

Table 9. Test of hyperstatistics of height and weight of ANSUR population.

$\mathrm{var}\left(S_p\right)=\frac{\mathrm{var}\left(m_p\right)}{{\mu }_2^k}+\frac{p^2{\mu }_p^2}{4{\mu }_2^{k+2}}\mathrm{var}\left(m_2\right)-\frac{p{\mu }_k}{{\mu }_2^{k+1}}\mathrm{cov}\left(m_k,m_2\right)$ (62)

in which

$\mathrm{var}\left(m_p\right)={\mu }_{2p}-{\mu }_p^2+p^2{\mu }_2{\mu }_{p-1}^2-2p{\mu }_{p-1}{\mu }_{p+1}$ (63)

and

$\mathrm{cov}\left(m_p,m_q\right)={\mu }_{p+q}-{\mu }_p{\mu }_q+pq{\mu }_2{\mu }_{p-1}{\mu }_{q-1}-p{\mu }_{p-1}{\mu }_{q+1}-q{\mu }_{p+1}{\mu }_{q-1}$. (64)

The standard error is then

$se\left(S_p\right)=\sqrt{\frac{\mathrm{var}\left(S_p\right)}{n}}$. (65)

Theoretical evaluations of Equation (65) by means of the univariate lognormal PDF are recorded in column 4.

The second approach was to evaluate the means of all moments in Equation (62) by their sample expectations given by Equation (61). The resulting empirical standard errors are recorded in column 3.

In the third approach three independent simulations of the male and female ANSUR populations were made with correlated lognormal RNGs representing the variables H and W, from which empirical values of the individual hyperstatistics were obtained. The three values for each hyperstatistic (p = 5, 6) of the two variables (H, W) for the two genders (M, F) are listed in column 5 in the cells for the mean of $S_p\left(X\right)$. The sample standard error of each set of 3 independent means was calculated from the relation

$se\left(S_p\left(X\right)\right)=\frac{1}{\sqrt{n\left(n-1\right)}}\sqrt{{\displaystyle \underset{k=1}{\overset{3}{\sum }}{\left(S_p{\left(X\right)}_k-\overline{S_p\left(X\right)}\right)}^2}}$ (66)

and recorded as the first entry in the cells for standard error in column 5. Note that n is Equation (66) is 3 and not the ANSUR population size of 4082 males or 1986 females. Also, because a sample size of 3 is statistically small, one uses the unbiased sample estimate of variance in Equation (66) in which the denominator is $n-1$ [39]. The second entry (separated from the first by a vertical bar) is the difference between the largest and smallest of the three estimated means of each $S_p\left(X\right)$.

In examining the three sets of estimated standard errors, one sees that they are approximately of the same magnitude, and that the empirical values of the hyperstatistics agree with theory within $\pm 2se$ for at least one of the three estimates. However, one glaring exception is the value of $se\left(S_6{\left(W\right)}_{\text{M}}\right)$ for the male cohort obtained by simulation, which is much larger than the other standard errors for the same statistic. This occurs because of what appears to be an exceptionally high value (~38) of the mean of $S_6{\left(W\right)}_{\text{M}}$ returned by one of the simulations. This occurrence raises an important conceptual issue that calls for caution when estimating standard errors under conditions where the exact statistical distribution is unknown.

As pointed out in Chapter 10 of Ref [28], the approximations based on error propagation theory (or some variant thereof) give a valid measure of precision provided that the distribution of the statistic approaches normality in the limit of a sufficiently large sample. This is not the case for higher orders of the hyperstatistics. The mere fact that a mean value of the statistic $S_6{\left(W\right)}_{\text{M}}$ can arise within just 3 simulations that is 6 times the theoretical and 17 times the empirical standard errors of error propagation theory shows that such outlying values occur with much higher probabilities than would be predicted by a normal distribution. Under such circumstances, the appropriate way to proceed is to determine whether the empirical mean values of the hyperstatistics fall within the range between the lowest and highest corresponding statistics obtained by simulation; in other words to rely on a kind of Monte Carlo validation. Evaluated this way, all the empirical hyperstatistics in Table 9 are seen to be consistent with theory.

Taken together, the chi-square tests and agreement of theoretical and empirical moments (or functions of moments) up to the 6th power of the variables support the propositions that H and W are marginally lognormal variables.

6. Implications for Body Mass Index (BMI)

The BMI, defined by the random variable B in Equation (1), has long been used as a measure of obesity and a risk factor for associated diseases under the assumption that it is correlated with weight but largely independent of height. In Section 2 the conditional expectation function of weight, given height, was derived on the basis of lognormal theory and shown to be very nearly a quadratic power law $W\left(H\right)\propto H^2$ (27) for the ANSUR male and female cohorts. If weight varies as the square of height, then the BMI (1) would be unaffected by variations in height, or, in other words, statistically independent of height.

Figure 14 provides additional justification of the BMI assumption. The left panels of the figure display scatter plots of the correlation of scaled BMI and scaled height for the ANSUR male cohort (top), ANSUR female cohort (middle), and RNG simulated female cohort (bottom). The isotropic patterns (apart from fluctuations) are very close to what are expected for the correlation of independent vectors, as shown in the first panel of Figure 5. Quantitatively, the lines of regression (dashed red for ANSUR, dashed blue for simulation) yield Pearson correlation coefficients $2.096\times {10}^{-3}$ (male), $5.406\times {10}^{-2}$ (female), and $6.700\times {10}^{-2}$ (simulation). The three sets of variables are not standard normal variables, so a null Pearson correlation does not necessarily imply total independence.

However, evaluation of the distance correlation of the scaled variables $\left(H,B\right)$ for the male and female cohorts (black points) respectively yielded by the method of repeated sampling empirical dCor values of 0.2475 ± 0.0116 and 0.2301 ± 0.0147, which are statistically equivalent to the dCor value for correlation of independent standard normal variables generated by computer simulation

Figure 14. Scatter plots of body mass index (BMI) with height (left panels) and weight (right panels). Patterns in black were calculated from the ANSUR population of male (top panels) and female (middle panels) cohorts. Patterns in red (bottom panels) were simulated for a population of 2000 (corresponding to the size of the female cohort) by means of correlated bivariate lognormal RNGs. The slopes of the lines of regression (dashed red or dashed blue) of the left panels are close to zero, signifying independence of BMI and height. The slope of the lines of regression of the right panels are close to 0.88 for the male cohort and 0.87 for the female cohort.

(red points). It is therefore reasonable to conclude that, if the ANSUR populations can serve as baselines, then the BMI and height of healthy adult male and female populations are largely statistically independent.

By contrast, the right panels of Figure 14 (black for empirical and red for simulation) show that BMI and weight are strongly linearly correlated, which was a desirable characteristic of the BMI. The Pearson correlation coefficients are respectively 0.8814 (male cohort), 0.8706 (female cohort), and 0.8755 (simulation of female cohort) as deduced algebraically from Equation (15) or from the slopes of the lines of regression.

7. Computer Simulation of Correlated Lognormal Random Vectors

Throughout this paper computer simulation of correlated random variables (RVs) has been employed for both analytical and graphical comparisons with corresponding empirical results. A brief description of the implementation of these simulations is given in this section.

The essential objective is the simulation of a pair of correlated lognormal RVs of specified parameters $\left(m_1,m_2,s_1,s_2,r\right)$. The starting point for the construction is the well-known algebraic identity for decomposition of a general normal variable [16]

$N_1\left(m_1,s_1^2\right)=m_1+s_1{N}_1\left(0,1\right)$ (67)

$N_2\left(m_2,s_2^2\right)=m_2+s_2{N}_2\left(0,1\right)$ (68)

where $N_1\left(0,1\right)$ and $N_2\left(0,1\right)$ are independent standard normal variables (ISNVs) of mean 0 and variance 1 that serve as basis states. Each ISNV represents a random number generator (RNG) of one’s mathematical software. Populations of size n are simulated by creating sets of n variates from each ISNV. The covariance $\langle \left(N_1\left(m_1,s_1^2\right)-m_1\right)\left(N_2\left(m_2,s_2^2\right)-m_2\right)\rangle$ is theoretically 0 (because $\langle N_1\left(0,1\right)N_2\left(0,1\right)\rangle$ is zero) and empirically should approach 0 numerically in the limit of increasing sample size n.

To create a normal RV $N_{2c}\left(m_2,s_2^2,r\right)$ correlated with the RV in Equation (67), one makes the following linear superposition

$N_{2c}\left(m_2,s_2^2,r\right)=m_2+r{s}_2{N}_1\left(0,1\right)+\sqrt{1-r^2}{s}_2{N}_2\left(0,1\right)$. (69)

Note that, according to the algebraic rules [16] that govern manipulation of independent normal RVs,

$a{N}_1\left(0,1\right)+b{N}_2\left(0,1\right)=N_1\left(0,a^2\right)+N_2\left(0,b^2\right)=N\left(0,a^2+b^2\right)$, (70)

where a and b are constants, one could combine the second and third terms in the right side of Equation (69) to recover the marginal distribution represented by Equation (68), since the correlation parameter r drops out. The set of RVs (67) and (69) then comprise a pair of correlated bivariate normal RVs, which, when implemented numerically on a computer, generate the respective variates $\left(y_{1,}{}_i,y_{2,}{}_i\right)$ for $i=1,\cdots ,n$. In the context of simulating samples of human height and weight, these variates are

$\begin{array}{l}{y}_{1,i}=\ln \left(h_i\right)\\ y_{2,i}=\ln \left(w_i\right)\end{array}$ (71)

Once one has created the sets of normal variates (71), it remains only to exponentiate them

$\begin{array}{l}{x}_{1,i}=\exp \left(y_{1,i}\right)=h_i\\ x_{2,i}=\exp \left(y_{2,i}\right)=w_i\end{array}$ (72)

to simulate a sample of correlated bivariate lognormal RVs. With some mathematical applications such as Maple, one can work with the abstract vectors and simply enter a statement like $X=\exp \left(Y\right)$, whereupon the application will know to exponentiate the individual variates as in relation (72). Thus, the bivariate lognormal vectors corresponding to Equations (67) and (69) generated by Maple were defined and implemented by the forms

${\Lambda }_1\left(m_1,s_1^2\right)=\exp \left(m_1+s_1{N}_1\left(0,1\right)\right)$ (73)

${\Lambda }_{2c}\left(m_2,s_2^2,r\right)=\exp \left(m_2+r{s}_2{N}_1\left(0,1\right)+\sqrt{1-r^2}{s}_2{N}_2\left(0,1\right)\right)$ (74)

For completeness, a final comment as to the actual nature of the RNGs employed in this paper is called for. All RNGs that employ a mathematical algorithm, in contrast to RNGs based on some random quantum process such as radioactive decay [40], generate pseudo-random numbers. These are sets of numbers that are generated reproducibly from a known starting point (seed value), yet nevertheless pass diverse statistical tests for randomness. The more stringent a test, and the more tests a RNG passes, the better is the RNG.

The MersenneTwister algorithm supplied by the Maple RandomTools Subpackage has passed the diehard tests of randomness [41] by G. Marsaglia as well as other tests, and provides numbers that can be considered cryptographically secure [42] [43]. It is safe to accept, therefore, that the independent normal basis states with which the simulation algorithm began and from which the correlated bivariate lognormal distributions were created were for all practical purposes sufficiently uncorrelated.

8. Conclusions and Interpretation

Knowledge of the exact distribution function of a random quantity provides the most complete statistical information attainable about that quantity. This is especially important in regard to anthropometric attributes the statistics of which are essential to clinical medicine and epidemiology.

The fundamental conclusion of this paper is that human height H and weight W in a population of healthy adults are statistically distributed as correlated bivariate lognormal random variables. Moreover, for all practical purposes, this distribution is thought not to be approximate, but empirically rigorous in samples of sufficient size. This means that five measurable parameters, comprising two means $\left(m_H,m_W\right)$, two variances $\left(s_H^2,s_W^2\right)$, and the Pearson linear correlation coefficient (r), suffice to determine all statistical attributes (probabilities, moments, correlations) regarding the relation of height and weight in a specified population.

In support of this conclusion, detailed statistical analyses of an extensive anthropometric data base of diverse individuals have shown the following:

• The variates of H and W of both gender cohorts satisfy chi-square tests of fitness to univariate lognormal distributions.

• The variates of lnH and lnW of both gender cohorts satisfy chi-square tests of fitness to univariate normal distributions.

• The sample moments of $\left(H,W\right)$ up to 6th order are consistent with the lognormal distribution.

• The sample moments of $\left(\ln H,\ln W\right)$ up to 6th order are consistent with the normal distribution.

• Theoretical correlation functions $R_{p,q}\left(r\right)$ of the normal distribution predict correctly the sample correlation coefficients $r_{p,q}$ of the variates of $\left(\ln H,\ln W\right)$.

• Theoretical correlation functions $C_{p,q}\left(r,s_H,s_W\right)$ of the lognormal distribution predict correctly the sample correlation coefficients $c_{p,q}$ of the variates of $\left(H,W\right)$.

• Computer simulations using correlated lognormal random number generators (RNGs) produce correlation functions and probability density plots of the variables $\left(H,W\right)$ that statistically match the corresponding empirical functions and plots.

• Computer simulations using correlated normal (RNGs) produce correlation functions and probability density plots of the variables $\left(\ln H,\ln W\right)$ that statistically match the corresponding empirical functions and plots.

• Empirically measured distance correlation (dCor) values of the variables $\left(H,W\right)$ agree with dCor values obtained from comparably sized populations simulated by correlated lognormal RNGs.

• Empirically measured distance correlation (dCor) values of the variables $\left(\ln H,\ln W\right)$ agree with dCor values obtained from comparably sized populations simulated by correlated normal RNGs.

• Removal of the line of regression in the empirical density plot of lnW vs lnH produces an isotropic density with null Pearson correlation coefficient. The density plot, null Pearson correlation coefficient, and dCor values statistically match the corresponding outcomes from two independent normal RNGs.

In short, taken altogether, the preceding extensive set of tests supports the proposition that the distribution and correlation of height and weight of a healthy adult human population are fully accounted for by a bivariate lognormal distribution.

A secondary point worth noting, given the importance of body mass index (BMI) to current medicine and epidemiology, is that the conditional expectation of weight, given height, theoretically derived from the lognormal distribution function yielded functional relations (23), (24) between weight and height. These functions, when evaluated with the lognormal parameters of the ANSUR male and female cohorts, led in both cases to a nearly exact quadratic power law (27), thereby justifying theoretically a long-held assumption underlying the use of BMI as a risk factor for obesity-related diseases.

In concluding this paper, it is useful to clarify what is meant by an “exact” statistical distribution. In the opinion of the author, who is an atomic and nuclear physicist, statistical distributions in science can arise, broadly speaking, in two ways.

The most fundamental way is as a consequence of a particular dynamical model. In physics, for example, the decay of radioactive nuclei is rigorously accounted for by a binomial probability function, based on a physical model of the independent decay of discrete, uncorrelated nuclei [44]. If the assumption of independence were found to be invalid—and there have been a considerable number of such challenges, only to have been debunked by more careful experiment and analysis [45] [46] [47] [48] —the discovery would have led to deep new insights into the structure and behavior of matter.

The second, less fundamental way, but nevertheless one of practical utility, is by empirical recognition and subsequent verification. To return to the previous physics example, suppose that the phenomenon of radioactive decay was discovered before there was any understanding or general acceptance of atoms as discrete units of matter³. Then radioactive decay would have been empirically observed to be a Poisson process, and, indeed, the Poisson distribution is widely depicted in books as a rigorous physical law. (See, for example [49].) However, in retrospect, a Poisson process can be interpreted as a degenerate case of a binomial process in the limit of a large number N of radioactive atoms with low probability p of decay, such that Np is the mean number of decays within a specified time interval.

The point of the foregoing examples is this: The rigorously exact distribution (binomial) revealed critical information about the constituents (discrete, independent) of the system. The apparently exact distribution (Poisson) was empirical and utilitarian, but revealed little about the system other than that the decay products were discrete. Under appropriately conceived radiation experiments, the difference between the binomial and Poisson distributions can be observed [50], and the fundamentality of the binomial distribution is established.

In regard to the statistical attributes of human height and weight, the consistency with a correlated bivariate lognormal distribution is, as shown in this paper, so extensive and close, that one must wonder whether it is a rigorously exact consequence of some biophysical mechanism or a limiting case of some other statistical process. How, for example, might a lognormal distribution arise from other distributions?

One such process might entail a random variable X comprising a product of some set of arbitrarily distributed random variables, in which case application of the Central Limit Theorem to lnX could result in a normal distribution. Then the parent variable X would itself be lognormal. It is difficult to conceive in detail, however, of mechanisms by which real biological processes responsible for human height and weight could engender such a hypothetical X as to produce a correlated bivariate lognormal distribution.

More generally, a lognormal distribution can also arise under circumstances where an intrinsically positive variable has a low mean and high variance, leading to a pronounced skewness. However, any of a large number of other skewed distributions could also arise, so the mechanism is not unique. Moreover, as demonstrated in this paper, whatever mechanism is invoked must produce not only the correct skewness, but also kurtosis and other hyperstatistics as well.

At this stage and until testable mechanisms are proposed, refutation of the exactness of the correlated bivariate lognormal distribution of human height and weight can only come from further detailed statistical analysis of larger populations. And if such future tests further confirm the exactness of the bivariate lognormal relation of height and weight, then, like the example of radioactivity cited above, this knowledge will have revealed something fundamental about the physical processes underlying human development.

Acknowledgements

The author thanks Trinity College for partial support through the research fund associated with the George A. Jarvis Chair of Physics.

Appendix—Glossary of Abbreviations

BMI—Body Mass Index

CDF—Cumulative Distribution Function

CF—Characteristic Function

d.o.f.—Degrees of Freedom

dCor—Distance Correlation

dCov—Distance Covariance

dVar—Distance Variance

ISNV—Independent Standard Normal Variable

PDF—Probability Density Function

RNG—Random Number Generator

RV—Random Variable

NOTES

¹These transformed functions are the characteristic function (CF), which is the Fourier transform of the PDF, and the cumulative distribution function (CDF), which is the integral of the PDF from some fixed point to the argument of the PDF. Thus, if $g\left(x\right)$ is the CDF, then the PDF $f\left(x\right)=\text{d}g\left(x\right)/\text{d}x$.

²Integration variables, in contrast to random variables, will be represented by lower case letters.

³This supposition is actually historically correct. Radioactivity was discovered by Henri Bequerel in 1896, whereas opposition to the existence of atoms by some leading scientists of the day lasted until about 1910.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Stigler, S.M. (1986) The History of Statistics: The Measurement of Uncertainty Before 1900. Harvard University Press, Cambridge, 265-361.
[2]	Bernstein, P.L. (1998) Against the Gods: The Remarkable Story of Risk. Wiley, New York, 152-171. https://doi.org/10.2307/2685740
[3]	Porter, T.M. (2004) Karl Pearson: The Scientific Life in a Statistical Age. Princeton University Press, Princeton, 235-239, 249-266. https://doi.org/10.5944/empiria.8.2004.989
[4]	Quetelet, L.A.J. (1835) A Treatisse on Man and the Development of His Faculties. Cambridge University Press, Cambridge. https://www.cambridge.org/core/books/treatise-on-man-and-the-development-of-his-faculties/AB13A647A6C8727C06AE5399D7422887
[5]	Sager, G. (1987) Relation between Body Height and Weight in Adult Humans. Gegenbaurs morphologisches Jahrbuch, 133, 563-571.
[6]	Rahmandad, H. (2014) Human Growth and Body Weight Dynamics: An Integrative Systems Model. PLOS ONE, 9, e114609. https://doi.org/10.1371/journal.pone.0114609 https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0114609
[7]	Lettre, G. (2011) Recent Progress in the Study of the Genetics of Height. Human Genetics, 129, 465-472. https://doi.org/10.1007/s00439-011-0969-x
[8]	Wikipedia (2022) Body Mass Index. https://en.wikipedia.org/wiki/Body_mass_index
[9]	World Health Organization (2021) Obesity and Overweight. https://www.who.int/news-room/fact-sheets/detail/obesity-and-overweight
[10]	Moody, J.N., et al. (2021) Body Mass Index and Polygenic Risk for Alzheimer’s Disease Predict Conversion to Alzheimer’s Disease. The Journals of Gerontology Series A Biological Sciences and Medical Sciences, 76, 1415-1422. https://doi.org/10.1093/gerona/glab117
[11]	Silverman, M.P. and Lipscombe, T.C. (2022) Exact Statistical Distribution of the Body Mass Index (BMI): Analysis and Experimental Confirmation. Open Journal of Statistics, 12, 324-356. https://doi.org/10.4236/ojs.2022.123022
[12]	Szekely, G.J., Rizzo, M.L. and Bakirov, N.K. (2007) Measuring and Testing Dependence by Correlation of Distances. The Annals of Statistics, 35, 2769-2794. https://doi.org/10.1214/009053607000000505
[13]	Szekely, G.J. and Rizzo, M.L. (2009) Brownian Distance Covariance. The Annals of Applied Statistics, 3, 1236-1265. https://doi.org/10.1214/09-AOAS312
[14]	Callahan, A. (2021) Is BMI a Scam? The New York Times. https://www.nytimes.com/2021/05/18/style/is-bmi-a-scam.html
[15]	Gordon, C.C., et al. (2014) 2012 Anthropometric Survey of U.S. Army Personnel: Methods and Summary Statistics. Technical Report Natick/TR-15/007, U.S. Army Natick Soldier Research and Engineering Center, Natick. https://www.openlab.psu.edu/ansur2
[16]	Silverman, M.P. (2014) A Certain Uncertainty: Nature’s Random Ways. Cambridge University Press, Cambridge, 511-514. https://doi.org/10.1017/CBO9781139507370
[17]	Forbes, C., Evans, M., Hastings, N. and Peacock, B. (2011) Statistical Distributions. 4th Edition, Wiley, New York, 131-134. https://doi.org/10.1002/9780470627242
[18]	A’Hearn, B., Peracchi, F. and Vecchi, G. (2009) Height and the Normal Distribution: Evidence from Italian Military Data. Demography, 46, 1-25. https://doi.org/10.1353/dem.0.0049
[19]	Diverse Populations Collaborative Group (2005) Weight-Height Relationships and Body Mass Index: Some Observations from the Diverse Populations Collaboration. The American Journal of Physical Anthropology, 128, 220-229. https://doi.org/10.1002/ajpa.20107
[20]	Johnson, W., et al. (2020) Differences in the Relationship of Weight to Height, and Thus the Meaning of BMI According to Age, Sex, and Birth Year Cohort. Annals of Human Biology, 47, 199-207. https://doi.org/10.1080/03014460.2020.1737731
[21]	Sperrin, M., Marshall, A.D., Higgins, V., Renehan, A.G. and Buchan, I.E. (2015) Body Mass Index Relates Weight to Height Differently in Women and Older Adults: Serial Cross-Sectional Surveys in England (1992-2011). Journal of Public Health, 38, 607-613. https://doi.org/10.1093/pubmed/fdv067
[22]	Benn, R.T. (1971) Some Mathematical Properties of Weight-for-Height Indices Used as Measures of Adiposity. Journal of Epidemiology & Community Health, 25, 42-50. https://doi.org/10.1136/jech.25.1.42
[23]	Rohrer, F. (1921) Der Index der Körperfülle als Maß des Ernährungszustandes [The Index of Corpulence as a Measure of Nutritional Condition]. Münchener Medizinische Wochenschrift, 68, 580-582.
[24]	Henneberg, M., Hugg, J. and Townsend, E.J. (1989) Body Weight/Height Relationship: Exponential Solution. American Journal of Human Biology, 1, 483-491. https://doi.org/10.1002/ajhb.1310010412
[25]	Cidras, M. (2015) Body Mass Exponential Index: An Age-Independent Anthropometric Nutritional Assessment. Open Access Library Journal, 2, 1-8. https://doi.org/10.4236/oalib.1101943
[26]	Trussell, J. and Bloom, D.E. (1979) A Model Distribution of Height or Weight at a Given Age. Human Biology, 51, 523-536.
[27]	Edwards, A.W.F. (1992) Likelihood. The Johns Hopkins University Press, Baltimore, 70-143.
[28]	Kendall, M.G. and Stuart, A. (1963) The Advanced Theory of Statistics Vol. 1: Distribution Theory. Hafner, New York, 94-119, 228-236.
[29]	Hotelling, H. (1953) New Light on the Correlation Coefficient and Its Transforms. Journal of the Royal Statistical Society: Series B, 15, 193-232. https://doi.org/10.1111/j.2517-6161.1953.tb00135.x
[30]	Mood, A.M., Graybill, F.A. and Boes, D.C. (1974) Introduction to the Theory of Statistics. 3rd Edition, McGraw-Hill, New York, 195-198, 233-236.
[31]	Chou, Y. (1969) Statistical Analysis: With Business and Economic Applications. Holt, Rinehart, and Winston, New York, 308-323.
[32]	Haldane, J.B.S. (1942) Moments of the Distributions of Powers and Products of Normal Variates. Biometrika, 32, 226-242. https://doi.org/10.1093/biomet/32.3-4.226
[33]	Arfken, G.B. and Weber, H.J. (2005) Mathematical Methods for Physicists. 6th Edition, Elsevier, New York, 83-87.
[34]	Wikipedia (2022) Bootstrapping (Statistics). https://en.wikipedia.org/wiki/Bootstrapping_(statistics)
[35]	Efron, B. (1979) Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics, 7, 1-26. https://doi.org/10.1214/aos/1176344552
[36]	Harding, B., Tremblay, C. and Cousineau, D. (2014) Standard Errors: A Review and Evaluation of Standard Error Estimators Using Monte Carlo Simulations. The Quantitative Methods for Psychology, 10, 107-123. https://doi.org/10.20982/tqmp.10.2.p107
[37]	Moment (Mathematics) (2022) Wikipedia. https://en.wikipedia.org/wiki/Moment_(mathematics)
[38]	Silverman, M.P., Strange, W. and Lipscombe, T.C. (2004) The Distribution of Composite Measurements: How to Be Certain of the Uncertainties in What We Measure. American Journal of Physics, 72, 1068-1081. https://doi.org/10.1119/1.1738426
[39]	Kendall, M.G. and Stuart, A. (1961) The Advanced Theory of Statistics Vol. 2: Inference and Relationship. Charles Griffin & Co., London, 1-8.
[40]	Walker, J. (1996) HotBits: Genuine Random Numbers, Generated by Radioactive Decay. https://www.fourmilab.ch/hotbits
[41]	Wikipedia (2022) Diehard Tests. https://en.wikipedia.org/wiki/Diehard_tests
[42]	Maplesoft.com (2022) Overview of the RandomTools [MersenneTwister] Subpackage. https://www.maplesoft.com/support/help/maple/view.aspx
[43]	Mapleprimes.com (2019) Are Maple’s Pseudo Random Number Generators Good Generators? Post by mmcdara 3900. https://mapleprimes.com/posts/211598-Are-Maples-Pseudo-Random-Number-Generators
[44]	Lapp, R.E. and Andrews, H.L. (1972) Nuclear Radiation Physics. Prentice-Hall, Englewood Cliffs, 36-40.
[45]	Silverman, M.P., Strange, W., Silverman, C.R. and Lipscombe, T.C. (1999) Tests of Alpha-, Beta-, and Electron Capture Decays for Randomness. Physics A, 262, 265-273. https://doi.org/10.1016/S0375-9601(99)00668-4
[46]	Silverman, M.P. and Strange, W. (2009) Search for Correlated Fluctuations in the Decay of Na-22. Europhysics Letters, 87, Article No. 32001. https://doi.org/10.1209/0295-5075/87/32001
[47]	Silverman, M.P. (2015) Search for Non-Standard Radioactive Decay Based on Distribution of Activities. Europhysics Letters, 110, Article No. 52001. https://doi.org/10.1209/0295-5075/110/52001
[48]	Silverman, M.P. (2016) Search for Anomalies in the Decay of Radioactive Mn-54. Europhysics Letters, 114, Article No. 62001. https://doi.org/10.1209/0295-5075/114/62001
[49]	Miller, D.G. (1972) Radioactivity and Radiation Detection. Gordon and Breach, New York, 88-99.
[50]	Foster, J., Kouris, K., Matthews, I.P. and Spyrou, N.M. (1983) Binomial vs Poisson Statistics in Radiation Studies. Nuclear Instruments and Method, 212, 301-305. https://doi.org/10.1016/0167-5087(83)90706-8