Choosing the Right Correlation Method: Theory and Rationale

Introduction

The Pearson product-moment correlation coefficient is the most widely used measure of linear association. However, it is only optimal when both variables are continuous and approximately normally distributed. In practice, researchers work with binary variables (yes/no, male/female), ordinal variables (Likert scales, education levels), and nominal categorical variables (colour, country), and each needs its own method. Pearson on such variables can be biased, attenuated, or simply meaningless.

smartcor selects a method from the detected variable types and from a few diagnostic tests run on the data. This vignette explains the statistical theory behind each method, why it is appropriate for a given variable-type combination, and what can go wrong when the wrong method is used. The decision framework article on the package website walks through the same rules with one example per pair.

Variable Types

smartcor first puts each variable into one of five types:

Type	Description	Examples
Continuous	Numeric with many distinct values	Income, temperature, test scores
Count	Non-negative integers with many distinct values	Number of citations, visits, or events
Binary	Exactly two distinct values	Male/female, yes/no, 0/1
Ordinal	Ordered categories, discrete	Likert scales (1–5), education level
Categorical	Unordered categories	Colour, country, blood type

The detect_type() function applies the following heuristic, in order:

1. Any variable with exactly 2 unique values is binary, whatever its class. This includes a two-level ordered factor.
2. Ordered factors with more than two levels are ordinal.
3. Unordered factors and character vectors with more than two unique values are categorical.
4. Logical vectors are binary.
5. Numeric vectors with 10 or fewer unique values are ordinal (configurable via ordinal_threshold).
6. Remaining non-negative, integer-valued numeric vectors are count (set detect_count = FALSE to skip this step).
7. All other numeric vectors are continuous.

So a 7-point Likert scale stored as integers (1, 2, …, 7) is detected as ordinal, which is how it should be treated. Count variables follow the continuous rules everywhere, diagnostics included; the separate label is only descriptive.

Step 0: Diagnostics

Before choosing a method, smart_cor() runs a small set of diagnostics on the pair. Which ones run depends on the variable types. All tests use the significance level alpha (default 0.05); the outlier rule has no significance level.

Outer-fence rule. detect_outliers() flags values strictly below $Q_1 - 3\,\mathrm{IQR}$ or strictly above $Q_3 + 3\,\mathrm{IQR}$, with type-7 quartiles. These are Tukey’s outer fences, so only far-out values are flagged. The rule is applied to both variables of a continuous-continuous pair and to the continuous variable of a continuous-ordinal pair. For continuous-binary and continuous-categorical pairs it is applied to the within-group residuals $x_i - \bar{x}_{g(i)}$, so that a difference between group means is not mistaken for outliers. It is never applied to an ordinal variable. When the IQR is zero (more than half of the values are tied) the rule is degenerate: nothing is flagged and the repeated values need separate inspection.

Mardia’s test checks bivariate normality for two continuous variables, through MVN::mardia(). It rejects when either the multivariate skewness or the multivariate kurtosis test rejects. Mardia’s test does not change the selection; it only qualifies Pearson’s inference (see below).

Jarque-Bera test. For one continuous variable, $$JB = \frac{n}{6}\left(S^2 + \frac{(K - 3)^2}{4}\right),$$ where $S$ and $K$ are the moment estimators of skewness and kurtosis with divisor $n$. Under normality $JB$ is approximately $\chi^2_2$. smartcor uses moments::jarque.test(). It is applied to the continuous variable of a continuous-ordinal pair and to the within-group residuals of continuous-binary and continuous-categorical pairs.

The $G^2$ threshold-model test. The latent-variable methods (biserial, polyserial, polychoric, tetrachoric) assume that the observed categories come from cutting a bivariate normal distribution at fixed thresholds. test_bivariate_normality() fits that model to the contingency table and compares observed and expected counts with the likelihood-ratio statistic $G^2 = 2 \sum_{ij} O_{ij} \log(O_{ij} / E_{ij})$, on $m_1 m_2 - m_1 - m_2$ degrees of freedom (Joreskog, 2005). A continuous variable is first binned into at most 10 quantile groups. For a continuous-binary pair the degrees of freedom are $m_1 - 2$. A 2 x 2 table is saturated (df = 0), so a binary-binary pair cannot be tested. This is the test behind assume_latent_normal = "auto".

Linearity. test_linearity() regresses each variable on the other, orders the residuals by the regressor, and checks them; residuals at tied regressor values are averaged. Box-Pierce, Ljung-Box, and BDS look for structure left in the ordered residuals; RESET looks for omitted nonlinear terms; Breusch-Pagan checks for heteroskedasticity. A direction rejects linearity when at least one structure test and RESET both reject; its p-value is $\max(\min(p_{BP}, p_{LB}, p_{BDS}), p_{RESET})$. The pair is treated as nonlinear only when both directions reject, and the combined p-value is the larger of the two directional p-values. The regression model assumes errors independent of the regressor, which can hold in only one direction unless the pair is bivariate normal, so a skewed margin alone can make one direction reject. Requiring both directions also makes the result independent of argument order.

test_bivariate_normality() picks the test from the types: Mardia for two continuous variables, Jarque-Bera plus $G^2$ for continuous-ordinal and continuous-binary pairs, and $G^2$ alone for ordinal and binary pairs. The results of every diagnostic that ran are returned in the diagnostics element of a smart_cor() result.

Correlation Methods by Variable-Type Combination

Continuous + Continuous: Pearson or Spearman

The Pearson correlation coefficient between continuous variables $X$ and $Y$ is defined as:

$$r = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^n (x_i - \bar{x})^2 \cdot \sum_{i=1}^n (y_i - \bar{y})^2}}$$

Why Pearson? For two continuous variables, Pearson is the maximum-likelihood estimator of the population correlation $\rho$ under bivariate normality. It captures linear association and has well-understood sampling properties (exact $t$-test for $H_0: \rho = 0$, Fisher’s $z$-transformation for confidence intervals).

When Pearson is not the best choice. Pearson is sensitive to outliers and non-linear relationships. When the relationship is monotonic but not linear, or when a few extreme values drive the estimate, Spearman’s rank correlation is more appropriate.

The rule is:

1. If the outer fence flags either variable, Spearman is selected.
2. Otherwise test_linearity() runs. If it rejects linearity, Spearman is selected.
3. Otherwise (linear, or the test is inconclusive) Pearson is selected.

When Pearson is selected and Mardia’s test rejects bivariate normality, Pearson is kept, because the relationship is linear and free of outliers. The rationale then says that the confidence interval and p-value are approximate and could be affected by non-normality. When the linearity rule finds residual structure or heteroskedasticity without nonlinearity, it recommends robust inference and bootstrap = "auto" replaces the Fisher-$z$ interval with a paired bootstrap interval. Kendall’s tau is always available as an alternative via method = "kendall".

Continuous + Binary: Biserial, Point-Biserial, or Rank-Biserial

The point-biserial correlation measures the association between a continuous variable and a 0/1 variable:

$$r_{pb} = \frac{\bar{X}_1 - \bar{X}_0}{s_X} \cdot \sqrt{\frac{n_1 n_0}{n^2}}$$

where $\bar{X}_1$ and $\bar{X}_0$ are the means of the continuous variable in each binary group, $s_X$ is the standard deviation of $X$ (divisor $n$), and $n_1$, $n_0$ are the group sizes.

The point-biserial correlation is identical to the Pearson correlation when the binary variable is coded as 0 and 1. This is an exact algebraic identity, not an approximation. The proof follows from substituting $Y \in \{0, 1\}$ into the Pearson formula and simplifying (see Appendix in Harshvardhan and Ranjan, 2026). Monte Carlo simulations confirm this: the maximum absolute difference between Pearson and point-biserial across 500 replications, multiple sample sizes, and multiple true correlations is exactly 0.

Biserial correlation. If the binary variable is a dichotomised version of a normal latent variable $Z^*$ (pass/fail from a test score, for example), the point-biserial coefficient understates the correlation between $X$ and $Z^*$. The biserial correlation estimates that latent correlation under the model that $(X, Z^*)$ is bivariate normal and $B = \mathbf{1}(Z^* > \tau)$. The classical formula is $r_b = r_{pb} \sqrt{pq} / \phi(z_p)$, where $p$ and $q$ are the group proportions and $\phi(z_p)$ is the normal density at the threshold. smartcor instead fits it as the polyserial correlation with a two-category ordinal variable, using polycor::polyserial() with its default two-step estimator; we recommend polyserial_ml = TRUE, which fits the full maximum-likelihood estimate. When latent normality holds the two agree closely.

Rank-biserial correlation. When the continuous variable has outliers or is far from normal within the groups, both the biserial and the point-biserial coefficients can be distorted. The rank-biserial correlation (Glass, 1966) uses only the ranks; it is described under Binary + Ordinal below.

The rule is:

1. If the outer fence flags the within-group residuals, or Jarque-Bera rejects their normality, rank-biserial is selected.
2. Otherwise the $G^2$ test runs on the binned continuous variable against the binary one. If it does not reject, biserial is selected.
3. If it rejects, point-biserial is selected.

assume_latent_normal = TRUE gives biserial and FALSE gives point-biserial, but only after step 1.

Binary + Binary: Phi, Tetrachoric, or Yule’s Q

When both variables are binary, the phi coefficient ($\phi$) measures their association. For a $2 \times 2$ contingency table:

$$\phi = \frac{ad - bc}{\sqrt{(a+b)(c+d)(a+c)(b+d)}}$$

where $a, b, c, d$ are the cell counts.

Phi = Pearson for 0/1 coding. Like the point-biserial case, $\phi$ is algebraically identical to the Pearson correlation when both variables are coded 0/1. Simulations confirm that the maximum difference is on the order of $10^{-14}$ (floating-point precision).

The attenuation problem. Phi measures the observed association between two binary variables. Often, though, each binary variable is a cut-off version of an underlying continuous (latent) variable. Pass/fail, for example, is a continuous test score cut at a threshold.

In this case, phi underestimates the latent correlation $\rho$. The attenuation depends on the marginal proportions (how far the cut points are from the median):

• When $\rho = 0.7$ and variables are dichotomised at the median, phi is approximately 0.45–0.50.
• The attenuation increases as cut points move away from the median.

The tetrachoric correlation corrects for this.

The tetrachoric correlation assumes that each observed binary variable arises from dichotomising an underlying continuous normal variable at some threshold. It estimates the correlation between the latent continuous variables, not the observed binary ones.

The model is:

$$(X^*, Y^*) \sim \text{Bivariate Normal}(\mu, \Sigma) \quad \text{with} \quad \Sigma = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}$$

The observed binary variables are $X = \mathbf{1}(X^* > \tau_x)$ and $Y = \mathbf{1}(Y^* > \tau_y)$, where $\tau_x, \tau_y$ are the thresholds estimated from the marginal proportions.

Use tetrachoric:

• When the binary variables are genuinely dichotomised continuous variables (e.g., pass/fail from test scores).
• When you want to recover the latent correlation for use in structural equation modelling or factor analysis.

Do not use tetrachoric:

• When the binary variable is inherently binary (e.g., male/female, alive/dead) with no meaningful continuous latent variable.
• When the latent normality assumption is implausible.

Yule’s Q. For two inherently binary variables, Yule’s Q $= (ad - bc)/(ad + bc)$ is a function of the odds ratio alone, so it does not depend on the marginal proportions. It is the 2 x 2 case of Goodman-Kruskal gamma.

A 2 x 2 table is always saturated (df = 0), so the $G^2$ test cannot run and the data cannot settle the latent-normality question. The rule is therefore set by assume_latent_normal:

1. "auto" (the default) gives phi.
2. TRUE gives tetrachoric.
3. FALSE gives Yule’s Q.

Continuous + Ordinal: Polyserial Correlation

The polyserial correlation is the continuous-ordinal analogue of polychoric. It assumes the ordinal variable arises from discretising a latent normal variable, while the continuous variable is observed directly. The model estimates $\rho$ between the continuous variable and the latent continuous variable underlying the ordinal one.

Simulations show polyserial is essentially unbiased (bias < 0.005) when the latent normality assumption holds. When the assumption is violated, a rank correlation is safer.

The rule is:

1. If the outer fence flags the continuous variable, or Jarque-Bera rejects its normality, a rank correlation is selected: Kendall’s tau-b if the ordinal variable has $m \le 5$ categories, Spearman otherwise.
2. Otherwise the $G^2$ test runs on the binned continuous variable against the ordinal one. If it does not reject, polyserial is selected.
3. If it rejects, Kendall’s tau-b ($m \le 5$) or Spearman is selected.

Ordinal + Ordinal: Polychoric Correlation

The polychoric correlation generalises the tetrachoric to ordinal variables with more than two categories. It assumes each observed ordinal variable is a discretisation of a latent normal variable:

$$X^* \sim N(0, 1), \quad X = k \iff \tau_{k-1} < X^* \leq \tau_k$$

The polychoric correlation $\hat{\rho}$ is the maximum-likelihood estimate of $\rho$ in the bivariate normal model, given the observed contingency table and estimated thresholds.

Advantages:

• Recovers the latent correlation without attenuation.
• Simulations show that polychoric is essentially unbiased (bias < 0.003) when the latent normality assumption holds, even with as few as 100 observations.
• More ordinal categories improve precision: the biggest improvement is from 2 to 3 categories, with diminishing returns beyond 5.

Disadvantages:

• Sensitive to violated normality. When the latent distribution is bimodal or heavily skewed, polychoric shows positive bias (up to +0.06 in simulations). The direction of bias is towards overestimation.
• Requires iterative optimisation; can fail to converge with very small samples or sparse tables.

Alternative: Kendall’s tau-b ($\tau_b$).

Kendall’s tau is defined in terms of concordant and discordant pairs:

$$\tau_b = \frac{C - D}{\sqrt{(C + D + T_X)(C + D + T_Y)}}$$

where $C$ = concordant pairs, $D$ = discordant pairs, $T_X$ = pairs tied on $X$ only, $T_Y$ = pairs tied on $Y$ only.

Simulations show that Kendall’s tau has a consistent negative bias relative to the true $\rho$ (it underestimates), but this bias is stable across all distribution types: normal, skewed, heavy-tailed, bimodal, and uniform. Goodman-Kruskal gamma, $(C - D)/(C + D)$, ignores ties altogether and is available as a further alternative.

The rule is:

1. If the $G^2$ test does not reject, polychoric is selected.
2. Otherwise Kendall’s tau-b is selected if either variable has $m \le 5$ categories, and Spearman if both have more.

Binary + Ordinal: Polychoric or Rank-Biserial

A binary variable can be treated as an ordinal variable with two categories, so the polychoric model applies here too. When it does not fit, the rank-biserial correlation (Glass, 1966) is used. It compares the mean ranks of the ordinal variable in each binary group:

$$r_{rb} = \frac{2(\bar{R}_1 - \bar{R}_0)}{n}$$

This is a nonparametric measure that respects the ordinal nature of one variable and the dichotomous nature of the other. It equals Cliff’s delta between the two groups.

The rule is:

1. If the $G^2$ test (df = $m - 2$) does not reject, polychoric is selected.
2. Otherwise rank-biserial is selected.

Kendall’s tau and Spearman are available as alternatives.

Continuous + Categorical: Correlation Ratio or Kruskal-Wallis $\eta^2_H$

A nominal variable has no order, so no signed correlation exists. The natural measure is the share of the variance of the continuous variable explained by the groups. That share is the correlation ratio squared:

$$\eta^2 = \frac{SS_{\text{between}}}{SS_{\text{total}}} = \frac{\sum_g n_g (\bar{x}_g - \bar{x})^2}{\sum_i (x_i - \bar{x})^2}$$

It is the $R^2$ of a one-way ANOVA and lies between 0 and 1. It reduces to the squared point-biserial correlation when there are two groups.

$\eta^2$ depends on group means and variances, so outliers and heavy tails affect it. The Kruskal-Wallis version replaces the values by their ranks. With the Kruskal-Wallis statistic $H$, $k$ groups, and $n$ observations,

$$\eta^2_H = \frac{H - k + 1}{n - k}.$$

It is reported as computed, so it can be slightly negative when the groups differ by less than chance would suggest.

The rule is:

1. If the outer fence flags the within-group residuals, or Jarque-Bera rejects their normality, $\eta^2_H$ is selected.
2. Otherwise $\eta^2$ is selected.

Cramer’s V on the continuous variable binned into quintiles is available via method = "cramers_v"; it discards the order and spacing of the continuous values.

Ordinal + Categorical: Freeman’s Theta

Cramer’s V treats an ordinal variable as nominal and loses its order. Freeman’s theta (Freeman, 1965) keeps the order. For each pair of groups $g$ and $h$, let $A_{gh}$ be the number of pairs of observations (one from each group) in which the member of $g$ has the higher ordinal value, and $B_{gh}$ the number in which the member of $h$ is higher. Then

$$\theta = \frac{\sum_{g < h} |A_{gh} - B_{gh}|}{\sum_{g < h} n_g n_h}.$$

Tied pairs count only in the denominator. Theta is 0 when every group has the same distribution on the ordinal variable and 1 when the groups are perfectly separated. It is the default for every ordinal-categorical pair; Cramer’s V is the alternative.

Binary + Categorical and Categorical + Categorical: Cramer’s V

When both variables are nominal (unordered categories), there is no concept of “direction” in the association. Cramer’s V is the standard symmetric measure:

$$V = \sqrt{\frac{\chi^2}{n \cdot (k - 1)}}$$

where $\chi^2$ is the Pearson chi-squared statistic and $k = \min(r, c)$ is the smaller of the number of rows and columns in the contingency table.

Cramer’s V ranges from 0 (no association) to 1 (perfect association) and is corrected for table dimensions. A binary variable paired with a categorical one is treated as a 2-level nominal variable, so Cramer’s V is used there as well.

Alternatives for two categorical variables:

• Theil’s U (uncertainty coefficient): An asymmetric measure based on information theory. $U(Y|X)$ measures the proportion of entropy in $Y$ that is explained by knowing $X$.
• Tschuprow’s T: Similar to Cramer’s V but uses the geometric mean of $(r-1)$ and $(c-1)$ instead of the minimum, making it stricter for non-square tables.

Why Not Always Use Pearson?

Why not use Pearson for everything? It depends on what you want the correlation to measure and what you will do with it.

1. Attenuation bias. When continuous variables are observed only as ordinal categories (e.g., Likert scales), Pearson/phi/point-biserial measures the association between the observed (discretised) values. This systematically underestimates the latent association. For example, with a true $\rho = 0.7$ and binary dichotomisation, phi gives approximately 0.45–0.50. Polychoric/tetrachoric recovers the correct 0.70.

This matters for factor analysis and structural equation modelling, where underestimated correlations lead to underestimated factor loadings and distorted model fit.

2. Scale appropriateness. Pearson assumes interval-scale data. Ordinal variables violate this: the “distance” between categories 1 and 2 is not necessarily the same as between 4 and 5. Rank-based measures (Kendall, Spearman) and latent-variable measures (polychoric) respect the ordinal scale.

3. Meaninglessness for nominal data. Computing Pearson on arbitrarily coded nominal variables (e.g., red = 1, blue = 2, green = 3) is meaningless because the numbers have no inherent order. Cramer’s V, $\eta^2$, and Freeman’s theta are built for nominal variables.

The Latent Normality Question

The choice between latent-variable methods (biserial, polyserial, polychoric, tetrachoric) and observed-scale methods (point-biserial, Kendall, Spearman, rank-biserial) depends on whether the underlying latent distribution is approximately normal.

When latent normality is reasonable:

• Psychological constructs measured on Likert scales (e.g., agree/disagree).
• Test scores dichotomised into pass/fail.
• Any variable that is conceptually continuous but observed in discrete categories.

When latent normality is doubtful:

• Inherently discrete variables (e.g., number of children).
• Variables with known multimodal or heavily skewed latent distributions.
• When you have no theoretical reason to assume normality.

In practice:

• If latent normality holds, the latent-variable methods are essentially unbiased and more efficient.
• If normality is violated, Kendall’s tau is a safe choice: it consistently underestimates the true $\rho$ but by a stable, predictable amount across all distribution types.
• When in doubt, compute both and compare with compare_methods(). If the estimates differ substantially, the normality assumption may be violated.

The default assume_latent_normal = "auto" lets the $G^2$ test decide, so latent-variable methods are selected whenever the data do not reject the threshold model. TRUE and FALSE fix the answer instead; the outlier and Jarque-Bera checks still run. Passing NULL asks interactively; in non-interactive sessions that falls back to TRUE. The smart_cor_df() convenience function uses "auto" and never prompts, so it is safe to use in scripts.

Summary of Method Selection Logic

The following table summarises smartcor’s decision matrix (Table 3 in Harshvardhan and Ranjan, 2026):

Variable types	Default method	Alternative
Correlation methods (signed, from -1 to +1)
Continuous + continuous	Pearson (linear, no outliers)	Spearman (outliers or nonlinear), Kendall’s tau
Continuous + binary	Biserial (if latent normality holds)	Point-biserial; rank-biserial (outliers or non-normal)
Binary + binary	Pearson = phi	Tetrachoric (latent normal), Yule’s Q (no latent continuum)
Continuous + ordinal	Polyserial (if latent normality holds)	Kendall’s tau ($m \le 5$) or Spearman
Ordinal + ordinal	Polychoric (if latent normality holds)	Kendall’s tau ($m \le 5$) or Spearman; Goodman-Kruskal gamma
Binary + ordinal	Polychoric (if latent normality holds)	Rank-biserial; Kendall’s tau, Spearman
Association measures (unsigned, from 0 to 1)
Continuous + categorical	Correlation ratio $\eta^2$	Kruskal-Wallis $\eta^2_H$ (outliers or non-normal); Cramer’s V (binned)
Binary + categorical	Cramer’s V	–
Ordinal + categorical	Freeman’s theta	Cramer’s V
Categorical + categorical	Cramer’s V	Theil’s U, Tschuprow’s T

$m$ is the number of categories of the ordinal variable (the smaller of the two for an ordinal pair). Count variables are treated as continuous, so each count combination follows the corresponding continuous row. The conditions in parentheses are the diagnostics of Step 0: the outer-fence rule, Jarque-Bera, the $G^2$ test, and test_linearity().

References

• Cramer, H. (1946). Mathematical Methods of Statistics. Princeton University Press.
• Freeman, L. C. (1965). Elementary Applied Statistics for Students in Behavioral Science. Wiley.
• Glass, G. V. (1966). Note on rank biserial correlation. Educational and Psychological Measurement, 26(3), 623–631.
• Harshvardhan, M. and Ranjan, P. (2026). smartcor: Intelligent Correlation Method Selection for Mixed Variable Types. arXiv preprint arXiv:2607.22285. https://doi.org/10.48550/arXiv.2607.22285
• Jarque, C. M. and Bera, A. K. (1980). Efficient tests for normality, homoscedasticity and serial independence of regression residuals. Economics Letters, 6(3), 255–259.
• Joreskog, K. G. (2005). Structural equation modeling with ordinal variables using LISREL. Scientific Software International.
• Kendall, M. G. (1938). A new measure of rank correlation. Biometrika, 30(1/2), 81–93.
• Kruskal, W. H. and Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis. Journal of the American Statistical Association, 47(260), 583–621.
• Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519–530.
• Olsson, U. (1979). Maximum likelihood estimation of the polychoric correlation coefficient. Psychometrika, 44(4), 443–460.
• Olsson, U., Drasgow, F., and Dorans, N. J. (1982). The polyserial correlation coefficient. Psychometrika, 47(3), 337–347.
• Pearson, K. (1900). Mathematical contributions to the theory of evolution. VII. On the correlation of characters not quantitatively measurable. Philosophical Transactions of the Royal Society A, 195, 1–47.
• Tomczak, M. and Tomczak, E. (2014). The need to report effect size estimates revisited. Trends in Sport Sciences, 21(1), 19–25.
• Tukey, J. W. (1977). Exploratory Data Analysis. Addison-Wesley.