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.
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:
ordinal_threshold).detect_count = FALSE to skip
this step).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.
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.
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:
test_linearity() runs. If it rejects
linearity, Spearman 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".
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:
assume_latent_normal = TRUE gives biserial and
FALSE gives point-biserial, but only after step 1.
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):
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:
Do not use tetrachoric:
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:
"auto" (the default) gives phi.TRUE gives tetrachoric.FALSE gives Yule’s Q.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:
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:
Disadvantages:
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:
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:
Kendall’s tau and Spearman are available as alternatives.
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:
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.
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.
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:
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 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:
When latent normality is doubtful:
In practice:
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.
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().