A Note on the Solution Multiplicity of the Vale–Maurelli Intermediate Correlation Equation

Abstract

The Vale and Maurelli algorithm is a widely used method that allows researchers to generate multivariate, nonnormal data with user-specified levels of skewness, excess kurtosis, and a correlation structure. Before obtaining the desired correlation structure, a transitional step requires the user to calculate the roots of a cubic polynomial referred to as the intermediate correlation equation. The Cardano method and a corollary of Rouché’s theorem were used to derive closed-form solutions to this equation. These solutions highlight the fact that three real-valued roots are possible, and solutions can either contain multiple roots within the allowable correlation range or can all be greater than 1 for many combinations of skewness and excess kurtosis. It was also found that large values of excess kurtosis exacerbate the issue of multiple solutions or solutions greater than 1, bringing further restrictions into the combinations of higher order moments that researchers can simulate from. A small computer study on the power of the t test for the correlation coefficient also uncovered that different values of the intermediate correlation can influence the results from Monte Carlos simulations. This note is intended to inform both researchers who conduct simulations with nonnormal data and users who inform their data analysis practice from simulation studies.

Keywords

computer applications computerized testing correlational analysis dimensionality factor analysis hierarchical linear modeling item response theory psychometrics research methodology validity/reliability

1. The Third-Order Polynomial Transformation

Within the social sciences, few algorithms have seen the level of popularity of the third-order polynomial method proposed by Fleishman (1978) in order to generate nonnormal data for Monte Carlo simulation studies. If a researcher is interested in generating nonnormal variables with a prespecified level of skewness $(γ_{1})$ or kurtosis¹ $(γ_{2})$ they may define a continuous random variable Y as:

Y = a + b Z + c Z^{2} + d Z^{3},

where a, b, c, and d are real-valued polynomial coefficients and $Z \sim N (0, 1)$ . In order to find the values of the coefficients needed to specify the degree of nonnormality, Fleishman (1978) derived the following system of polynomial equations:

a = - c .

b^{2} + 6 b d + 2 c^{2} + 15 d^{2} = 1.

2 c (b^{2} + 24 b d + 105 d^{2} + 2) = γ_{1} .

24 (b d + c^{2} [1 + b^{2} + 28 b d] + d^{2} [12 + 48 b d + 141 c^{2} + 225 d^{2}]) = γ_{2},

so that for values of $(γ_{1}, γ_{2})$ , one only needs to find the roots of the system (Equations 2 –5) and substitute the solutions back in Equation 1 to obtain a random variable with the desired skewness and kurtosis values.

The Fleishman method was originally derived to generate univariate distributions, and modifications are needed to generalize it into a multivariate setting. If variables following a specific covariance structure are Fleishman-transformed, the final nonnormal variables will have different correlations from those originally intended. If Fleishman-transformed values are rotated to follow a covariance structure by some procedures such as the Kaiser and Dickman (1962) or Cholesky decomposition, the higher order moments will not match the initial ones. A solution to this problem was proposed by Vale and Maurelli (1983), where they conceived a second correlation matrix, referred to as the “intermediate correlation matrix,” which would impose an initial covariance structure on the data. Once the variables were correlated, each one would be individually subjected to the third-order polynomial transformation, so that the intermediate correlation matrix would be modified into the correlation structure initially intended by the researcher. Consider again a random variable Y, which has been Fleishman-transformed from a standard normal variable defined as above, $Z \sim N (0, 1)$ . Y can be expressed as

Y = w^{'} z = [a, b, c, d] [\begin{matrix} 1 \\ \begin{matrix} Z \\ \begin{matrix} Z^{2} \\ Z^{3} \end{matrix} \end{matrix} \end{matrix}] .

Since Y is standardized by construction, the cross product of any two variables $(Y_{1}, Y_{2})$ that have underwent the third-order polynomial transformation is their correlation coefficient. This correlation can be expressed as

r_{Y_{1} Y_{2}} = E [Y_{1} Y_{2}] = E [w_{1}^{'} z_{1} z_{2}^{'} w_{2}] = w_{1}^{'} E [z_{1} z_{2}^{'}] w_{2} = w_{1}^{'} R w_{2},

where the expectation of R is given by

R = [\begin{matrix} 1 & 0 & 1 & 0 \\ 0 & r_{Z_{1} Z_{2}} & 0 & 3 r_{Z_{1} Z_{2}} \\ 1 & 0 & 2 r_{Z_{1} Z_{2}}^{2} + 1 & 0 \\ 0 & 3 r_{Z_{1} Z_{2}} & 0 & 6 r_{Z_{1} Z_{2}}^{3} + 9 r_{Z_{1} Z_{2}}^{} \end{matrix}] .

An algebraic rearrangement of Equations 7 and 8 yields the final intermediate correlation formula as conceived by Vale and Maurelli (1983):

r_{Y_{1} Y_{2}} = r_{Z_{1} Z_{2}} (b_{1} b_{2} + 3 b_{1} d_{2} + 3 d_{1} b_{2} + 9 d_{1} d_{2}) + r_{Z_{1} Z_{2}}^{2} (2 c_{1} c_{2}) + r_{Z_{1} Z_{2}}^{3} (6 d_{1} d_{2}) .

2. Limitations of the Vale–Maurelli Transformation

The algorithm suggested by Vale and Maurelli (1983) exhibits certain limitations, many of which stem from the way in which the Fleishman (1978) method operates. Given that the focus of this note is the Vale–Maurelli approach, the interested reader can refer to Astivia and Zumbo (2015), Headrick (2002), and Headrick (2004) for an overview of the limits to the Fleishman (1978) method.

The first limitation of the Vale–Maurelli method is the restriction concerning the types of multivariate nonnormal data that it can generate. Foldnes and Grønneberg (2015) derived the probability density of Vale–Maurelli-transformed variables and found them to be contained within the family of Gaussian copula distributions. This brings into question whether or not the Monte Carlo studies that have relied on this method are simulating nonnormal data that are truly encountered in the real world or if they are restricted to a specific subset of multivariate distributions. A second limitation of this method pertains to the manner in which the intermediate correlation matrix is assembled. Equation 9 implies that the correlations must be calculated in a pairwise fashion, and there is no guarantee that this matrix will be positive definite once compiled. For such cases, Demirtas, Hedeker, and Mermelstein (2012) recommend calculating the “nearest matrix” as outlined in Higham (2002) and utilize said matrix instead.

Another limitation is that for certain values of skewness and kurtosis, the Vale–Maurelli transformation may not span the full [−1, +1] correlation range, limiting the effect size values that researchers can choose from. Headrick (2002, 2010) has shown that the type of coefficients obtained from the Fleishman (1978) method may place undue boundaries on the solution to Equation 9 and a higher order transformation, such as the fifth-order polynomial he derived, would be needed to alleviate this issue.

3. Multiple Solutions to the Vale–Maurelli Transformation: A Mathematical Description

An important feature of the Vale–Maurelli algorithm that has gone mostly unnoticed within the published literature is that the intermediate correlation shown in Equation 9 is defined as a third-degree polynomial, which, as per the Fundamental Theorem of Algebra, must have three roots (including those with multiplicities greater than 1). There is currently no description of these solutions, and there are no guidelines concerning which ones should be preferred over others when designing a Monte Carlo study. The purpose of this note is, therefore, to provide a full mathematical account of the multiple solutions to the Vale–Maurelli intermediate correlation equation, explore their properties, and offer recommendations to researchers who may be interested in using this method to generate multivariate, nonnormal data for simulation purposes.

Define $(Y_{1}, Y_{2})$ as two continuous, random variables that have been Fleishman-transformed from two standard normal variables $(Z_{1}, Z_{2})$ , respectively. Let $z : = r_{Z_{1} Z_{2}}, q : = - r_{Y_{1} Y_{2}}, δ : = (b_{1} + 3 d_{1}) (b_{2} + 3 d_{2}), β : = 2 c_{1} c_{2}, α : = 6 d_{1} d_{2}$ , so that Equation 9 simplifies to

α z^{3} + β z^{2} + δ z + q = 0.

Note that since q is chosen by the researcher, it is always a predetermined constant.

Using the Cardano method (see Mollin, 1947, for a full description), it is possible to find closed-form expressions for the cubic polynomial roots. For this specific instance, the roots are

z_{1} = S + T - \frac{β}{3 α},

z_{2} = - \frac{S + T}{2} - \frac{β}{3 α} + \frac{i \sqrt{3}}{2} (S - T),

z_{3} = - \frac{S + T}{2} - \frac{β}{3 α} - \frac{i \sqrt{3}}{2} (S - T),

where

S = \sqrt[3]{R + \sqrt{Q^{3} + R^{2}}},

T = \sqrt[3]{R - \sqrt{Q^{3} + R^{2}}},

with

Q = \frac{3 α δ - β^{2}}{9 α^{2}}, R = \frac{9 α β δ - 27 α^{2} q - 2 β^{3}}{54 α^{3}} .

The expression $Q^{3} + R^{2}$ is referred to as the discriminant, and it describes the nature of the solutions as falling in one of the three possible cases.

If $Q^{3} + R^{2} = 0$ :

Q^{3} + R^{2} = 0 \underset{}{\Rightarrow} S = T \underset{}{\Rightarrow} \frac{i \sqrt{3}}{2} (S - T) = 0,

so that $z_{1}, z_{2}, and z_{3}$ are all real-valued numbers, with $z_{2} = z_{3}$ .

If $Q^{3} + R^{2} > 0 :$

\begin{array}{l} Q^{3} + R^{2} > 0 \underset{}{\Rightarrow} \sqrt{Q^{3} + R^{2}} \in ℝ, \\ \underset{}{\Rightarrow} S \neq T, S, T \in ℝ, \\ \underset{}{\Rightarrow} \frac{i \sqrt{3}}{2} (S - T) \neq 0, \frac{i \sqrt{3}}{2} (S - T) \in ℂ, \end{array}

so that z₁is the only real-valued root, with $z_{2}, z_{3}$ being complex conjugates of one another.

If $Q^{3} + R^{2} < 0 :$

\begin{array}{l} Q^{3} + R^{2} < 0 \underset{}{\Rightarrow} \sqrt{Q^{3} + R^{2}} \in ℂ, \\ \underset{}{\Rightarrow} S \neq T, S, T \in ℂ, \\ \underset{}{\Rightarrow} \frac{i \sqrt{3}}{2} (S - T) \neq 0, \frac{i \sqrt{3}}{2} (S - T) \in ℝ, \end{array}

which implies $z_{1}, z_{2}, and z_{3}$ are all real-valued roots and all three are different.

A back substitution to the original polynomial coefficients from Equation 9 results in:

R = \frac{9 α β δ - 27 α^{2} q - 2 β^{3}}{54 α^{3}} = \frac{54 c_{1} c_{2} d_{1} d_{2} (b_{1} + 3 d_{1}) (b_{2} + 3 d_{2}) - 486 d_{1}^{2} d_{2}^{2} r_{Y_{1} Y_{2}} - 8 c_{1}^{3} c_{2}^{3}}{5832 d_{1}^{3} d_{2}^{3}} .

Q = \frac{3 α δ - β^{2}}{9 α^{2}} = \frac{9 d_{1} d_{2} (b_{1} + 3 d_{1}) (b_{2} + 3 d_{2}) - 2 c_{1}^{2} c_{2}^{2}}{162 d_{1}^{2} d_{2}^{2}} .

The analytical form of the roots obtained via the Cardano method makes evident some of the conditions required to be able to successfully use the Vale–Maurelli transformation. The first condition is that the roots must be real-valued numbers. For the case where the discriminant is greater than 0, choosing the real solution would be the only sensible path for researchers wishing to use the intermediate correlation equation in their simulation methodology. For the other two cases, a second reasonable requirement would be to only select the solution (or solutions) that exist within the range of the Pearson correlation coefficient [−1, +1]. Determining which intermediate correlations are “valid” is equivalent to finding boundary conditions of the solutions to the cubic polynomial described in Equation 9 outside of the interval [−1, +1].

Consider the polynomial $p (x) = a_{n} x^{n} + a_{n - 1} x^{n - 1} + a_{n - 2} x^{n - 2} + . . . + a_{k} x^{k} + . . . + a x + a_{0}$ , where all the coefficients a_i are real numbers. A corollary of Rouché’s theorem (cf. Saff & Snider, 2003, for a detailed description) states that if the polynomial coefficient index $k = n$ (for the upper bound) or $k = 0$ (for the lower bound), then there always exist a positive real number w such that

| a_{k} | w^{k} \geq | a_{0} | + | a | w + . . . + | a_{k - 1} | w^{k - 1} + | a_{k + 1} | w^{k + 1} + . . . + | a_{n} | w^{n} .

The constant w takes the following forms depending on whether it is an upper or a lower bound. For the upper bound,

w = max (1, \frac{1}{| a_{n} |} (| a_{0} | + | a_{1} | + | a_{2} | + . . . + | a_{n - 1} |)),

and for the lower bound,

w = \frac{| a_{0} |}{max (| a_{0} |, | a_{1} | + {| a_{2} |}_{.} + . . . + | a_{n - 1} | + | a_{n} |)} .

Taking the case of the lower bound, by substituting the coefficients of Equation 9, one obtains

w = \frac{| r_{Y_{1} Y_{2}} |}{max (| r_{Y_{1} Y_{2}} |, | (b_{1} + 3 d_{1}) (b_{2} + 3 d_{2}) | + | 2 c_{1} c_{2} | + | 6 d_{1} d_{2} |)},

which implies that if $| r_{Y_{1} Y_{2}} | \geq | (b_{1} + 3 d_{1}) (b_{2} + 3 d_{2}) | + | 2 c_{1} c_{2} | + | 6 d_{1} d_{2} |$ , then the smallest possible solution (in absolute value) to Equation 9 is 1, and the other two solutions must be greater than 1 so that a valid intermediate correlation does not exist.

4. Multiple Solutions to the Vale–Maurelli Transformation: A Computer Simulation Exploration

The analytical description of the intermediate correlation equation discussed in the previous section shines a light on two important special cases. The first is that more than one real-valued solution is possible for Equation 9. Although it appears that researchers are more likely to find themselves in the instance where the discriminant $Q^{3} + R^{2} > 0$ , because both R² and the denominator of Q³ are always positive, there is nothing within the solutions to the system described in Equations 2 through 5 that would guarantee this. The second is that, for the cases where all three solutions are real numbers, multiple solutions can either exist within the [−1, +1] range or they or can all be strictly greater than 1.

To the authors’ knowledge, there are no investigations within the published literature that describe the prevalence of these cases or that attempt to link them to any specific values of skewness and kurtosis. These properties of the Vale–Maurelli transformation can have several unintended consequences that could influence the results and conclusions of Monte Carlo simulations that rely on this algorithm. To begin with, they place practical restrictions on the combinations of higher order moments that researchers can simulate from beyond the quadratic relationship that limits the third and fourth moments. If certain values of skewness and kurtosis yield solutions to Equation 9 that are strictly greater than 1, then these values cannot be used in a simulation since the intermediate correlation matrix would be nonsensical. For the cases where more than one value within the valid correlation range is possible, it is becomes unclear which solution among the possible candidates should be selected. Take, for instance, the field of structural equation modeling (SEM) that has relied considerably on the Fleishman (1978) method and the Vale–Maurelli transformation to explore the robustness of various estimators, such as normal theory maximum likelihood (e.g., Curran, West, & Finch, 1996; Falk, 2017; Flora & Curran, 2004). Foldnes and Grønneberg (2017) derived closed-form expressions of the population Γ matrix of fourth-order moments (which is used in almost all robust corrections for SEM estimators) for Vale–Maurelli-transformed variables and discovered that the values of this Γ matrix depend on the intermediate correlations. If different solutions to Equation 9 imply different population Γ matrices, it becomes unclear whether or not certain solutions would cause this Γ matrix to behave in different and unexpected ways.

Unfortunately, it is impossible to find closed-form expressions for the skewness/kurtosis values that would yield either multiple valid solutions or solutions greater than 1 due to the fact that algebraically solving for the system (Equations 2 –5) yield fifth- and sixth-degree polynomials for which a closed-form solution cannot be found as per the Abel–Ruffini Theorem. To tackle this issue, an empirical computer investigation was conducted in the statistical programming language R (R Core Team, 2017).

A grid-search method was used over a large collection of possible combinations of skewness/kurtosis values and researcher-specified final correlation values of 0.1 and 0.5, which correspond to a small and large effect size (Cohen, 1992). The first step was generating a potential list of values, one for skewness and another one for the kurtosis, from −2 to 101.5 in intervals of 0.005 to mimic a large range of options for the third and fourth moments. The values of these lists were paired with one another to ensure every possible combination between the two were present, creating a long series of potential candidate values for the $(γ_{1}, γ_{2})$ pairs associated with each $(Y_{1}, Y_{2})$ . Once the list of potential skewness and kurtosis combinations was built, it needed to be further refined so only sensible values would remain. Kraatz (2011) found that the Fleishman (1978) method could only generate real-valued solutions to the system of polynomial Equations 2 through 5, whenever $γ_{2} \geq 1.588, γ_{1}^{2} - 1.139$ and if kurtosis values were less than approximately 101.38. All values that did not satisfy these two conditions were removed from the list.

Fleishman polynomial coefficients were obtained for the $(γ_{1}, γ_{2})$ in the reduced list and substituted in the intermediate correlation equation, isolating all solutions which were either greater than |1| (meaning all the real-valued roots where |1| or greater so researchers could not choose any other root) or if multiple solutions existed within the allowable range of the correlation coefficient.

Figure 1 shows the positive section of the skewness/kurtosis parabola defined by the Fleishman (1978) method. For pairs of variables $(Y_{1}$ , $Y_{2})$ , the intermediate correlation is the same irrespective of whether Y₁ or Y₂ have specific values of skewness or kurtosis within any given combination of the two, which is why only one set of (skewness, kurtosis) is presented in the figure. For example, if Y₁ has a skewness of 2 and a kurtosis of 7 and Y₂ has a skewness of 3 and a kurtosis of 21, the same intermediate correlation would be found if Y₁ has a skewness of 3 and a kurtosis of 21 and Y₂ has a skewness of 2 and a kurtosis of 7. Within Figure 1 the gray circles represent sets of unique (skewness, kurtosis) set of Y₁ and the black triangles are the unique (skewness, kurtosis) of Y₂, which, when chosen together to solve for the Fleishman coefficients, result in an intermediate correlation outside of the allowable [−1, +1] range.

FIGURE 1.

Combinations of skewness and kurtosis values for final correlations of 0.1 and 0.5 that yield solutions to the intermediate correlations outside of the [−1, +1] range.

The first noticeable difference between both panels is that the number of instances where the solutions are greater than |1| is considerably larger for the effect size of 0.5 than for 0.1. In both cases, it appears that the combination of low skewness and high kurtosis values results in intermediate correlations that violate the allowable range, although for the case of r = .5, the differences between the kurtosis values of Y₁ and Y₂ need to be substantially larger. For instance, the first pairs of $(γ_{1}, γ_{2})$ represented in the left panel of Figure 1 resulting in an intermediate correlation outside the allowable range are $(1, 1)$ for Y₁ and $(4.685, 25.110)$ for Y₂, therefore a difference of almost 25 points between one kurtosis and the other. The minimal overlap between the gray and black areas on the right-most panel of Figure 1 presents further evidence of this, highlighting the fact that simulating nonnormal variables with different marginal distributions tend to generate more solutions outside [−1, +1]. The pattern is somewhat different for r = .1, where the overlap between the skewness/kurtosis combinations is particularly pronounced, so only two clouds of points are visible, one at lower kurtosis values in the parabola and other at much higher ones. The substantial overlap between the gray and black areas on the leftmost panel also points toward the fact that, contrary to the case of r = .5, if the final correlation is r = .1, then generating data with the same marginal distributions result in a larger number of cases where the solutions to the intermediate correlation equation are all strictly greater than 1.

Figure 2 shows the second special case of the intermediate correlation equation where either two or three solutions were found within the [−1, +1] range. In this instance, the pattern is somewhat reversed from what was found in Figure 1, given that a larger combination of skewness and kurtosis values resulted in multiple valid solutions for the effect size 0.1 as opposed to the effect size of 0.5. It also seems that, whereas the combination of higher order moments for r = .1 was more spread out along the Fleishman (1978) parabola, for r = .5, these combinations of moments tended to cluster themselves at both the lower and higher ends of the kurtosis dimension, where combinations of lower $(γ_{1}, γ_{2})$ values would need to be paired with larger ones in order to obtain multiple valid solutions to the intermediate correlation equation.

FIGURE 2.

Combinations of skewness and kurtosis values for final correlations of 0.1 and 0.5 that yield multiple valid solutions.

5. Illustrative Example From the Published Literature, $γ_{1} = 3, γ_{2} = 21$

Unidimensional marginal distributions with a skewness of 3 and kurtosis of 21 have been used within the social and behavioral sciences literature as a condition in Monte Carlo simulation studies testing “extreme” violations of the normality assumption (e.g., Berkovits, Hancock, & Nevitt, 2000; Curran et al., 1996; Nevitt & Hancock, 2000; Vallejo, Gras, & Garcia, 2007). In each case, multivariate nonnormal data were generated through the third-order polynomial method and correlated using the Vale–Maurelli transformation. This particular combination of values for higher order moments, however, falls within one of the special cases described in the previous sections where more than one intermediate correlation exists within the [−1, +1] range. More specifically, solving for the Fleishman coefficients and substituting them in Equation 1 define a random variable Y:

Y = - 0.6371182 - (0.6816322) Z + (0.6371182) Z^{2} + (0.148741) Z^{3},

so that Y now follows a nonnormal distribution with population skewness of 3 and kurtosis of 21. For illustration purposes, take the case of a large effect size of $r_{Y_{1} Y_{2}}$ = 0.5. Substituting these coefficients in the discriminant described in the third section yields the case where all the three roots are distinct and real. Following the use of this skewness/kurtosis combination within the published literature, the same marginal distributions were specified for Y₁ and $Y_{2}$ . Therefore, by substituting these repeating values in Equation 9, one obtains:

0.5 = r_{Z_{1} Z_{2}} ({0.6816322}^{2} + 6 (0.6816322) (0.148741) + 9 ({0.148741}^{2})) + 2 r_{Z_{1} Z_{2}}^{2} ({0.6371182}^{2}) + 6 r_{Z_{1} Z_{2}}^{3} ({0.148741}^{2}),

which, when solved for, yields the approximate roots $r_{1 Z_{1} Z_{2}} = - 5.9387499,$ $r_{2 Z_{1} Z_{2}} = - 0.8898729,$ and $r_{3 Z_{1} Z_{2}} = 0.7127582$ . The polynomial solutions $r_{2 Z_{1} Z_{2}}$ and $r_{3 Z_{1} Z_{2}}$ are valid intermediate correlations that satisfy the population correlation value initially set by the researcher.

Figure 3 shows a small demonstration, similar to the study presented by Beversdorf and Sa (2011), where the power of the t test for the Pearson correlation is compared across conditions of bivariate normality and nonnormality $(γ_{1} = 3, γ_{2} = 21)$ using the intermediate correlations of −0.8898729 and 0.7127582 for the Vale–Maurelli transformation. To minimize Monte Carlo error, the same standard normal variables were used in each power calculation, so that the only changes were the Fleishman coefficients used to induce nonnormality and the value of the intermediate correlation equation with 100,000 replications across sample sizes, N = 10, 11, 12, 13,…, 98, 99, 100. The test showed its greatest power when the assumption of normality was satisfied but, when it is not, there is a differential effect depending on which intermediate correlation is used.

FIGURE 3.

Power curves of the t test for the Pearson correlation coefficient under bivariate normality (none condition) and two intermediate correlation equations (−0.8898 and 0.7127). Population effect size (r) is set at 0.5 with sample sizes from 10 to 100.

Overall, the negative intermediate correlation influenced the simulation results so that the t test showed less power across all sample sizes when compared to its positive counterpart. The power differences across sample sizes were sometimes close to 0.1 between the negative and positive intermediate correlations, particularly around the sample size of 50 where the curve begins to form. The results also demonstrate that, even at the largest sample size of N = 100, the data generated using the negative intermediate correlation resulted in a power curve that does not reach the upper bound of 1, in contrast with the condition of bivariate normality or nonnormality using the positive intermediate correlation. No small sample bias of the correlation estimate was found in any of the nonnormal cases.

6. Conclusion

The Vale and Maurelli (1983) multivariate extension of the Fleishman (1978) method is one of the most widely used simulation algorithms within the simulations within the social and behavioral sciences. In spite of its popularity, the solution multiplicity of power polynomial transformations within the context of data-generating algorithms has received very little attention within the published literature. Although alluded to the studies by Steiger (2014), Astivia and Zumbo (2015), and Foldnes and Grønneberg (2015), the only in-depth treatment the authors have been able to find the empirical impact that these solutions have on simulation studies is in the doctoral work of Kraatz (2011), Luo (2011), and the article published by Astivia and Zumbo (2018). To the authors’ knowledge, a description or even the mere acknowledgment that the intermediate correlation equation also has multiple solutions has mostly gone unnoticed within the published literature. This note, therefore, provided a detailed account of this property and offers the potential reader some guidance concerning how to proceed when using this method.

For the cases, such as the one presented in the fourth section, where multiple valid solutions are possible, the general recommendation would be to investigate both cases and report any potential differences that may be found. Since the solution multiplicity of the Vale–Maurelli approach is a feature, and not an error, of the method itself, there are no immediate criteria to judge whether one solution is more proper or “valid” than another. For cases where all solutions to the intermediate correlation equation are outside the [−1, +1] range, it is perhaps advisable to make a note of it and use a different combination of marginal skewness and kurtosis values, since the intermediate correlation equation will no longer be a valid correlation equation and cannot be used in the initial step where multivariate normal data are generated.

Although other data generation methods are available in the literature, no one method can possibly encompass all types of nonnormal data and various options exist for applied researchers (Astivia & Zumbo, 2017). Methods such as Ruscio and Kaczetow’s (2008) iterative algorithm, Kowalchuk and Headrick’s (2010) multivariate g-and-h distributions, Mair, Satorra, and Bentler’s (2012) copula approach, or Foldnes and Olsson’s (2016) independent generator rely on known distributions with well-documented theoretical properties that circumvent the issue of needing to transform the data from normality to nonnormality. Although they are an important contribution to the literature, many of these methods do not have well-documented properties. To the best of the authors’ knowledge, for instance, the probability density function of these methods is not known. They can also be computationally intensive, such as the Ruscio and Kaczetow’s (2008) iterative algorithm or Kowalchuk and Headrick’s (2010) multivariate g-and-h distributions. Users of these methods in simulation studies need to consider some of these properties in their designs. Methods that rely on power transformations such as the Fleishman (1978) method, the Vale and Maurelli (1983) extension commented on here, or Headrick’s (2002) fifth-order approach are generally straightforward to implement, but since they rely on higher order polynomials, both for the marginal skewness and kurtosis as well as for the intermediate correlation, they can offer multiple solutions which may or may not provide sensible results.

There are several future research avenues presented in this note. A particularly important one is documenting the empirical impact that the multiple solutions to the intermediate correlation equation can have on applied Monte Carlo simulation research as shown in the fifth section. Yuan, Bentler, and Zhang (2005) have shown analytically how the fourth-order moment influences the standard error of quadratic statistics such as variances and covariances. For this particular case, the data were generated from Fleishman-transformed distributions, whose fourth moments are known, thanks to the work of Foldnes and Grønneberg (2017) in terms of intermediate correlations. By using different intermediate correlations, a differential effect of the power of the t test for the product–moment correlation became apparently contingent on which intermediate correlation was used. Extending this type of work to more complex scenarios involving multivariate statistical techniques would highlight the fact that, when using the third-order polynomial transform, researchers need to be aware of the solution multiplicity property of both of the Fleishman system shown in Equations 2 through 5 and the intermediate correlation Equation 9. As computer power becomes more accessible to researchers, simulation studies are becoming more prevalent within the published literature. It is, therefore, of utmost importance to understand the properties of these algorithms and how they can potentially influence a simulation design and the conclusions obtained from it.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Note

References

Astivia

O. L. O.

Zumbo

B. D.

(2015). A cautionary note on the use of the Vale and Maurelli method to generate multivariate, nonnormal data for simulation purposes. Educational and Psychological Measurement, 75, 541–567. doi:10.1177/0013164414548894

Astivia

O. L. O.

Zumbo

B. D

. (2017). Population models and simulation methods: The case of the Spearman rank correlation. British Journal of Mathematical and Statistical Psychology. Advance online publication. doi:10.1111/bmsp.12085

Astivia

O. L. O.

Zumbo

B. D.

(2018). On the solution multiplicity of the Fleishman method and its impact in simulation studies. British Journal of Mathematical and Statistical Psychology. doi:10.1111/bmsp.12126

Berkovits

Hancock

G. R.

Nevitt

. (2000). Bootstrap resampling approaches for repeated measure designs: Relative robustness to sphericity and normality violations. Educational and Psychological Measurement, 60, 877–892. doi:10.1177/00131640021970961

Beversdorf

(2011). Tests for correlation on bivariate non-normal data. Journal of Modern Applied Statistical Methods, 10, 699–709. doi:10.22237/jmasm/1320121680

Cohen

(1992). A power primer. Psychological Bulletin, 112, 155–159. doi:10.1037/0033-2909.112.1.155

Curran

P. J.

West

S. G.

Finch

J. F.

(1996). The robustness of test statistics to nonnormality and specification error in confirmatory factor analysis. Psychological Methods, 1, 16–29. doi:10.1037/1082-989X.1.1.16

Demirtas

Hedeker

Mermelstein

R. J.

(2012). Simulation of massive public health data by power polynomials. Statistics in Medicine, 31, 3337–3346. doi:10.1002/sim.5362

Falk

C. F

. (2017). Are robust standard errors the best approach for interval estimation with nonnormal data in structural equation modeling? Structural Equation Modeling: A Multidisciplinary Journal, 1–23. Advance online publication. doi:10.1080/10705511.2017.1367254

10.

Fleishman

A. I.

(1978). A method for simulating non-normal distributions. Psychometrika, 43, 521–532. doi:10.1007/BF02293811

11.

Flora

D. B.

Curran

P. J.

(2004). An empirical evaluation of alternative methods of estimation for confirmatory factor analysis with ordinal data. Psychological Methods, 9, 16–29. doi:10.1037/1082-989X.1.1.16

12.

Foldnes

Grønneberg

(2015). How general is the Vale–Maurelli simulation approach? Psychometrika, 80, 1066–1083. doi:10.1007/s11336-014-9414-0

13.

Foldnes

Grønneberg

(2017). The asymptotic covariance matrix and its use in simulation studies. Structural Equation Modeling: A Multidisciplinary Journal, 1–16. Advance online publication. doi:10.1080/10705511.2017.1341320

14.

Foldnes

Olsson

U. H.

(2016). A simple simulation technique for nonnormal data with prespecified skewness, kurtosis, and covariance matrix. Multivariate Behavioral Research, 51, 207–219. doi:10.1080/00273171.2015.1133274

15.

Headrick

T. C.

(2002). Fast fifth-order polynomial transforms for generating univariate and multivariate non-normal distributions. Computational Statistics and Data Analysis, 40, 685–711. doi:10.1016/S0167-9473(02)00072-5

16.

Headrick

T. C.

(2004). Transformations for simulating multivariate distributions. Journal of Modern Applied Statistical Methods, 3, 65–71. doi:10.22237/jmasm/1083370080

17.

Headrick

T. C.

(2010). Statistical simulation: Power method polynomials and other transformations. Boca Raton, FL: CRC Press.

18.

Higham

N. J.

(2002). Computing the nearest correlation matrix—A problem from finance. IMA Journal of Numerical Analysis, 22, 329–343. doi:10.1093/imanum/22.3.329

19.

Kaiser

H. F.

Dickman

(1962). Sample and population score matrices and sample correlation matrices from an arbitrary population correlation matrix. Psychometrika, 27, 179–182. doi:10.1007/BF02289635

20.

Kowalchuk

R. K.

Headrick

T. C.

(2010). Simulating multivariate g-and-h distributions. British Journal of Mathematical and Statistical Psychology, 63, 63–74. doi:10.1348/000711009X423067

21.

Kraatz

(2011). Investigating the performance of procedures for correlations under nonnormality (Unpublished doctoral dissertation). Vanderbilt University, Nashville, TN.

22.

Luo

(2011). Some aspects on confirmatory factor analysis of ordinal variables and generating non-normal data (Unpublished doctoral dissertation). Uppsala University, Uppsala, Sweden.

23.

Mair

Satorra

Bentler

P. M.

(2012). Generating nonnormal multivariate data using copulas: Applications to SEM. Multivariate Behavioral Research, 47, 547–565. doi:10.1080/00273171.2012.692629

24.

Mollin

R. A.

(1947). Algebraic number theory. Boca Raton, FL: CRC press.

25.

Nevitt

Hancock

G. R.

(2000). Improving the root mean square error of approximation for nonnormal conditions in structural equation modeling. Journal of Experimental Education, 68, 251–268. doi:10.1080/00220970009600095

26.

R Core Team. (2017). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from http://www.R-project.org/

27.

Ruscio

Kaczetow

(2008). Simulating multivariate nonnormal data using an iterative algorithm. Multivariate Behavioral Research, 43, 355–381. doi:10.1080/00273170802285693

28.

Saff

E. B.

Snider

A. D.

(2003). Fundamentals of complex analysis with applications to Engineering and Science. Upper Saddle River, NJ: Prentice Hall.

29.

Steiger

J. H.

(2014, October). Still crazy after all these years: Complexity, principles, and practice in multivariate statistics. Paper presented at the Society of Multivariate Experimental Psychology, Nashville, TN.

30.

Vale

C. D.

Maurelli

V. A.

(1983). Simulating multivariate nonnormal distributions. Psychometrika, 48, 465–471. doi:10.1007/BF02293687

31.

Vallejo

Gras

J. A.

Garcia

M. A.

(2007). Comparative robustness of recent methods for analyzing multivariate repeated measures designs. Educational and Psychological Measurement, 67, 410–432. doi:10.1177/0013164406294777

32.

Yuan

K. H.

Bentler

P. M.

Zhang

(2005). The effect of skewness and kurtosis on mean and covariance structure analysis: The univariate case and its multivariate implication. Sociological Methods & Research, 34, 240–258. doi:10.1177/0049124105280200

A Note on the Solution Multiplicity of the Vale–Maurelli Intermediate Correlation Equation

Abstract

Keywords

1. The Third-Order Polynomial Transformation

2. Limitations of the Vale–Maurelli Transformation

3. Multiple Solutions to the Vale–Maurelli Transformation: A Mathematical Description

4. Multiple Solutions to the Vale–Maurelli Transformation: A Computer Simulation Exploration

5. Illustrative Example From the Published Literature, γ 1 = 3 , γ 2 = 21

6. Conclusion

Footnotes

Declaration of Conflicting Interests

Funding

Note

References

5. Illustrative Example From the Published Literature, $γ_{1} = 3, γ_{2} = 21$