1.2.1 Distinguishing Between Cause and Effect

Does money make people happy, or do happier people make more money (Diener & Biswas-Diener, Reference Diener and Biswas-Diener2002)? Do violent video games cause consumers to be aggressive, or are people with aggressive tendencies more attracted to violent video games (Gentile et al., Reference Gentile, Lynch, Linder and Walsh2004)? Does social media use reduce the well-being of users, or are unhappy people more likely to engage in social media (Kross et al., Reference Kross, Verduyn, Demiralp, Park, Lee, Lin, Shablack, Jonides and Ybarra2013)? Does counterproductive work behavior lead to sleep disturbances or do people with poor sleep quality tend to be counterproductive in their work behavior (Shi et al., Reference Shi, Fairchild and Wiedermann2023)? Does smoking contribute to the development of depression, or are people with depressive symptoms more likely to engage in health-damaging behavior (Munafò & Araya, Reference Munafò and Araya2010)?

These are typical questions raised by researchers in the social sciences when relying on observational data to test their hypotheses (e.g., due to ethical or financial constraints to perform experiments). For each question, it is, in principle, conceivable that at least two causal explanations exist for a variable association. In other words, in each case the hypothesized causal mechanism is, in theory, reversible. That is, for two variables x and y, a causal model of the form x → y (x causes y) or the causally competing model y → x (y causes x) can be entertained to explain the link between x and y. Unfortunately, standard covariance-based methods to evaluate variable associations (such as correlation analysis and linear regression models) do not help one to distinguish causally competing models (Chapter 2 introduces correlation and regression in their symmetric forms).

In contrast, methods of DDA have specifically been developed for the task of analyzing reversible data scenarios. In contrast to standard methods of association, DDA makes use of higher than second moments (e.g., skewness, kurtosis, co-skewness, and co-kurtosis) to test causally competing models. Provided that certain data conditions are fulfilled, the framework of DDA enables one to probe which causal explanation (x → y or y → x) finds more support from the data. Three DDA components are available: (1) distributional properties of observed variables; (2) distributional properties of errors of competing models; and (3) independence properties of predictors and errors. Chapters 3, 4, and 5 introduce the three components in detail. Specifically, each chapter (a) presents test statistics compatible with principles of direction of dependence, (b) summarizes guidelines and decision rules for causal model selection, and (c) illustrates the application of DDA components using synthetic and real-world data examples.

1.2.2 Identifying the Presence of Hidden Confounders

When relying on observational (nonexperimental) data, confounders (i.e., common causes of x and y; sometimes called lurking or omitted variables) are likely to be present which affect the causal relation between focal variables. In the worst case, an observed association between two variables may entirely be attributable to the presence of a confounder. For example, common genetic factors may (in parts) be responsible for depressive symptomology and health-damaging behavior (Tsuang et al., Reference Tsuang, Francis, Minor, Thomas and Stone2012), or parenting style may affect both violent video game exposure and aggression (Cote et al., Reference Cote, Coles and Dal Cin2021). Furthermore, such common causes may not be available to the researchers.

In observational studies, hidden confounders are omnipresent. Here, the main issue is that causal effect estimates are biased whenever common factors are present that are not statistically accounted for. Depending on the sign of the confounding effects, the bias introduced by confounders can go in either direction, that is, the causal effect estimate of interest can be smaller or larger than the true population effect. Thus, including confounders is crucial to guarantee consistent causal effect estimation. In practice, since researchers almost never have access to all influential common factors, the question is whether the available set of covariates is appropriate to sufficiently de-confound the relation between a focal predictor and a focal outcome. The methods presented in this book can, under certain data conditions, be used to answer this question. Here, the key ingredient for confounder detection is that unconsidered common factors create systematic dependence between predictors and model errors. Provided that the variables under study deviate from the Gaussian distribution, non-independence between predictors and errors can be detected. Chapter 6 discusses the identification of confounders in detail.

The task of proper confounder control is further complicated by the fact that availability alone is not an appropriate criterion for covariate selection. Covariates can, for example, take the role of a collider variable (i.e., a variable that is caused by the predictor and outcome under study). Here, it is well-known that adjusting for a collider variable induces bias into the causal effect estimate of interest (Elwert & Winship, Reference Elwert and Winship2014). Distinguishing eligible (confounders) and non-eligible covariates (colliders) is, thus, a critical element in causal modeling. Chapter 7 discusses covariate selection in the context of direction of dependence, shows how one can empirically identify eligible covariates, and integrates the components of DDA into a cohesive framework.

1.2.3 Identifying Data Sectors of Predominant Causal Models

Consider the development of brain and behavior. This development can be described by processes in a dynamic network of components (see, e.g., Fischer & van Geert, Reference Fischer, van Geert, Molenaar, Lerner and Newell2014) which are subject to change over time. These changes can involve the strength as well as the causal orientation of the connection between components. For example, let x be a child’s skill of adding numbers and y the child’s abstract understanding of addition. Since a certain level of skill of addition is required to develop an abstract understanding of addition, in early developmental stages, x can be conceptualized as a causal precursor of y (or, in path notation, x → y). Gaining an abstract understanding of what addition means, however, will lead to improved addition skills. In other words, in later stages, a causal connection of the form y → x is established. The stage of development, thus, takes the role of a causal path modifying factor. Now, suppose one uses the child’s age as a proxy for the stage of development. Thus, conditioned on the age of the child, the causal path may either go from x to y or from y to x. The methods presented in Chapter 9 (i.e., extensions of DDA to conditional causal models) can help evaluating such structural changes through identifying data sectors for which x → y (or y → x) constitutes the predominant causal model.

Causal heterogeneity can come in many different facets. In addition to potential heterogeneity of the magnitude and the direction of causation, the influence of hidden common causes can also vary from subpopulation to subpopulation. For example, in studying the causal connection between a child’s skill in adding numbers and the child’s abstract understanding of addition, factors that affect both (e.g., skill acquisition in other areas of numerical cognition) can become more prominent in later stages. In addition, contextual factors (e.g., additional educational experiences that facilitate mathematics comprehension) can change over time. In other words, depending on the developmental stage of the child, a confounder may be more or less active. Such subtle differences in the causal structure are, however, important when one aims at identifying the most plausible explanation for an underlying data-generating mechanism – in particular, when information on confounding factors is not available to the researcher. Methods of conditional DDA, introduced in Chapter 9, can help to answer for whom and under what conditions one can expect hidden confounding to be most pronounced.

1.2.4 Identifying Mechanisms of Mediation

Statistical mediation analysis (see, e.g., Baron & Kenny, Reference Baron and Kenny1986; MacKinnon, Reference MacKinnon2008; Wiedermann & von Eye, Reference Wiedermann and von Eye2015a; Wright, Reference Wright1934) constitutes a valuable tool in testing causal pathways of variable relations. Here, in addition to a tentative predictor and a tentative outcome, a third variable (the mediator) is considered which is hypothesized to transmit the effect from the predictor to the outcome. In other words, one decomposes the effect of the predictor on the outcome into a direct and indirect effect component. To interpret direct and indirect effects as causal, a series of assumptions must be met: First, the mediational structure must be correctly specified. For three variables (the predictor, the mediator, and the outcome), however, various acyclic models can be specified and, from a purely statistical perspective, each model fits the data equally well (known as Markov equivalence). Second, confounders must be properly accounted for. In other words, even when the mediation model is correctly specified with respect to the direction of causation, unconfoundedness assumptions are required to endow direct and indirect effects with causal meaning. The methods presented in this book are designed to critically evaluate both absence of reverse causation and confounding biases. Chapter 9 outlines statistical principles and decision guidelines of direction of dependence in mediation models.

1.2.5 How Robust Are the Causal Claims?

The output of any mathematical model of empirical data depends to some degree on the characteristics of the input. Consequently, a critical examination of the stability and robustness of model results constitutes an indispensable element of any data analysis. In the context of direction of dependence, stability analysis can be used to address questions concerning the robustness of the discerned causal claims against sample composition. Furthermore, sensitivity analysis addresses questions concerning the robustness of DDA results against violated assumptions. Chapter 8 presents approaches to test the stability and sensitivity in the context of probing the causal direction of variable relations.

1.3.1 Types of Moments of Variables

In general, moments describe the distribution of probability mass. This distribution is determined by parameters that describe the location, scale, and shape of a distribution. Specifically, suppose that x is a random variable with density function $p_x\left(x\right)$. The k-th moment of x around a constant (nonrandom) value c can be written as

$E\left[{\left(x-c\right)}^k\right]={\int }_{\mathit{-}\mathit{\infty }}^{\mathit{\infty }}{\left(x-c\right)}^k{p}_x\left(x\right)\mathit{dx},$(1.5)

with E denoting the so-called expected value operator. Depending on the constant c, one can distinguish between raw moments and central moments. Raw moments describe the moments of a variable around the origin of zero (i.e., c = 0), whereas central moments describe moments around the arithmetic mean of the distribution, i.e., $c=E\left[x\right]={\mathrm{\mu }}_x$. Thus, the k-th central moment of a variable can be expressed as

$m_k=E\left[{\left(x-{\mathrm{\mu }}_x\right)}^k\right]={\int }_{\mathit{-}\mathit{\infty }}^{\mathit{\infty }}{\left(x-{\mathrm{\mu }}_x\right)}^k{p}_x\left(x\right)\mathit{dx}.$(1.6)

Central moments have the advantage that, due to subtracting the mean ${\mathrm{\mu }}_x$, statements about the probability distribution of x are mean (location)-invariant. Suppose, for example, one wants to compare the variability of daily temperature for two seasons, winter and summer. Because, in summer, the earth is tilted toward the sun and, in winter, it is tilted away from the sun, one expects (on average) higher temperatures in the summer and colder temperatures in the winter. To compare the fluctuation in temperature for the two seasons one must, therefore, account for the average temperature. Here, central moments serve this purpose.

In addition to central moments around the mean, one can define so-called standardized moments. The k-th standardized moment is given by

${\overline{m}}_k=\frac{m_k}{{\mathrm{\sigma }}^k}=E\left[{\left(\frac{x-{\mathrm{\mu }}_x}{{\mathrm{\sigma }}_x}\right)}^k\right],$(1.7)

with $m_k=E\left[{\left(x-{\mathrm{\mu }}_x\right)}^k\right]$ and ${\mathrm{\sigma }}^k$ denoting the standard deviation of x (i.e., ${\mathrm{\sigma }}_x=\sqrt{E\left[{\left(x-{\mathrm{\mu }}_x\right)}^2\right]}$; for further details see Section 1.3.3) taken to the k-th power. Dividing ($x-{\mathrm{\mu }}_x$) by the standard deviation, i.e., $\left(x-{\mathrm{\mu }}_x\right)/{\mathrm{\sigma }}_x$, leads to standardized scores of x. As a result, k-th standardized moments are both, location (mean)- and scale (variance)-invariant, which allows statements about the distribution of the probability mass along the tails of the distribution, regardless of the location and the spread of the x scores. In other words, third-, fourth-, and higher than fourth standardized moments quantify the tailedness of the distribution while accounting for the center of the distribution (the mean) and the variation around the center (variance). The following sections discuss the first four moments (mean, variance, skewness, and kurtosis) in detail.

1.3.2 Mean

The first raw moment, also known as the expectation of x, describes the center of mass of a probability distribution and is given by

$E\left[x\right]={\mathrm{\mu }}_x={\int }_{\mathit{-}\mathit{\infty }}^{\mathit{\infty }}x{p}_x\left(x\right)\mathit{dx}.$(1.8)

In other words, the first raw moment describes the one point where the probability distribution is perfectly balanced. Alternatively, we can interpret the first moment as the distance of the center of the distribution from the origin of zero. The corresponding sample estimate is given by the well-known arithmetic mean equation

${\widehat{\mathrm{\mu }}}_x=\frac{1}{n}{\sum }_i{x}_i$(1.9)

with i = 1, …, n where n denotes the sample size.

To illustrate, consider Galton’s (Reference Galton1886) classical height data, which were used to perform one of the first regression analyses in the history of statistics. The dataset contains height measures (in inches) from 928 (adult) children and the average heights of the two parents (“midparents”). The present illustration only focuses on the height of the (adult) children. The left panel of Figure 1.1 gives the density distribution of children’s heights. The vertical line gives the first raw moment, that is, the mean of children’s heights of 68.09 inches (172.95 cm). In other words, at a height of 68.09 inches, the probability distribution is perfectly balanced. Next, consider the right panel of Figure 1.1. Here heights of the 928 children were subtracted by the mean prior visualization. That is, instead of x, we plot the quantities x – 68.09. Because subtracting a constant (here, the mean) constitutes a linear transformation, the overall shape of the distribution is not affected by this subtraction. However, this transformation systematically changes the center mass of the distribution. The mean of the distribution is shifted to the origin of zero. In other words, the first centered moment of ($x-{\mathrm{\mu }}_x$) is of lower interest because it always takes the value zero, that is,

$\begin{array}{l}E\left[\left(x-{\mathrm{\mu }}_x\right)\right]=E\left[x\right]-E\left[{\mathrm{\mu }}_x\right]\\ ={\mathrm{\mu }}_x-{\mathrm{\mu }}_x\\ =0,\end{array}$(1.10)

where we make use of the fact that the expected value of a constant equals the constant itself ($E\left[{\mathrm{\mu }}_x\right]={\mathrm{\mu }}_x$). Section 1.3.3 focuses on the variance of a probability distribution.

Figure 1.1 Density of Galton’s data on the height of 928 children. Left panel: Density distribution of the original data. Right panel: Density distribution of mean-standardized scores.

Figure 1.1Long description

The left plot is labelled Children’s Height on the x-axis and Density on the y-axis, with values on the x-axis ranging from 60 to 75. It displays a grey-shaded kernel density curve with a vertical line near 68, indicating the sample mean. The boundary of the shaded area passes through: (60, 0), (65, 0.06), (70, 0.17), and (75, 0). The right plot is labelled Mean-Standardised Children’s Height on the x-axis and Density on the y-axis, with x-values ranging from minus 5 to plus 5. It presents a similarly shaped grey-shaded density curve with a vertical line at 0, representing the mean of the standardised data. The boundary of the shaded area passes through (negative 8, 0.025), (negative 4, 0.05), (1, 0.017), and (7, 0.025). Both plots display a right-skewed distribution. All values are estimated.

1.3.3 Variance

The second central moment describes the variance of a probability distribution and measures the spread of the distribution about the mean. Formally, the variance can be expressed as

$m_2={\mathrm{\sigma }}_x^2=E\left[{\left(x-{\mathrm{\mu }}_x\right)}^2\right]={\int }_{\mathit{-}\mathit{\infty }}^{\mathit{\infty }}{\left(x-{\mathrm{\mu }}_x\right)}^2{p}_x\left(x\right)\mathit{dx}.$(1.11)

In words, the variance is defined as the expected value of the squared deviation of x about the mean. Note that the variance is a measure of the spread of the distribution in squared units of the data. Taking the square root of the variance leads to the standard deviation of x,

${\mathrm{\sigma }}_x=\sqrt{E\left[{\left(x-{\mathrm{\mu }}_x\right)}^2\right]},$(1.12)

which serves as a measure of the spread of the distribution in the measurement unit of x. For both measures it holds that larger values indicate larger variability of the scores about the center (mean) of the probability distribution. In contrast, a variance (standard deviation) of zero implies that all x scores take the same value. The corresponding sample estimates of the variance and the standard deviation are given by

$\begin{array}{l}{\widehat{\mathrm{\sigma }}}_x^2=\frac{1}{n}{\sum }_i{\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)}^2\\ {\widehat{\mathrm{\sigma }}}_x=\sqrt{\frac{1}{n}{\sum }_i{\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)}^2}.\end{array}$(1.13)

Reconsider Galton’s height data of the 928 (adult) children. The left panel of Figure 1.2, again, gives the density distribution of the raw data. However, now, we focus on both the center and the spread of the distribution. Here, a variance of ${\widehat{\mathrm{\sigma }}}_x^2$ = 6.34 and, thus, a standard deviation of ${\widehat{\mathrm{\sigma }}}_x=\sqrt{6.34}$ = 2.52 are observed. In other words, on average, children’s heights deviate by 2.52 inches from the mean (visualized through horizonal lines in Figure 1.2). Next, we focus on the height distribution given in the right panel of Figure 1.2. Here, prior data visualization, data points have been standardized using $\left(x-{\widehat{\mathrm{\mu }}}_x\right)/{\widehat{\mathrm{\sigma }}}_x$ with ${\widehat{\mathrm{\mu }}}_x$ = 68.09 and ${\widehat{\mathrm{\sigma }}}_x$ = 2.52. Thus, the resulting probability distribution again exhibits a mean of zero. In addition, the probability distribution shows a variance (and standard deviation) of one, which is a necessity due to the fact that

$\begin{array}{l}{\overline{m}}_2=\frac{m_2}{{\mathrm{\sigma }}^2}\\ =\frac{E\left[{\left(x-{\mathrm{\mu }}_x\right)}^2\right]}{{\left(\sqrt{E\left[{\left(x-{\mathrm{\mu }}_x\right)}^2\right]}\right)}^2}\\ =\frac{E\left[{\left(x-{\mathrm{\mu }}_x\right)}^2\right]}{E\left[{\left(x-{\mathrm{\mu }}_x\right)}^2\right]}\\ =1\end{array}$(1.14)

Figure 1.2 Density of Galton’s height data obtained from 928 (adult) children. Left panel: Density distribution of the raw data. Right panel: Density distribution of standardized scores. Vertical lines give the center (mean) of the distributions, horizontal lines indicate the spread (standard deviations) of the distributions.

Figure 1.2Long description

The left plot is labelled Children’s Height on the x-axis, ranging from 60 to 75, and Density on the y-axis, ranging from 0 to 0.15. A grey-shaded density curve illustrates the distribution of heights, with a vertical line near 68 indicating the mean, and a horizontal bar centred below this line displays standard deviation. The horizontal bar spans the distance 64 to 71. The curve that forms the boundary of the shaded area passes through the points (60, 0), (65, 0.065), (70, 0.16), and (75, 0). The right plot is labelled Standardised Children’s Height on the x-axis, ranging from minus 3 to plus 3, and Density on the y-axis peaking near 0.4. The curve that forms the boundary of the shaded area, passes through the following points (negative 3, 0), (negative 1, 0.14), (0.5, 0.4), and (3, 0). A vertical line at 0 marks the mean, and a horizontal bar underneath spans the distance negative 1 to 1. All values are estimated.

Note, again, however, that the overall distributional shape is not altered by this standardization. Therefore, higher than second moments (e.g., skewness and kurtosis discussed in detail in Sections 1.3.4 and 1.3.5) are not affected by data standardization. Section 1.3.4 moves toward higher than second moments of probability distributions and focuses on the third moment as a measure of skewness.

1.3.4 Skewness

The third standardized moment serves as a measure of the skewness and quantifies the relative size of the two tails of a probability distribution. Formally, the skewness of x (${\mathrm{\gamma }}_x$) can be expressed as

${\overline{m}}_3={\mathrm{\gamma }}_x=\frac{m_3}{m_2^{3/2}}=E\left[{\left(\frac{x-{\mathrm{\mu }}_x}{{\mathrm{\sigma }}_x}\right)}^3\right],$(1.15)

and is commonly interpreted as a measure of the symmetry of the probability distribution. The measure in Eq. (1.15) can be interpreted as the expected value of a standardized variable taken to the third power. Because the cubic function preserves the sign of the differences, one obtains (1) a value of zero for symmetric distributions, that is, cases where, relative to the mean of the distribution, the two tails are equal in size, (2) a value larger than zero for right skewed (or positively skewed) distributions, and (3) a value smaller than zero for left (negatively) skewed distributions (note, however, that special cases exist in which a skewed distribution can result in a skewness of zero; cf. Dey et al., Reference Dey, Al-Zahrani and Basloom2017). The sample estimate of the skewness is given by

${\widehat{\mathrm{\gamma }}}_x=\frac{{\widehat{m}}_3}{{\widehat{m}}_2^{3/2}}=\frac{{\sum }_i{\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)}^3/n}{{\left[{\sum }_i{\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)}^2/n\right]}^{3/2.}}$(1.16)

Note that additional multiplicative factors, resulting in the estimates ${\widehat{\mathrm{\gamma }}}_x\left[\sqrt{n\left(n-1\right)}/\left(n-2\right)\right]$ and ${\widehat{\mathrm{\gamma }}}_x{\left[\left(n-1\right)/n\right]}^{3/2}$, can be imposed on Eq. (1.16) to adjust for biases when estimating the population value (Joanes & Gill, Reference Joanes and Gill1998; Wright & Herrington, Reference Wright and Herrington2011). As n approaches infinity (i.e., n → ∞), these alternative skewness estimates approach ${\widehat{\mathrm{\gamma }}}_x$, given in Eq. (1.16).

To illustrate, consider the two artificial distributions given in Figure 1.3. Here, the left panel gives an example of a right (positively) skewed distribution exhibiting a (population) skewness of $\mathrm{\gamma }$ = 1.07. Note that, according to Eq. (1.15), data points within one standard deviation from the mean (center) of the distribution can be expected to contribute little to the overall skewness measure. The reason for this is that deviations smaller than one in absolute values (i.e., $\mid x-{\mathrm{\mu }}_x\mid$ < 1) will be even smaller when taken to the third power. In contrast, absolute deviations larger than 1 ($\mid x-{\mathrm{\mu }}_x\mid$ > 1) tend to become larger when cubed. Therefore, data points outside one standard deviation contribute more to the overall skewness than data points within the range of one standard deviation. In the case of a right skewed distribution, the right tail of the distribution tends to be longer, resulting in a skewness larger than zero.

Figure 1.3 Artificial examples of a right (positively) skewed probability distribution (left panel) and a left (negatively) skewed probability distribution (right panel). Vertical lines give the center (mean) of the distributions, horizontal lines indicate the spread (standard deviations) of the distributions.

Figure 1.3Long description

Two side-by-side area plots, both filled with grey shading under the curve. The left plot displays a right-skewed distribution with a peak toward the left side and a long tail extending to the right. A vertical line marks the mean, positioned to the right of the peak, and a horizontal bar below indicates the range of standard deviation. The right plot shows a left-skewed distribution, with a peak toward the right and a long tail on the left. A vertical line denotes the mean, positioned left of the peak, with a horizontal bar below the x-axis showing standard deviation. Both plots are unlabelled on the axes but visually emphasise the impact of skewness on the position of the mean and standard deviation.

In contrast, consider the distribution presented in the right panel of Figure 1.3. Here, the overall skewness will be smaller than zero (specifically, the distribution exhibits a (population) skewness of $\mathrm{\gamma }$ = –1.07) because, relative to the center (mean) of the distribution, the left tail of the distribution is longer than the tail to the right. Thus, more observations outside the one standard deviation interval will produce large negative values, resulting in a negative skewness estimate.

Next, we reconsider Galton’s height data with respect to its distributional symmetry. Figure 1.4 again shows the distribution of standardized children’s heights. However, this time, we focus on the tails of the distribution with respect to the one standard deviation interval. Observations outside the one standard deviation interval are marked dark gray, observations within one standard deviation are marked light gray. Clearly, the majority, that is, ~ 69% of the sample exhibit body heights within one standard deviation of the mean. The remaining 31% fall outside the one standard deviation interval. However, we now inspect the contribution of each observation within and outside one standard deviation to the estimated skewness (a similar approach has been carried out by Westfall (Reference Westfall2014) in the context of the kurtosis, discussed in Section 1.3.5). For this purpose, we start with standardizing the observations using $z_i=\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)/{\widehat{\mathrm{\sigma }}}_x$ = (x_i – 68.09)/2.52. Next, we separate the computation of the skewness for data below, within, and above one standard deviation, which results in

$\begin{array}{l}{\widehat{\mathrm{\gamma }}}_x^{\left(\mathit{\text{below}}\right)}=\frac{1}{928}\sum _{z_i\le -1}{z}_i^3=-0.749,\\ {\widehat{\mathrm{\gamma }}}_x^{\left(\mathit{\text{within}}\right)}=\frac{1}{928}\sum _{-1<z_i<1}{z}_i^3=0.019,\\ {\widehat{\mathrm{\gamma }}}_x^{\left(\mathit{\text{above}}\right)}=\frac{1}{928}\sum _{z_i\ge 1}{z}_i^3=0.643.\end{array}$(1.17)

Figure 1.4 Standardized density of Galton’s height data obtained from 928 (adult) children. Areas shaded in dark gray represent data outside one standard deviation of the mean.

Figure 1.4Long description

The x-axis is labelled Standardised Children’s Height and ranges from negative 3 to positive 3, and the y-axis is labelled Density ranging from 0 to 0.4. A curve that forms the boundary of the shaded area, passes through the following points (negative 3, 0), (negative 1, 0.14), (0.5, 0.4), and (3, 0). The curve is filled with dark grey, with a lighter grey shaded region between minus 1 and plus 1 standard deviation. A vertical line at 0 marks the mean, and a horizontal black bar beneath the curve shows the plus or minus 1 standard deviation range. The area on either side of the negative 1 and the positive 1 mark, is shaded darker than the area between them.

Clearly, observations within one standard deviation contribute little to the overall skewness, whereas observations below and above the one standard deviation interval are considerably larger. Furthermore, the contribution of observations below one standard deviation is slightly larger (in absolute values) than the contribution of observations above one standard deviation. We can, thus, expect that the left tail contributes more to the computation of the skewness than the right tail. In other words, we expect a skewness smaller than zero. A skewness estimate of $\widehat{\mathrm{\gamma }}$ = –0.088 confirms this conjecture and we can conclude that the distribution of body heights is slightly left skewed. Note, however, that the skewness estimate does not reach statistical significance (D’Agostino z = –1.10, p = 0.272). Section 1.3.5 focuses on the fourth higher moment to characterize the shape of probability distributions, the kurtosis.

1.3.5 Kurtosis

Section 1.3.4 introduced the third central moment of a probability distribution (the skewness) as a measure that quantifies the relative size of the tails of the distribution. The fourth central moment of a probability distribution, as a measure of kurtosis, can be introduced in a similar fashion. However, in contrast to the third moment (quantifying the relative tailedness around the mean), the fourth moment quantifies the total tailedness of a distribution around the mean. Formally, the kurtosis (${\mathrm{\kappa }}_x$) can be defined as

${\overline{m}}_4={\mathrm{\kappa }}_x=\frac{m_4}{m_2^2}=E\left[{\left(\frac{x-{\mathrm{\mu }}_x}{{\mathrm{\sigma }}_x}\right)}^4\right].$(1.18)

In contrast to the skewness, for decades, the interpretation of kurtosis was subject to considerable debate. Some authors suggested that kurtosis can be interpreted as a measure of the peakedness of a distribution (see, e.g., Cooligan, Reference Cooligan2013; Katz et al., Reference Katz, Elmore, Wild and Lucan2013; Lee et al., Reference Lee, Lee and Lee2013), while others suggested that kurtosis is correctly interpreted as a measure of tail heaviness (e.g., Ali, Reference Ali1974). That kurtosis is measuring both, peakedness and tail heaviness of a distribution has been suggested by, for example, Ruppert (Reference Ruppert1987) and DeCarlo (Reference DeCarlo1997). Alternative interpretations included kurtosis as a measure of bimodality (Darlington, Reference Darlington1970) and kurtosis as a measure to detect the presence of outliers (Livesey, Reference Livesey2007). Overall, however, agreement has been reached that the correct interpretation of kurtosis is that of a measure of the tail heaviness or tail extremity only, and that kurtosis does not convey any useful information regarding the peakedness of a distribution (Westfall, Reference Westfall2014).

Related to the kurtosis given in Eq. (1.18) is the so-called excess-kurtosis

${\mathrm{\delta }}_x={\mathrm{\kappa }}_x-3,$(1.19)

where the additional “– 3” accounts for the fact that the Gaussian (normal) distribution exhibits a kurtosis of 3. In other words, the excess-kurtosis measures the tail extremity of a distribution relative to the Gaussian distribution. Here, following Pearson’s (Reference Pearson1905) terminology, one can distinguish “platykurtic,” “mesokurtic,” and “leptokurtic” distributions. A platykurtic distribution exhibits a kurtosis that is smaller than the one observed for a Gaussian distribution, that is, an excess-kurtosis smaller than zero (${\mathrm{\delta }}_x$ < 0; also known as a sub-Gaussian distribution). A distribution is mesokurtic when the kurtosis equals 3 (or, equivalently, when the excess-kurtosis is zero), that is, when the tail heaviness of a distribution is similar to the one of the Gaussian distribution. Lastly, a leptokurtic distribution exhibits a kurtosis larger than 3 (or, equivalently, an excess-kurtosis larger than zero; also known as a super-Gaussian distribution). Note, however, that platykurtic and leptokurtic distributions are not directly comparable. The reason for this is that the positive excess-kurtosis values of leptokurtic distributions are unbounded (i.e., in principle, the excess-kurtosis can reach infinity). In contrast, for platykurtic distributions, the minimum possible (excess-kurtosis) value is ${\mathrm{\delta }}_x$ = –2, which is, for example, attained by the Bernoulli distribution with a probability of 50% (see, e.g., Hyvärinen et al., Reference Hyvärinen, Karhunen and Oja2001). Furthermore, it is worth noting that the skewness and the (excess-)kurtosis are not independent of each other. For arbitrary standardized distributions and relative to the skewness, the lower bound of the excess-kurtosis in relation to the skewness is given by ${\mathrm{\delta }}_x\ge {\mathrm{\gamma }}_x^2-2$ (see, e.g., Pearson, Reference Pearson1916; Teuscher & Guiard, Reference Teuscher and Guiard1995).

Next, we focus on two artificial examples of platy- and leptokurtic distributions given in Figure 1.5. The distribution in the left panel of Figure 1.5 exhibits a kurtosis of ${\mathrm{\kappa }}_x$ = 2 which is equivalent to an excess-kurtosis of ${\mathrm{\delta }}_x$ = –1. Thus, relative to the Gaussian curve, the present distribution has lighter tails; hence, the distribution is platykurtic. In contrast, the right panel of Figure 1.5 gives an example of a leptokurtic distribution with a kurtosis of ${\mathrm{\kappa }}_x$ = 9 (or an excess-kurtosis of ${\mathrm{\delta }}_x$ = 6), exhibiting heavier tails than that of a Gaussian distribution.

Figure 1.5 Artificial examples of a platykurtic (left panel) and leptokurtic (right panel) distribution. Vertical lines give the center (mean) of the distributions, horizontal lines indicate the spread (standard deviations) of the distributions.

Figure 1.5Long description

Two side-by-side density plots with grey shading under the curves. The left plot depicts a symmetric semicircular distribution with a flat base and a curved top, resembling a dome. A vertical line is centred at the peak, indicating the mean, with a horizontal bar beneath the x-axis showing the range of standard deviation. The right plot displays a sharply peaked symmetric distribution with steep slopes on either side of the narrow peak. A vertical line is aligned at the peak representing the mean, and a horizontal bar below marks the standard deviation range.

Sample estimates of the kurtosis and the excess-kurtosis are available via

$\begin{array}{l}{\widehat{\mathrm{\kappa }}}_x=\frac{{\widehat{m}}_4}{{\widehat{m}}_2^2}=\frac{{\sum }_i{\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)}^4/n}{{\left[{\sum }_i{\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)}^2/n\right]}^2}\\ {\widehat{\mathrm{\delta }}}_x={\widehat{\mathrm{\kappa }}}_x-3.\end{array}$(1.20)

Here, adjustments based on the sample size (n), e.g., using $\left({\widehat{\mathrm{\delta }}}_x+3\right){\left(\left[n-1\right]/n\right)}^2-3$, are recommended to account for biases in the population estimate (Wright & Herrington, Reference Wright and Herrington2011). However, as n → ∞, the adjusted excess-kurtosis estimate approaches ${\widehat{\mathrm{\delta }}}_x$.

To illustrate the sample computation of the (excess-)kurtosis, we, again, revisit Galton’s height data of (adult) children. Using Eq. (1.20), we obtain a kurtosis of $\widehat{\mathrm{\kappa }}$ = 2.656, which corresponds to an excess-kurtosis of $\widehat{\mathrm{\delta }}$ = 2.656 – 3 = –0.344. In other words, the height distribution is platykurtic. Furthermore, according to the Anscombe–Glynn test (Anscombe & Glynn, Reference Anscombe and Glynn1983), we obtain a z-score of –2.499 which comes with a p-value of 0.012. It is therefore concluded that the distribution significantly deviates from the Gaussian distribution with respect to its kurtosis. Next, to emphasize the interpretation as a measure of tail heaviness, we focus on the relative contributions the tails of the distribution have on the computation of its kurtosis (see, also, Westfall, Reference Westfall2014). For this purpose, we use children’s heights in standardized form, $z_i=\left(x_i-{\widehat{\mathrm{\mu }}}_x\right)/{\widehat{\mathrm{\sigma }}}_x$ = (x_i – 68.09)/2.52 (see Figure 1.4), and separate the computation of its kurtosis by data within and outside one standard deviation of the mean. Here, we obtain

$\begin{array}{l}{\widehat{\mathrm{\kappa }}}_x^{\left(\mathit{\text{within}}\right)}=\frac{1}{928}\sum _{\left|z_i\right|<1}{z}_i^4=\mathrm{0.102,}\\ {\widehat{\mathrm{\kappa }}}_x^{\left(\mathit{\text{outside}}\right)}=\frac{1}{928}\sum _{\mid z_i\mid \ge 1}{z}_i^4=2.548,\end{array}$(1.21)

suggesting that about 2.548 / (0.102 + 2.548) × 100 = 96.2% of the kurtosis statistic is determined by observations outside one standard deviation of the mean, confirming that the kurtosis measures the extremity of the tails of the distribution.

1.3.6 Higher than Fourth Moments

Sections 1.3.4 and 1.3.5 introduced third and fourth higher moments to characterize the shape of probability distributions. Both higher moments come with clear interpretations and are useful to quantify the degree of non-Gaussianity. One natural question is whether higher than fourth moments (e.g., fifth, sixth, …, moments) exist and whether these higher moments can also contribute to the characterization of the shape of probability distributions.

From the definition of the k-th standardized moment, i.e., ${\overline{m}}_k=m_k/{\mathrm{\sigma }}^k$ with $m_k=E\left[{\left(x-{\mathrm{\mu }}_x\right)}^k\right]$, it follows that standardized moments can be taken to any higher powers (e.g., k = 5, 6, 7, …). Here, odd-powered standardized moments (e.g., using k = 5, 7, …) quantify the relative tailedness, whereas even-powered standardized moments (e.g., using k = 6, 8, …) quantify the total tailedness of a probability distribution. The fifth and sixth standardized moments, for example, are referred to as the hyper-skewness and the hyper-kurtosis (Bazavov et al., Reference Bazavov, Bollweg, Ding, Enns, Goswami, Hegde, Kaczmarek, Karsch, Larsen, Mukherjee, Ohno, Petreczky, Schmidt, Sharma and Steinbrecher2020), suggesting that third and fifth as well as fourth and sixth standardized moments quantify similar characteristics of a probability distribution. Such similarities exist for any even and odd higher powers. To visualize these similarities, consider Figures 1.6 and 1.7 (see also Gundersen, Reference Gundersen2020).

Figure 1.6 Artificial right skewed distribution with the functions f(x) = x³ and f(x) = x⁴ superimposed.

Figure 1.6Long description

The graph plots two mathematical functions and a shaded density area. The solid curve represents the function f of x equals x cubed, which is antisymmetric about the origin and crosses through it, with negative values to the left and positive values to the right. The dashed curve shows f of x equals x to the power 4, forming a symmetric parabolic shape with its minimum at the origin and increasing steeply on both sides. A smooth shaded area is plotted between these two curves, concentrated around the interval from minus 0.5 to 1 along the x-axis, peaking slightly above 1 on the y-axis. A horizontal line is drawn at y equals 0.

Figure 1.7 Artificial right skewed distribution with the functions f(x, k) = x^k (k = 5, 6, 7, 8) superimposed.

Figure 1.7Long description

The graph displays the curves of four functions as follows. f of x equals x superscript 5 end superscript, f of x equals x superscript 6 end superscript, f of x equals x superscript 7 end superscript, and f of x equals x superscript 8 end superscript, each represented by a distinct line style, solid, dotted, dashed, and dash-dot, respectively. The x-axis ranges from minus 1.5 to 1.5, and the y-axis ranges from minus 1.5 to 2. A shaded region lies between the curves of f of x equals x superscript 5 end superscript and f of x equals x superscript 8 end superscript in the interval around x equals 0. The functions are plotted symmetrically around the origin.

Figure 1.6 shows an artificial positively skewed distribution with the functions f(x) = x³ and f(x) = x⁴ superimposed. In line with the properties of skewness and kurtosis described in Sections 1.3.4 and 1.3.5, values close to the center of the distribution can be expected to contribute little to the statistics, however, for larger negative deviations from the center, large negative values are observed for f(x) = x³ and large positive values are obtained for f(x) = x⁴. In contrast, for large positive deviations from the center, both f(x) = x³ and f(x) = x⁴ attain large positive values.

Next, consider Figure 1.7 which gives the same probability distribution with functions involving higher than fourth powers superimposed. Specifically, f(x, k) = x^k with k = 5, 6, 7, and 8. Clearly, the selected odd powers (k = 5, 7) behave in a way similar to raising x to the third power, and even powers (k = 6, 8) behave in a way similar to raising x to the fourth power. In other words, higher than fourth standardized moments, in essence, replicate the information gathered by third and fourth standardized moments. Since it is easier to reliably estimate lower moments than higher moments from empirical data, third and fourth moments are preferred over higher than fourth moments.

1.3.7 Cumulants

So far, this chapter has focused on moments of variables to characterize their distributions. Moments, however, are not the only constants that can be used for this purpose. Cumulants constitute an alternative set of constants to characterize distributional features of variables, which are, from a theoretical perspective, sometimes more useful than moments (Chapter 4 will make use of cumulants to derive direction of dependence statistics based on error distributions of causally competing linear models). Formally, cumulants describe the coefficients of the Taylor series expansion of the second characteristic function (or cumulant generating function). The second characteristic function is given by the natural logarithm of the first characteristic function (or moment generating function) which itself is defined as the continuous Fourier transform of the probability density function $p_x\left(x\right)$ (Hyvärinen et al., Reference Hyvärinen, Karhunen and Oja2001; Stuart & Ord, Reference Stuart and Ord1994).

Let $\mathit{cu}{m}_k\left(x\right)$ be the k-th cumulant of x. When the expected value of x is nonzero, i.e., E[x] ≠ 0, the first four cumulants can be expressed as (see, e.g., Hyvärinen et al., Reference Hyvärinen, Karhunen and Oja2001)

$\begin{array}{l}\mathit{cu}{m}_1\left(x\right)=E\left[x\right]\\ \mathit{cu}{m}_2\left(x\right)=E\left[x^2\right]-E{\left[x\right]}^2\\ \mathit{cu}{m}_3\left(x\right)=E\left[x^3\right]-3E\left[x^2\right]E\left[x\right]+2E{\left[x\right]}^3\\ \mathit{cu}{m}_4\left(x\right)=E\left[x^4\right]-3E{\left[x^2\right]}^2-4E\left[x^3\right]E\left[x\right]+12E\left[x^2\right]E{\left[x\right]}^2-6E{\left[x\right]}^4\end{array}$(1.22)

which, for a zero-mean variable (i.e., E[x] = 0) reduce to

$\begin{array}{l}\mathit{cu}{m}_1\left(x\right)=0\\ \mathit{cu}{m}_2\left(x\right)=E\left[x^2\right]\\ \mathit{cu}{m}_3\left(x\right)=E\left[x^3\right]\\ \mathit{cu}{m}_4\left(x\right)=E\left[x^4\right]-3E{\left[x^2\right]}^2.\end{array}$(1.23)

In other words, a simple correspondence between cumulants and central moments exists for the first three moments. That is, the first three cumulants equal the first three central moments. However, this simple relation no longer holds for the fourth moment (see also Dodge & Rousson, Reference Dodge and Rousson1999).

As mentioned at the beginning of this section, cumulants have several advantageous properties. First, if x is a random variable and c is an additive constant, then the cumulants of x + c are given by

$\mathit{cu}{m}_k\left(x+c\right)=\mathit{cu}{m}_k\left(x\right)$(1.24)

for k > 1. For k = 1, one obtains $\mathit{cu}{m}_1\left(x+c\right)=\mathit{cu}{m}_1\left(x\right)+c$. Second, if c is a multiplicative constant, then

$\mathit{cu}{m}_k\left(\mathit{cx}\right)=c^k\mathit{cu}{m}_k\left(x\right),$(1.25)

Third, for p independent random variables $x_1$, $x_2$, …, $x_p$, it holds that

$\mathit{cu}{m}_k\left(x_1+x_2+\dots +x_p\right)=\mathit{cu}{m}_k\left(x_1\right)+\mathit{cu}{m}_k\left(x_2\right)+\dots +\mathit{cu}{m}_k\left(x_p\right).$(1.26)

In words, the k-th cumulant of the sum of independent variables equals the sum of the k-th cumulants of the summands. Note that this additivity property also holds for the first three moments, but does not hold for the fourth moment (Dodge & Rousson, Reference Dodge and Rousson1999).

The three properties listed in Eqs. (1.24), (1.25), and (1.26) have immediate consequences when working with cumulants in the context of linear models. To illustrate, consider the simple linear function $y=a+\mathit{bx}+z$ with x, y, and z being continuous random variables, and $a$ and $b$ denoting constants. Making use of Eqs. (1.24), (1.25), and (1.26) implies that the k-th cumulant of y (with k > 1) is given by

$\begin{array}{l}\mathit{cu}{m}_k\left(y\right)=\mathit{cu}{m}_k\left(a+\mathit{bx}+z\right)\\ =\mathit{cu}{m}_k\left(\mathit{bx}\right)+\mathit{cu}{m}_k\left(z\right)\\ =b^k\mathit{cu}{m}_k\left(x\right)+\mathit{cu}{m}_k\left(z\right),\end{array}$(1.27)

that is, the k-th cumulant of y equals the sum of the weighted k-th cumulant of x and the k-th cumulant of z. These properties will play a vital role in the development of distributional measures of direction of dependence.

1.3.8 Joint Moments

So far, this chapter has focused on moments and cumulants of single random variables. Both sets of quantities, however, can be extended to the multivariate case. For simplicity, we focus on the bivariate case with the understanding that the presented concepts can be extended to more than two variables. Let x and y be two random variates. The bivariate moment $m_{\mathit{kl}}$ about the origins ${\mathrm{\mu }}_x$ and ${\mathrm{\mu }}_y$ is defined as (cf. Stuart & Ord, Reference Stuart and Ord1994)

$m_{\mathit{kl}}=E\left[{\left(x-{\mathrm{\mu }}_x\right)}^k{\left(y-{\mathrm{\mu }}_y\right)}^l\right]={\int }_{\mathit{-}\mathit{\infty }}^{\mathit{\infty }}{\int }_{\mathit{-}\mathit{\infty }}^{\mathit{\infty }}{\left(x-{\mathrm{\mu }}_x\right)}^k{\left(y-{\mathrm{\mu }}_y\right)}^l\mathit{dF},$(1.28)

where $m_{k0}$ denotes the k-th moment of the marginal distribution of x and $m_{0l}$ denotes the l-th moment of the marginal distribution of y. When k ≠ 0 and l ≠ 0, $m_{\mathit{kl}}$ is known as the product moment of x and y. For example, when k = l = 1, the quantity m₁₁ refers to the covariance of x and y,

$m_{11}=cov\left(x,y\right)=E\left[\left(x-{\mathrm{\mu }}_x\right)\left(y-{\mathrm{\mu }}_y\right)\right].$(1.29)

When x and y have zero means and unit variances (which can be achieved through variable standardization), k = l = 1 leads to the Pearson correlation ${\mathrm{\rho }}_{\mathit{xy}}$ (for further details see Chapter 2),

$cov\left(x,y\right)=E\left[\mathit{xy}\right]={\mathrm{\rho }}_{\mathit{xy}}.$(1.30)

In a similar way, we can define higher-order cross-products, sometimes called higher-order covariances (Dodge & Rousson, Reference Dodge and Rousson2001; Wiedermann et al., Reference Wiedermann, Wiedermann, Kim, Sungur and von Eye2021),

$m_{\mathit{kl}}=cov{\left(x,y\right)}_{\mathit{kl}}=E\left[{\left(x-{\mathrm{\mu }}_x\right)}^k{\left(y-{\mathrm{\mu }}_y\right)}^l\right],$(1.31)

and higher-order correlations,

$\mathrm{cor}{\left(x,y\right)}_{\mathit{kl}}=\frac{cov{\left(x,y\right)}_{\mathit{kl}}}{{\mathrm{\sigma }}_x^k{\mathrm{\sigma }}_y^l}.$(1.32)

Using Eqs. (1.31) and (1.32) one can define higher-order co-moments. For example, using {k, l} = {1, 2} and {2, 1} one obtains two estimates of the co-skewness of x and y,

$\begin{array}{l}\mathrm{cor}{\left(x,y\right)}_{12}=\frac{cov{\left(x,y\right)}_{12}}{{\mathrm{\sigma }}_x^1{\mathrm{\sigma }}_y^2}\\ \mathrm{cor}{\left(x,y\right)}_{21}=\frac{cov{\left(x,y\right)}_{21}}{{\mathrm{\sigma }}_x^2{\mathrm{\sigma }}_y^1}.\end{array}$(1.33)

In a similar way, measures of the co-kurtosis of x and y can be defined when {k, l} = {1, 3}, {2, 2}, and {3, 1},

$\begin{array}{l}\mathrm{cor}{\left(x,y\right)}_{13}=\frac{cov{\left(x,y\right)}_{13}}{{\mathrm{\sigma }}_x^1{\mathrm{\sigma }}_y^3}\\ \mathrm{cor}{\left(x,y\right)}_{22}=\frac{cov{\left(x,y\right)}_{22}}{{\mathrm{\sigma }}_x^2{\mathrm{\sigma }}_y^2}\\ \mathrm{cor}{\left(x,y\right)}_{31}=\frac{cov{\left(x,y\right)}_{31}}{{\mathrm{\sigma }}_x^3{\mathrm{\sigma }}_y^1}.\end{array}$(1.34)

Co-skewness and co-kurtosis are natural extensions of the covariance to higher than second moments. Here, the co-skewness is a measure of distributional asymmetry that preserves the association between the two variables. Similarly, co-kurtosis is a measure of tailedness that preserves the x–y relation. Chapter 3 shows that these co-moment measures exhibit properties that are useful in determining the direction of causation in linear models.

To illustrate the application of higher co-moment measures, we, once again, revisit Galton’s height data. However, this time, we consider the heights of both (adult) children and their parents. In Galton’s (Reference Galton1886) original study, parents’ heights were averaged (forming midparent’s heights), with female heights being multiplied by a factor of 1.08 (for further details see Hanley, Reference Hanley2004). Figure 1.8 shows the scatterplot of the standardized children’s and midparent’s heights (data points were slightly jittered to avoid overlap). Here, the LOWESS fit suggests a slightly nonlinear relationship (see also Wachsmuth et al., Reference Wachsmuth, Wilkinson and Dallal2003).

Figure 1.8 Scatterplot of standardized Galton’s height data with LOWESS line superimposed.

Figure 1.8Long description

The scatter plot displays individual data points representing the relationship between standardised midparental height on the x-axis and standardised children's height on the y-axis, with both axes ranging approximately from minus 2.5 to 2.5. Each black circle denotes an observation. A smooth fitted curve passes through the cloud of points, indicating a positive non-linear trend between the two variables. The data points are spread with considerable vertical variation for each horizontal value. The best fit curve passes through the following points. (negative 2, negative 1), (negative 1, negative 0.5), (0, 0), (1, 0.5), and (2, 1). All values are estimated.

While Galton’s data have repeatedly been used to illustrate linear regression and correlation, the present re-analysis focuses on higher moment characteristics of the data using the higher-order correlation formula in Eq. (1.32). Recall that $\mathrm{cor}{\left(x,y\right)}_{11}$ leads to the standard Pearson correlation coefficient, and $\mathrm{cor}{\left(x,y\right)}_{k0}$ and $\mathrm{cor}{\left(x,y\right)}_{0l}$ can be used to describe higher moments of the marginal distributions. Table 1.1 summarizes the higher-order correlation estimates for several sets of power values. In addition, 95% bootstrap percentile confidence limits (based on 2,000 resamples) are reported for statistical inference.

Table 1.1Higher-order correlation estimates $\mathit{cor}{\left(x,y\right)}_{\mathit{kl}}$ for different power values (x = children’s height, y = mid-parent’s height) together with the 95% percentile bootstrap confidence intervals (based on 2,000 resamples).

Table 1.1Long description

The column-headers are labelled: K, l, estimate, and 95% C I. The last column is divided further into lower and upper. The statistics are labelled: correlation, skewness of child, skewness of mid-parent, co-skewness, kurtosis of child, kurtosis of mid-parent, and co-kurtosis. The listed values are as follows. Row number 1: Correlation, 1, 1, 0.459, 0.406, 0.513. Row number 2: Skewness of Child, 3, 0, negative 0.088, negative 0.189, 0.025. Row number 3: Skewness of Mid-Parent, 0, 3, negative 0.035, negative 2.8E minus 17, 2.9E minus 17. Row number 4: Co-Skewness, 1, 2, 0.074, negative 0.029, 0.176. Row number 5: Co-Skewness, 2, 1, 0.095, negative 0.001, 0.189. Row number 6: Kurtosis of Child, 4, 0, 2.656, 2.507, 2.825. Row number 7: Kurtosis of Mid-Parent, 0, 4, 3.058, 2.873, 3.270. Row number 8: Co-Kurtosis, 1, 3, 1.537, 1.333, 1.722. Row number 9: Co-Kurtosis, 3, 1, 1.298, 1.113, 1.485. Row number 10: Co-Kurtosis, 2, 2, 1.516, 1.367, 1.660.

First, using k = l = 1 gives a Pearson correlation estimate of 0.459 (95% CI = [0.406, 0.513]), suggesting that children’s and midparent’s heights are positively correlated. Next, to estimate the skewness of the marginal height distributions, we use k = 3 and l = 0 as well as k = 0 and l = 3, leading to skewness estimates of children’s and midparent’s heights. Both estimates are close to zero, suggesting that both distributions are symmetric. In addition, the two co-skewness estimates (i.e., using k = 2 and l = 1 as well k = 1 and l = 2) are close to zero as well. The last step focuses on the tailedness of the height distributions and estimates the kurtosis of the marginal distributions (with {k, l} = {4, 0} and {0, 4}) and the three co-kurtosis measures (with {k, l} = {3, 1}, {2, 2}, and {1, 3}). Results suggest that children’s heights tend to be platykurtic (with a kurtosis value significantly smaller than 3, which replicates results presented above) and midparent’s heights tend to be mesokurtic (with a kurtosis estimate close to 3). Finally, to interpret results of the three co-kurtosis estimates, it is helpful to consider the theoretical values one would expect when both variables were Gaussian. In this case, the co-kurtosis reduces to

$\begin{array}{l}\mathrm{cor}{\left(x,y\right)}_{22}=1+2{\mathrm{\rho }}_{\mathit{xy}}^2\\ \mathrm{cor}{\left(x,y\right)}_{13}=\mathrm{cor}{\left(x,y\right)}_{31}=3{\mathrm{\rho }}_{\mathit{xy}}.\end{array}$(1.35)

In other words, under the assumption that children’s and midparent’s heights are Gaussian, one would expect $\mathrm{cor}{\left(x,y\right)}_{22}=1+2\bullet {0.459}^2$ = 1.421 and $\mathrm{cor}{\left(x,y\right)}_{13}=\mathrm{cor}{\left(x,y\right)}_{31}=3\bullet 0.459$ = 1.377. Considering the point estimate together with the 95% CIs suggests that $\mathrm{cor}{\left(x,y\right)}_{22}$ is not significantly larger with $\mathrm{cor}{\left(x,y\right)}_{22}$ = 1.516 (95% CI = [1.367, 1.660]). In addition, 95% CIs of $\mathrm{cor}{\left(x,y\right)}_{13}$ and $\mathrm{cor}{\left(x,y\right)}_{31}$ include the theoretical value of 1.377. While, in the present application, we used estimates of higher-order correlations to characterize marginal and joint height distributions, Chapters 3 and 4 will make use of higher-order correlations to construct test statistics that can be informative in discerning the direction of causation.