Formulas for Pearson correlations with and without weights

The weighted Pearson correlation is computed using the formula $$r_{Pearson}=\frac{\sum_{i=1}^n \left[ w_i (x_i-\bar{x})(y_i-\bar{y}) \right]}{\sqrt{\sum_{i=1}^n \left( w_i (x_i-\bar{x})^2 \right)\sum_{i=1}^n \left( w_i (y_i-\bar{y})^2 \right) }} $$

where w_i is the weights, x̄ is the weighted mean of the X variable ($\bar{x}=\frac{1}{\sum_{i=1}^n w_i}\sum_{i=1}^n w_i x_i$), ȳ is the weighted mean of the Y variable ($\bar{y}=\frac{1}{\sum_{i=1}^n w_i}\sum_{i=1}^n w_i y_i$), and n is the number of elements in X and Y.³

The unweighted Pearson correlation is calculated by setting all of the weights to one.

Formulas for Spearman correlations with and without weights

For the Spearman correlation coefficient the unweighted coefficient is calculated by ranking the data and then using those ranks to calculate the Pearson correlation coefficient–so the ranks stand in for the X and Y data. Again, similar to the Pearson, for the unweighted case the weights are all set to one.

For the unweighted case the highest rank receives a value of 1 and the second highest 2, and so on down to the nth value. In addition, when data are ranked, ties must be handled in some way. The chosen method is to use the average of all tied ranks. For example, if the second and third rank units are tied then both units would receive a rank of 2.5 (the average of 2 and 3).

For the weighted case there is no commonly accepted weighted Spearman correlation coefficient. Stata does not estimate a weighted Spearman and SAS does not document their methodology in either of the corr or freq procedures.

The weighted case presents two issues. First, the ranks must be calculated. Second, the correlation coefficient must be calculated.

Calculating the weighted rank for an individual level is done via two terms. For the jth element the rank is

rank_j = a_j + b_j

The first term a_j is the sum of all weights W less than or equal to this value of the outcome being ranked (ξ_j)

$$a_j = \sum_{i=1}^n w_i \mathbf{1}\left( \xi_i < \xi_j \right)$$

where 1(⋅) is the indicator function that is one when the condition is true and 0 when the condition is false, w_i is the ith weight and ξ_i and ξ_j are the ith and jth value of the vector being ranked, respectively.

The term b_j then deals with ties. When there are ties each unit receives the mean rank for all of the tied units. When the weights are all one and there are n tied units the vector of tied ranks would be v = (a_j + 1, a_j + 2, …, a_j + n). The mean of this vector (here called rank¹ to indicate it is a specific case of rank when the weights are all one) is then

$$rank_j^1=\frac{1}{n} \sum_{i=1}^n \left(a_j + i \right)$$ $$=\frac{1}{n} \left( n a_j + \frac{n(n+1)}{2} \right)$$ $$=a_j + \frac{n+1}{2}$$

thus

$$b_j^1=\frac{n+1}{2}$$

where the superscript one is again used to indicate that this is only for the unweighted case where all weights are set to one.

For the weighted case this could be $\mathbf{v}=\left(a_j+w_1', a_j+w_1'+w_2', \dots, a_j+\sum_{k=1}^n w_k' \right)^T$ where W’ is a vector containing the weights of the tied units. It is readily apparent that the mean of this vector value will depend on the ordering of the weights. To avoid this, the overall mean of all possible permutations of the weights is calculated. The following formula does just that

$$b_j = \frac{n+1}{2}\bar{w}_j$$

where w̄_j is the mean weight of all of the tied units. It is easy to see that when the weights are all one w̄_j = 1 and b_j = b_j¹. The latter (more general) formula is used for all cases.

After the X and Y vectors are ranked they are plugged into the weighted Pearson correlation coefficient formula shown earlier.

##Formulas for polyserial correlation with and without weights For the polyserial correlation, it is again assumed that there are two continuous variables X and Y that have a bivariate normal distribution.⁴

$$\begin{pmatrix} X \\ Y \end{pmatrix} \sim N \left[ \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix}, \boldsymbol{\Sigma} \right]$$

where N(A, Σ) is a bivariate normal distribution with mean vector A and covariance matrix Σ. For the polyserial correlation, Y is discretized into the random variable M according to

$$m_i= \begin{cases} 1 \quad \mathrm{if} \theta_2 < y_i < \theta_3 \\ 2 \quad \mathrm{if} \theta_3 < y_i < \theta_4 \\ \vdots \\ t \quad \mathrm{if} \theta_{t+1} < y_i < \theta_{t+2} \end{cases}$$

where θ indicates the cut points used to discretize Y into M, and t is the number of bins. For notational convenience, θ₂ ≡ −∞ and θ_t + 2 ≡ ∞.⁵

To give a concrete example, the following figure shows the density of Y when the cuts points are, for this example, θ = (−∞, −2, −0.5, 1.6, ∞). In this example, any value of −2 < y_i < −0.5 would have m_i = 2.

Figure 1 :. Density of Y for cutpoints θ = (−∞, −2, −0.5, 1.6, ∞).

hi

Figure 1 :. Density of Y for cutpoints θ = (−∞, −2, −0.5, 1.6, ∞).

Notice that μ_y is not identified (or is irrelevant) because, for any a ∈ ℝ, setting μ̃_y = μ_y + a and $\tilde{\boldsymbol{\theta}}=\boldsymbol{\theta} + a$ lead to exactly the same values of M and so one of the two must be arbitrarily assigned. A convenient decision is to decide μ_y ≡ 0. A similar argument holds for σ_y so that σ_y ≡ 1.

For X, Cox (1974) observes that the MLE mean and standard deviation of X are simply the average and (population) standard deviation of the data and do not depend on the other parameters.⁶ This can be taken advantage of by defining z to be the standardized score of x so that $z \equiv \frac{x- \bar{x}}{ \hat\sigma_x}$.

Combining these simplifications, the probability of any given x_i, m_i pair is

Pr(ρ = r, Θ = θ; Z = z_i, M = m_i) = ϕ(z_i)∫_θmi + 1^{θ_mi + 2}            f(y|Z = z_i, ρ = r)dy

where Pr(ρ = r, Θ = θ; Z = z_i, M = m_i) is the probability of the event ρ = r and the cuts points are θ, given the ith data point z_i and m_i; ϕ(⋅) is the standard normal; and f(Y|Z, ρ) is the distribution of Y conditional on Z and ρ. Because Y and Z are jointly normally distributed (by assumption) $$f(Y|Z=z_i,\rho=r) =N\left(\mu_y + \frac{\sigma_y}{\sigma_z}r(z_i-\mu_z), (1-r^2){\sigma_y}^2 \right)$$ because both Z and Y are standard normals f(y|Z = z_i, ρ = r) = N(r ⋅ z_i, (1 − r²)) Noticing that $\frac{y-r\cdot z}{\sqrt{1-r^2}}$ has a standard normal distribution $$\mathrm{Pr}\left( \rho=r, \boldsymbol{\Theta}=\boldsymbol{\theta} ; Z=z_i, M=m_i \right) = \phi(z_i) \left[ \Phi\left( \frac{\theta_{m_i+2} - r \cdot z_i}{\sqrt{1-r^2}} \right) - \Phi \left( \frac{\theta_{m_i+1} - r \cdot z_i}{\sqrt{1-r^2}} \right) \right]$$

where Φ(⋅) is the standard normal cumulative density function. Using the above probability function as an objective, the log-likelihood is then maximized. $$\ell(\rho=r, \boldsymbol{\Theta}=\boldsymbol{\theta};\mathbf{Z}=\mathbf{z},\mathbf{M}=\mathbf{m}) = \sum_{i=1}^n w_i \ln\left[ \mathrm{Pr}\left( \rho=r, \boldsymbol{\Theta}=\boldsymbol{\theta} ; Z=z_i, M=m_i \right) \right]$$ where w_i is the weight of the i^th members of the vectors Z and Y. For the unweighted case, all of the weights are set to one.

The value of the nuisance parameter θ is chosen to be

θ̂_j + 2 = Φ⁻¹(n/N)

where n is the number of values to the left of the jth cut point (θ_j + 2 value) and N is the number of data points overall. Here two is added to j to make the indexing of θ agree with Cox (1974) as noted before. For the weighted cause n is replaced by the sum of the weights to the left of the jth cut point and N is replaced by the total weight of all units

$$\hat{\theta}_{j+2} = \Phi^{-1}\left( \frac{\sum_{i=1}^N w_i \mathbf{1}(m_i < j) }{\sum_{i=1}^N w_i} \right)$$

where 1 is the indicator function that is 1 when the condition is true and 0 otherwise.

###Computation of polyserial correlation For the polyserial, derivatives of ℓ can be written down but are not readily computed. When the ML argument is set to FALSE (the default), a one dimensional optimization of ρ is calculated using the optimize function in the stats package and the values of θ from the previous paragraph. When the ML argument is set to TRUE, a multi-dimensional optimization is done for ρ and θ using the bobyqa function in the minqa package. See the wCorr Arguments vignette for a comparison of these two methods.

Because the numerical optimization is not perfect when the correlation is in a boundary condition (ρ ∈ {−1, 1}), a check for perfect correlation is performed before the above optimization by simply examining if the values of X and M have agreeing order (or opposite but agreeing order) and then the MLE correlation of 1 (or -1) is returned.

##Methodology for polychoric correlation with and without weights

Similar to the polyserial correlation, the polychoric correlation is a simple case of two continuous variables X and Y that have a bivariate normal distribution. In the case of the polyserial correlation the continuous (latent) variable Y was observed as a discretized variable M. For the polychoric correlation, this is again true but now the continuous (latent) variable X is observed as a discrete variable P according to

$$p_i= \begin{cases} 1 \quad \mathrm{if} \theta'_2 < x_i < \theta'_3 \\ 2 \quad \mathrm{if} \theta'_3 < x_i < \theta'_4 \\ \vdots \\ t \quad \mathrm{if} \theta'_{t'+1} < x_i < \theta'_{t'+2} \end{cases}$$

where θ remains the cut points for the distibution defining the transformation of Y to M, θ′ is the cut points for the transformation from X to P, and t′ is the number of bins for P. Similar to θ, θ′ has θ′₂ ≡ −∞ and θ′_t′ + 2 ≡ ∞.

As in the polyserial correlation, μ_y is not identified (or is irrelevant) because, for any a ∈ ℝ, setting μ̃_y = μ_y + a and $\tilde{\boldsymbol{\theta}}=\boldsymbol{\theta} + a$ lead to exactly the same values of M and so one of the two must be arbitrarily assigned. The same is true for μ_x. A convenient decision is to decide μ_y = μ_x ≡ 0. A similar argument holds for σ_y and σ_x so that σ_y = σ_x ≡ 1

Then the probability of any given m_i, p_i pair is

Pr(ρ = r, Θ = θ, Θ′ = θ′; P = p_i, M = m_i) = ∫_{θ′pi + 1}^{θ′_pi + 2}∫_θmi + 1^{θ_mi + 2}            f(x, y|ρ = r)dydx

where ρ is the correlation coefficient.

Using this function as an objective, the log-likelihood is then maximized. $$\ell(\rho=r, \boldsymbol{\Theta}=\boldsymbol{\theta}, \boldsymbol{\Theta}'=\boldsymbol{\theta}';\mathbf{P}=\mathbf{p},\mathbf{M}=\mathbf{m}) = \sum_{i=1}^n w_i \ln\left[\mathrm{Pr}\left( \rho=r, \boldsymbol{\Theta}=\boldsymbol{\theta}, \boldsymbol{\Theta}'=\boldsymbol{\theta}' ; P=p_i, M=m_i \right) \right] $$

This is the weighted log-likelihood function. For the unweighted case all of the weights are set to one.

###Computation of polychoric correlation This again mirrors the treatment of the polyserial. The derivatives of ℓ for the polychoric can be written down but are not readily computed. When the ML argument is set to FALSE (the default), a one dimensional optimization of ρ is calculated using the optimize function from the stats package and values of θ and θ′ are computed using the last equation in the section titled, “Formulas for polyserial correlation with and without weights”. When the ML argument is set to TRUE a multi-dimensional optimization is done for ρ, θ, and θ′ using the bobyqa function in the minqa package. See the wCorr Arguments vignette for a comparison of these two methods.

Because the optimization is not perfect when the correlation is in a boundary condition (ρ ∈ {−1, 1}), a check for perfect correlation is performed before the above optimization by simply examining if the values of P and M have a Goodman-Kruskal correlation coefficient of -1 or 1. When this is the case, the MLE of -1 or 1, respectively, is returned.

#Simulation evidence on the correctness of the estimating methods

It is easy to prove the consistency of the θ for the polyserial correlation and θ and θ′ for the polychoric correlation using the non-ML case. Similarly, for ρ, because it is an MLE that can be obtained by taking a derivative and setting it equal to zero, the results are asymptotically unbiased and obtain the Cramer-Rao lower bound.

This does not speak to the small sample properties of these correlation coefficients. Previous work has described their properties by simulation; and that tradition is continued below.⁷

Simulation study of unweighted correlations

In what follows, when the exact method of selecting a parameter (such as n) is not noted in the above descriptions it is described as part of each simulation.

Across a number of iterations (the exact number of times will be stated for each simulation), the following procedure is used:

• select a true Pearson correlation coefficient ρ;
• select the number of observations n;
• generate X and Y to be bivariate normally distributed using a pseudo-Random Number Generator (RNG);
• using a pseudo-RNG, select the number of bins for M and P (t and t′) independantly from the set {2, 3, 4, 5};⁸
• select the bin boundaries for M and P (θ and θ′) by sorting the results of (t − 1) and (t′ − 1) draws, respectively, from a normal distribution using a pseudo-RNG;
• confirm that at least 2 levels of each of M and P are occupied (if not, return to previous step); and
• calculate and record relevant statistics.

Bias, and RMSE of the unweighted correlations

This sections shows the bias of the correlations as a function of the true correlation coefficient, ρ. To that end, a simulation was done at each level of the cartesian product of ρ ∈ (−0.99, −0.95, −0.90, −0.85, ..., 0.95, 0.99), and n ∈ {10, 100, 1000}. For precision, each level of ρ and n was run fifty times. The bias is the mean difference between the true correlation coefficient (ρ_i) and estimate correlation coefficient (r_i). The RMSE is the square root of the mean squared error.

$$\mathrm{RMSE}= \sqrt{ \frac{1}{n} \sum_{i=1}^n \left( r_i - \rho_i \right)^2 }$$

And the bias is given by

$$bias=\frac{1}{n}\sum_{i=1}^n \left(r_i - \rho_i \right)$$

Figure 2 : shows the bias as a function of the true correlation ρ. Only the polyserial shows no bias at any level of n, shown by no clear deviation from 0 at any level of ρ. For the Pearson correlation there is bias when n = 10 that is not present when n = 100 or 1,000. This is a well known property of the estimator.⁹ Similarly, the polychoric shows bias when n = 10.

The Spearman correlation shows bias at all of the tested levels of n. The bias is zero when the true correlation is 1, 0, or -1; is positive when ρ is below 0 (negative correlation); and is negative when ρ is above 0 (positive correlation). In this section, the Spearman correlation coefficient is compared with the true Pearson correlation coefficient. When this is done, the bias is expected because the Spearman correlation is not intended to recover a Pearson type correlation coefficient; it is designed to measure a separate quantity.

Figure 2 :. Bias Versus ρ for Unweighted Correlations.

Figure 3 : shows the RMSE as a function of ρ. All of the correlation coefficients have a uniform RMSE as a function of ρ near ρ = 0 that decreases near |ρ| = 1. All plots also show a decrease in RMSE as n increases. This plot shows that there is no appreciable RMSE differences as a functions of ρ. In addition, it show that our attention to the MLE correlation of -1 or 1 at edge cases did not make the RMSE much worse in the neighborhood of the edges (|ρ| ∼ 1).

Figure 3 :. Root Mean Square Error Versus ρ for Unweighted Correlations.

Simulation study of weighted correlations

When complex sampling (other than simple random sampling with replacement) is used, unweighted correlations may or may not be consistent. In this section the consistency of the weighted coefficients is examined.

When generating simulated data, decisions about the generating functions have to be made. These decisions affect how the results are interpreted. For the weighted case, if these decisions lead to something about the higher weight cases being different from the lower weight cases then the test will be more informative about the role of weights. Thus, while it is not reasonable to always assume that there is a difference between the high and low weight cases, the assumption (used in the simulations below) that there is an association between weights and the correlation coefficients serves as a more robust test of the methods in this package.

##Results of weighted correlation simulations

Simulations are carried out in the same fashion as previously described but include a few extra steps to accommodate weights. The following changes were made:

• Weights are assigned according to w_i = (x − y)² + 1, and the probability of inclusion in the sample was then $Pr_i = \frac{1}{w_i}$.
• For each unit, a uniformly distributed random number was drawn. When that value was less than the probability of inclusion (Pr_i), the unit was included.

Units were generated until n units were in the sample.

Two simulations were run. The first shows the mean absolute deviation (MAD)

$$MAD=\frac{1}{n}\sum_{i=1}^n |r_i - \rho_i|$$

as a function of ρ and was run for n = 100 and ρ ∈ (−0.99, −0.95, −0.90, −0.85, ..., 0.95, 0.99), with 100 iterations run for each value of ρ.

The following plot shows the MAD for the weighted and unweighted results as a function of ρ when n = 100. This shows that for values of ρ near zero, under our simulation assumptions (for all but the Spearman correlation) the weighted correlation performs better than (that is, has lower MAD than) the unweighted correlation for all correlation coefficients. Over the entire range, the difference between the two is never such that the unweighted has a lower MAD. Thus, under the simulated conditions at least, the weighted correlation has lower or approximately equal MAD for every value of the true correlation coefficient (ρ).

Figure 6 :. Mean Absolute Deviation Versus ρ (Weighted).

The second simulation (shown in Figure 7 :) used the same values of ρ and used n ∈ {10, 100, 1000, 10000} and shows how RMSE and sample size are related. In particular, it shows first-order convergence of the weighted Pearson, polyserial, and polychoric correlation coefficient.

For the previous plots the calculated Spearman correlation coefficient was compared to the generating Pearson correlation coefficient. For this plot only, the Spearman correlation coefficient to the true Spearman correlation coefficient. This is because the Spearman coefficient is not attempting to estimate the Pearson correlation. To do this the simulation is modified slightly. A population of data is generated and the true Spearman correlation coefficient then is calculated as the population coefficient.¹⁰ Then, a sample from the population with varying probability as described in the weighted simulation section is used to calculate sample Spearman correlation coefficient. Then the root mean squared difference between the sample and population coefficients are calculated as with the Pearson–except that the population Spearman correlation coefficient is used in place of the Pearson correlation coefficient (ρ).

Thus, the results in Figure 7 : show that, when compared to the true Spearman correlation coefficient, the weighted Spearman correlation coefficient is consistent.

In all cases the RMSE is lower for the weighted than the unweighted. Again, the fact that the simulations show that the unweighted correlation coefficient is not consistent does not imply that it will always be that way–only that this is possible for these coefficients to not be consistent.

Figure 7 :. Root Mean Square Error Versus ρ (Polyserial, Pearson, Polychoric panels) or Population Spearman correlation coefficient (Spearman panel) for Weighted Correlations

#Conclusion Overall the simulations show first order convergence for each unweighted correlation coefficient with an approximately linear computation cost. Further, under our simulation assumptions, the weighted correlation performs better than (has lower MAD or RMSE than) the unweighted correlation for all correlation coefficients.

We show the first order convergence of the weighted Pearson, polyserial, and polychoric correlation coefficient. The Spearman is shown to not consistently estimate the population Pearson correlation coefficient but is shown to consistently estimate the population Spearman correlation coefficient–under the assumptions of our simulation.