Normrank Correlations for Testing Associations and for Use in Latent Variable Models

1 Simulation 1: Which Values of δ to Use?

The blom function in rcompanion (Mangiafico, 2020) calculates equation (3) with δ values for 0 (Van der Waerden’s), π ∕ 8 ≈ 0.393 (Elving’s), 1 ∕ 3 ≈ 0.333 (Tukey’s), 3 ∕ 8 = 0.375 (Blom’s), and 1 ∕ 2 = 0.5 (the rankit method). Here, the values in δ ∈ [ 0 , 1 ) are examined. Values of δ from 0 to 0.99 in steps of 0.01, from 0.99001 to 0.99999 in steps of 0.00001, and 0.999999999, and used. Values very near one are used because δ = 1 produces ∞ for the rank n value so there is interest in the transformation’s behavior near this threshold. The code for this and the other simulations are presented in Appendix A.

Five levels of discreteness in the data are examined. The first are data drawn from a uniform distribution with a range of (0, 10), though any continuous distribution of size n will produce the same transformed values. The next level uses integers from 1 to 20, created by doubling the first variable, adding 0.5, and rounding. Analogous approaches were used to create ten point, five point, and three point integer scales, similar to those often used in surveys and personality research. This approach will produce roughly equal numbers of values in each bin. Sample sizes of n = 50 , n = 100 , and n = 500 were used. Thus, there are 5 × 3 = 15 conditions. The W from the Shapiro–Wilk (Shapiro & Wilk, 1965) test, as implemented in shapiro.test in R using Royston’s algorithm (1995), is used as a measure to compare distance from normality. W is the square of a Pearson correlation so cannot exceed one.

The values are stored, plotted in Figure 1, and the maximum and minimum W values, along with their δ values, are presented in Table 2. Note that for the case of continuous data, individual W values are not shown because each replication produces the same transformed data and therefore the same W . The results are presented as follows:

1. Little effect in the changes in sample size between n = 50 , n = 100 , and n = 500 ,

2. The more continuous a variable is, the more closely it can approximate the normal distribution, which is itself continuous,

3. There appears very little differences among the W s for the δ s ∈ ( 0 , 1 ] for all the discrete variables,

4. There is an effect of δ for the continuous variable with the W reaching a maximum around 0.7 (with a maximum W rounding to 1.000) and then dropping off when δ approaches 1. The dropping off does not happen for the discrete distributions because several values have the extreme bins, and therefore, their mid-rank is well above 1 or below n . If samples drawn from populations with fewer values that would be in the extreme bins, they would drop off too as sometimes only one value would be in the extreme bins, and therefore would have a rank of 1 or n , and a transformed value of it being infinite.

5. The δ that produces the highest W with discrete values tend to be higher than 0.7, particularly as the sample size increases.

Figure 1

Finding W values with continuous and discrete scales for different δ values, for three sample sizes. The black line is the continuous data, the red line is for the 20-point scale, the green line for the 10-point scale, the darker blue line for the 5-point scale, and the light blue line for the 3-point scale.

2 How are Data Typically Distributed?

Simulation is used in this article to explore the value of the normrank correlation measure. Simulations are done to examine testing for an association, exploratory factor analysis (EFA), and SEM. A selection of the possible data creation scenarios are reported. The aim is to show a selection of uses of this normrank procedure to illustrate its potential, rather than to provide an exhaustive set of possible data creation scenarios.

Before discussing the main simulations reported in this article, it is necessary to examine what types of distributions are typical in education and allied disciplines. The distribution of variables used in correlations varies by discipline and by many other contextual details. The interest here is with kurtosis (for discussion of this statistic, see DeCarlo, 1997; Wright & Herrington, 2011) because of its effect on power (Wilcox, 2017). There are several types of kurtosis based on the fourth moment and more robust alternatives (e.g., Moors, 1988). These are reviewed by Joanes and Gill (1998). Type 3 kurtosis is the default of the function kurtosis found in the e1071 package (Meyer, Dimitriadou, Hornik, Weingessel, & Leisch, 2018) and is reported here.

Micceri (1989) requested data from several education sources, test publishers, and authors of psychology and education journal articles. He analyzed 440 separate distributions with respect to kurtosis. He classified their kurtosis as follows: similar to the uniform distribution (3.2% of his sample), less than the normal distribution (14.8%), similar to the normal distribution (15.8%), similar to a moderate contaminated normal distribution (17.7%), similar to an extreme contaminated normal distribution (32.5%), and similar to a double exponential distribution (16.6%). He defined a moderate contaminated distribution as 95% drawn from a normal distribution with σ = 1 and 5% drawn from a normal distribution with σ = 2 . Extreme contaminated normal was defined as 85% drawn from a normal distribution with σ = 1 and 15% drawn from a normal distribution with σ = 3 . Thus, about half of the distributions he examined had kurtosis at least that of what he called the extreme contaminated normal distribution. He found approximately 18% to have kurtosis less than normal.

Paraphrasing Poincaré, everybody believes the normal distribution assumption is justified: the mathematicians – because they think the scientists say all their data are normally distributed; the scientists – because they think the mathematicians say it is justified because of the central limit theorem. Micceri (1989) titled his paper “The unicorn, the normal curve, and other improbably creatures” to stress that real data are seldom normally distributed. Thus, scientists do not say that all their data are usually normally distributed. Tukey (1960) showed that nonnormally distributed data do affect statistical results, countering the idea that the central limit theory makes this assumption justified in all circumstances. Tukey’s advice was to assume that “all assumptions are wrong” (Tukey, 1986, p. 72), but then to explore how problematic it is to make them and what alternatives exist. One of the main goals of this article is to examine the costs of assuming data are normally distributed and discuss an alternative.

Four distributions are considered here to illustrate how the three correlations perform: the uniform distribution, the normal distribution, a contaminated normal distribution with 15% with σ = 3 , and the χ 2 with df = 1 , which is equivalent to the square of a normal distribution. The normal and contaminated normal distributions are shown in Figure 2. They can appear very similar in histograms, but the tails of the contaminated normal distribution are larger (right panel). Tukey showed that this can have a large impact on the power of the tests as the standard deviations can become much larger. The kurtosis for these four distributions are: uniform equals 1.80, normal equals 0, χ 1 2 equals 12, and the contaminated normal used equals approximately 5.06. The χ 1 2 distribution has skewness of 8 ≈ 2.83 . The others are symmetric. An example of how to perform a normrank correlation is in the first section of Appendix B.

Figure 2

The normal distribution ( μ = 0 , σ = 1 ) and a contaminated normal distribution that includes 15% with σ = 3 . The right panel focuses on the tails with x > 2 . This makes the heavier tail of the contaminated distribution clearer.

3 Simulation 2: Testing for an Association

Pearson’s correlation, Spearman’s correlation, and the normrank correlation are compared for whether they maintain the stated Type 1 error rate, here 5%, and for their power to detect an association. Other tests of association exist (e.g., polychoric correlations), but the desire was to compare these because they all use the same procedure after transforming the original data (with no transformation, ranking, or with equation (3) using δ = 0.7 ) and the popularity of the first two procedures. A concern for Pearson’s approach is that it is greatly influenced by extreme values (e.g., Wilcox, 2017). The ranking transformation used with Spearman’s correlation directly addresses this so that all values are equally spaced, provided that there are no ties (ties, and discrete variables generally, are not considered in this article). Compared with ranking, equation (3) changes this so that more extreme rank values are spaced more than less extreme rank values. Whether these extreme values are spaced further apart than the original data depends on the distribution of the original data (Table 1).

First, the performance of these measures are considered, when there is no association. Sample sizes of 40, 60, and 100 were used. The distributions-uniform, normal, contaminated normal, and χ 1 2 – were the same for both variables. The χ 2 distribution with d f = 1 , which is equivalent to squaring a normally distributed variable with μ = 0 and σ = 1 , was included so that one of the variables is skewed (skewness is 8 ≈ 2.83 ). There were 10,000 replications per condition for the remaining simulations reported in this article. The R code for these is in Appendix A and can be adapted for different sample sizes, distributions, and other parameters. The p -values are found for the Pearson, Spearman, and normrank correlations using R’s cor.test function. The nominal α is 0.05, so the proportion of p < 0.05 should be near this value. The proportions of p < 0.05 results and the 95% confidence intervals for these proportions, found with the Wilson method, are reported. Because the number of replications is large, the widths of these intervals are small. Differences of 0.01 between the correlation procedures are large enough that they are discussed, and differences greater than 0.02 may be considered substantial in some contexts.

Table 3, shows for each procedure, for each distribution, and for each sample size, the proportion of p -values less than 0.05 and the confidence intervals of these proportions. All are near the nominal level of 0.05. None of the individual cell proportions is greater than 0.01 away from the nominal level. Therefore, it is concluded that all are adequate at controlling Type 1 error.

This simulation was repeated, but with the null hypothesis being false. The variables were created using the three distributions but for one of the variables, 0.32 times the other was added. This results in a Pearson correlation of approximately r = 0.30 , what Cohen (1992) called a medium sized effect, for the normally distributed variables. The correlations are slightly different for the other distributions and other measures. The primary interest here is whether the different statistical procedures vary in how well they detect an association for the same data, not differences across the distributions.

Table 4 shows the proportion of the replications where the null hypothesis was correctly rejected for each procedure, for each distribution, and for each sample size. As expected, these proportions increase with sample size (i.e., the power increases with the sample size). For the normal distribution, Pearson’s correlation correctly rejects H o about 1% more often than the normrank correlation and about 4–5% more than Spearman’s correlation. Thus, in this situation, Spearman’s is the least powerful and Pearson’s is the most powerful. With both the uniform and contaminated normal distributions, the normrank procedure had about a 3–4% advantage over Pearson’s with Spearman’s performing between the two. Based on Micceri (1989) findings, where most distributions he analyzed had greater than normal kurtosis and a small proportion less than normal, it is likely that education and other social science researchers will face circumstances where the normrank is more powerful than Pearson’s procedure. In none of the situations examined was Spearman’s procedure the most powerful procedure. It is recommended not to use Spearman’s correlation, and it will not be further examined in this article.

4 Simulation 3: EFA

In this simulation, Pearson’s correlation and the normrank correlation are compared for EFA when there are no underlying factors and when there are 1–4 underlying factors with 12 observed variables. Spearman’s is not included as it performed worse than others for all distributions tested in the previous simulation.

EFA is one of the most used forms of the latent variable model particularly within the educational and social sciences. It is usually the first latent variable model that students are taught. If you have some number of observed variables and believe that this observed set is caused by a smaller number of latent variables, but you do not have particular relationships hypothesized, EFA can be an appropriate tool. EFA produces two main results. First, it can suggest how many latent variables are appropriate for modeling the data. This is often done using a related technique: calculating the eigenvalues of the correlation matrix. The second result is that it shows which observed variables are related to the same latent variable as others. There are several sources for EFA, for example, in increasing level of assumed mathematical knowledge: Bartholomew, Knott, and Moustaki (2011), Loehlin and Beaujean (2017), Mair (2018), Mulaik (2010), and Wright and Wells (2020).

The simulation in this section will examine whether using a Pearson correlation matrix or a matrix composed using the normrank transformation is better able both to determine the correct number of latent variables and to identify the loadings correctly. Velicer, Eaton, and Fava (2000) (see also Auerswald & Moshagen, 2019) review different methods to determine the number of factors. Based on this, the empirical Kaiser criterion is used (Braeken & van Assen, 2017). It involves comparing the observed eigenvalues with those observed from a random data matrix. It is implemented in the EFA.dimensions package (O’Connor, 2021). To estimate whether the pattern of factor loadings is accurate, factor analysis was conducted assuming the correct number of factors. The highest loading for each observed variable was recorded. For the complete pattern of loadings to be deemed accurate, all observed variables which were influenced by the same latent variable during data creation had to have their highest loading for the same latent variable. This was done with the default settings of the R function factanal.

Only a single sample size is used here. There is no agreed upon minimum sample size for EFA (de Winter, Dodou, & Wieringa, 2009; MacCallum, Widaman, Preacher, & Hong, 2001), but several rules of thumb exist. Loehlin and Beaujean (2017) say “with high commonalities, a small number of factors, and a relatively large number of indicators per factor, N s of 100 or or less yield good recovery of factors” (p. 207). A sample size of 100 will be used with 12 observed variables since it is near the minimum that might be considered adequate (though larger samples are recommended). Each observed variable will be created as:

where x is a scalar to control the amount of noise (see below).

The latent variables and the noise variables are each drawn from one of several contaminated normal distributions. Because Micceri (1989) found the preponderance of the variables he examined to have greater than normal kurtosis, this simulation will focus on distributions with varying amounts of contamination for the contaminated normal distribution. In addition, results for the χ 1 2 distribution are presented for comparison to show a skewed distribution.

Contamination for contaminated normal distributions can be created in different ways. Here, the variability of the two distributions are fixed at σ 1 = 1 and σ 2 = 3 and the proportion from each of these distributions is varied between 0 and 1, with the endpoints being normal distributions. This is done in 101 steps of 0.01. The kurtosis for these different contaminated normal distributions is shown in Figure 3. The maximum kurtosis is when about 10% is drawn from a normal distribution with σ = 3 . The uniform distribution is presented for comparison.

Figure 3

Kurtosis for normal distributions with σ = 1 for proportions from 0 to 1 with the remainder being drawn from normal distributions with σ = 3 , as well as the kurtosis for the uniform distribution.

Different amounts of noise are used for the different numbers of latent variables to prevent results being at ceiling (correctly identifying that there are four latent variables is more difficult than identifying that there is one). The scalars for noise are x = 2.5 , x = 2 , x = 1.5 , and x = 1.2 for one, two, three, and four latent variables, respectively.

Figures 4 and 5 show that when the data are all produced by noncontaminated normal distributions (mixtures of 0 and 1), Pearson’s and the normrank correlation perform similarly for both correctly identifying the number of factors and identifying the correct pattern of loadings. They perform similarly until the distribution has less than about 50% drawn from the one with the larger standard deviation. This is also where the kurtosis becomes noticeable larger (Figure 3). For these distributions, the normrank procedure was more accurate both for identifying the correct number of factors and for identifying the correct pattern of factor loadings.

Figure 4

Proportion of replications identifying the correct number of latent factors. Note that y -axes are different for the five plots. The proportion from the population with σ = 3 increased from 0 to 1 in steps of 0.01, with 10,000 replications at each step.

Figure 5

Proportion of time identifying the correct loadings (i.e., the highest loading being the same as the others influenced by the same factor). The proportion from the population with σ = 3 increased from 0 to 1 in steps of 0.01, with 10,000 replications at each step. Note that the y -axis scales are different.

Table 5 shows the results when the data are distributed χ 1 2 . When there is more than one factor used to create the data, the normrank correlation is more accurate than Pearson’s correlation for both identifying the number of factors and the factor loading pattern.

There are many different parameters that could be used to construct data for EFA and different ways to assess the accuracy of the results. The results presented here are just to illustrate the use of the two correlation procedures. For these, and consistent with the correlation simulation, when the data derive from normal distributions the Pearson’s and the normrank correlation perform similarly, but for nonnormal distributions, both those with less than normal kurtosis (the uniform distributions) and those with greater than normal kurtosis (the contaminated normal distributions), the normrank performed better. Given Micceri’s (1989) finding that most of the distributions he analyzed in psychology and education were not similar to the normal distribution, and most of these had greater than normal kurtosis, the recommendation here is that the normrank correlation should be preferred over Pearson’s correlation when conducting EFA for many psychology and education data sets.

5 Simulation 4: SEM

SEM can take numerous forms. It allows the researcher to hypothesize relationships among latent variables. However, it is a complex analysis that should be used cautiously.

When we come to models for relationships between latent variables we have reached a point where so much has to be assumed that one might justly conclude that the limits of scientific usefulness have been reached, if not exceeded.

Bartholomew et al. (2011, p. 228)

A good applied introduction on SEM using R is Beaujean (2014), with more coverage in Loehlin and Beaujean (2017). The lavaan manual Rosseel (2012), one of the main R SEM packages, is a good resource. Here, only a single example is shown to show how the normrank correlation compares with Pearson’s correlation using three distributions (uniform, normal, contaminated normal) for a relatively simple SEM model. A single situation is used to explore if the findings are consistent with the earlier simulations. There are numerous variations in SEM models that could be explored further. As with the other simulations, within each data creation scenario, all the random variables will have the same distributions: uniform (range 0 to 1), normal ( μ = 0 , σ = 1 ), and contaminated normal (85% drawn from a normal with μ = 0 and σ = 1 and 15% drawn from a normal with μ = 0 and σ = 3 ).

The data will be created from the model depicted in Figure 6. There are three latent variables, l v 1 − l v 3 , each influencing three observed variables, labeled o 1 − o 9 . l v 1 influences l v 2 , and l v 2 influences l v 3 . For half of the simulations, l v 1 also influences l v 3 . The statistical models compare the data creation models with and without this effect and evaluated for whether this effect is detected at p < 0.05 . Given the nature of SEM modeling, sometimes models fail to converge. Often alternative estimation procedures and different control settings can be used to achieve convergence. Here, the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm is used, but some errors persisted and are treated as missing. Only 0.12% of the 20,000 replications failed to produce a result (Table 6).

Figure 6

Causal model used to create the data for SEM simulation. The effect corresponding with the arrow surrounded by the dashed ellipse is included for half the samples.

The primary SEM results can be summarized as follows: Pearson’s and the normrank produce similar findings when the distributions are uniform or normal, but the normrank approach performs better in terms for both not rejecting the true null (by a couple of percent) and correctly rejecting a false null hypothesis (by about seven percent) when the distributions have high kurtosis. Large differences were found for the different distributions for incorrectly rejecting a true null (i.e., the Type 1 error proportions). In the uniform conditions, the true null was incorrectly rejected less than the nominal level of 0.05, while in the normal conditions, it was rejected slightly more than the nominal level, and in the contaminated normal conditions, it was rejected much more. These differences highlight the importance of considering how the distributions underlying the data (real and simulate) can affect statistical models.

6 Summary

Pearson’s product moment correlation is one of the most used statistics in all of science. It is highly influenced by extreme values, which makes this procedure less reliable for some distributions. Spearman’s rank correlation is a popular alternative when there are extreme outliers and the analyst wants to limit their impact. The normrank correlation, using what is sometimes called the rank based inverse normal transformation, results in data more closely resembling a normal distribution.

Four simulations were conducted. The first was to show which variation of the normrank transformation to use. For continuous distributions, δ = 0.7 appeared to perform well. This was used in the remaining simulations. The second simulation examined how well Pearson’s, Spearman’s, and the normrank correlation performed detecting an association. Spearman’s was never the best of these three. The recommendation is not to use it. Pearson’s was more powerful than the normrank correlation with the normally distributed data, but the difference was fairly small. The normrank correlation was more powerful than Pearson’s correlation when the distribution had sub- and super-normal kurtosis. Because Micceri (1989) found most correlations in psychology and education had greater than normal kurtosis, the normrank correlation is likely to be more powerful and should be considered. However, it is important to stress that Pearson’s approach was better for normally distributed data, as has been shown in other contexts (e.g., Beasley et al., 2009), using the normrank transformation will not always produce more powerful statistical tests. A further reservation for all procedures involving ranking is that information about the original variable metric is lost.

The final two simulations compared the use of Pearson’s correlation with the normrank correlation in factor analysis and SEM. The results were similar to the second simulation. When the data are drawn for normal distributions, the procedures perform similarly, but when they are drawn from nonnormal distributions, with kurtosis much different from the normal distribution, the normrank performed better. Overall,

It is concluded, therefore that:

1. Spearman’s procedure should not be used.

2. The normrank should be considered unless the researcher has strong reasons to believe their data are drawn for near-normal distributions, but researchers may also wish to consider numeric transformations.

It is important to stress that this article only addresses one alternative to Pearson’s correlation on untransformed variables. The article shows that the normrank transformation can improve the subsequent statistics, but there are other approaches. First, if there is a particular substantive model that is consistent with a particular transformation (e.g., count data, which according to Mosteller & Tukey, 1977, may profit from the ln transformation), this should be considered. Second, the Box-Cox transformation (Box & Cox, 1964) usually makes data more similar to the normal distribution and is often used. If you have a particular model, whether linear or otherwise, about the relationship between the variables no transformation involving ranking should be used and the applicable model should be evaluated. If the worry is about the impact of extreme values more robust loss functions could be used (e.g., Hoaglin et al., 2000) and if the worry is about not meeting the distributional assumptions bootstrapping could be used (e.g., Efron & Tibshirani, 1993; Wright, London, & Field, 2011). In short, there are many alternatives to conducting Pearson’s correlation on untransformed variables, and in some situations, the normrank transformation data should be considered.