Researchers compute correlations all the time, but the right method depends on the types of the two variables. Most software defaults to Pearson, which suits two continuous variables with a linear relationship. For binary, ordinal, or categorical variables, other methods (point-biserial, phi, polychoric, polyserial, Cramer’s V, Kendall’s tau, and others) are better suited.
Published research shows how often this is ignored. For the accompanying paper (Harshvardhan and Ranjan, 2026), we traced correlation usage through the replication supplements of published economics articles: across 197 papers issuing at least one correlation command, 92.3% of the calls were Pearson, and 88% of the papers used Pearson and nothing else, whatever the variable types involved. We also recomputed 3,322 variable pairs from those supplements. The type-appropriate method differs from the authors’ choice in 8.1% of the pairs with no ordinal variable, and in all 732 pairs that involve an ordinal variable. Where the methods differ, the median change in the coefficient is small (\(|\Delta r| \approx 0.017\)), but it goes as high as \(0.47\), and in 36 pairs the change moves the result across \(p < 0.05\).
smartcor does this for you. It detects the variable types, runs a few diagnostics, picks a method, and tells you why it picked it.
smartcor is available on CRAN:
A Python implementation, pysmartcor, is available
on PyPI (pip install pysmartcor) and is kept at feature
parity with the R package.
The main function is smart_cor(). You pass it two
vectors:
# continuous + continuous: Pearson unless there are outliers or nonlinearity
smart_cor(mtcars$mpg, mtcars$wt)
#>
#> ── Variable types ──
#>
#> ℹ Detected as continuous (numeric with 25 unique values).
#> mtcars$mpg: "continuous"
#> ℹ Detected as continuous (numeric with 29 unique values).
#> mtcars$wt: "continuous"
#>
#> ── Method selection ──
#>
#> ✔ Selected: Pearson Correlation
#>
#> ── Smart Correlation ───────────────────────────────────────────────────────────
#> Estimate: -0.8677
#> Method: Pearson Correlation
#> Variables: mtcars$mpg ("continuous") x mtcars$wt ("continuous")
#> N: 32
#> p-value: < 0.001
#> Test: t-test on r (cor.test)
#> H0: rho = 0
#> Small p-values (e.g., p < 0.05) indicate evidence against H0.
#> 95% CI (bootstrap): [-0.9252, -0.7971]
#> Source: Percentile bootstrap
#>
#> Diagnostics
#> Outliers (outer fence): none
#> Normality (Mardia): p = 0.0178 (rejected)
#> Linearity: p = 0.0562 (no nonlinearity detected)
#>
#> Sorted residual diagnostics (y on x): Pearson (robust/bootstrap SE) The
#> interval uses a paired bootstrap; the p-value remains the standard Pearson
#> t-test. Mardia's test rejects bivariate normality, so the confidence interval
#> and p-value are approximate and could be affected by non-normality.
#>
#> ℹ Alternatives: "spearman" and "kendall" (pass `method = "..."` to use)The rationale for this pair says that Mardia’s test rejects bivariate normality. Pearson is kept because the relationship is linear and free of outliers, but its interval and p-value are approximate.
# continuous + binary: biserial when the diagnostics support latent normality
smart_cor(mtcars$mpg, mtcars$vs)
#>
#> ── Variable types ──
#>
#> ℹ Detected as continuous (numeric with 25 unique values).
#> mtcars$mpg: "continuous"
#> ℹ Detected as binary (numeric with 2 unique values: 0, 1).
#> mtcars$vs: "binary"
#>
#> ── Method selection ──
#>
#> ℹ The LR test does not reject latent bivariate normality (p = 0.228).
#> ✔ Selected: Biserial Correlation
#>
#> ── Smart Correlation ───────────────────────────────────────────────────────────
#> Estimate: 0.8194
#> Method: Biserial Correlation
#> Variables: mtcars$mpg ("continuous") x mtcars$vs ("binary")
#> N: 32
#> p-value: < 0.001
#> Test: Wald test on Fisher-z scale
#> H0: rho = 0
#> Small p-values (e.g., p < 0.05) indicate evidence against H0.
#> 95% CI: [0.5596, 0.9326]
#> Source: Olsson, Drasgow, and Dorans (1982)
#>
#> Diagnostics
#> Outliers (outer fence): none
#> Normality (Jarque-Bera): p = 0.429 (not rejected)
#> Latent normality (LR chi-square = 10.55, df = 8): p = 0.228 (not rejected)
#>
#> No outliers and no evidence against normality of the within-group residuals.
#> The LR test does not reject latent bivariate normality (p = 0.228); Biserial
#> Correlation selected.
#>
#> ℹ Alternatives: "point_biserial" and "rank_biserial" (pass `method = "..."` to use)Before choosing a method, smart_cor() runs a few checks
on the pair. Which ones run depends on the variable types, and each is
also available as a function.
detect_outliers() applies the outer-fence rule: values
below \(Q_1 - 3\,\mathrm{IQR}\) or
above \(Q_3 + 3\,\mathrm{IQR}\) are
flagged.
detect_outliers(mtcars$mpg)
#> No outliers outside [-6.7, 44.92]
detect_outliers(c(mtcars$mpg, 120))
#> 1 outlier outside [-6.4, 44.7]: 120test_bivariate_normality() picks the normality test from
the types. Two continuous variables get Mardia’s test. A continuous
variable paired with a binary or ordinal one gets the Jarque-Bera test
plus the likelihood-ratio (\(G^2\))
test of the bivariate normal threshold model. Two ordinal or binary
variables get the \(G^2\) test
alone.
test_bivariate_normality(mtcars$mpg, mtcars$wt)
#>
#> ── Bivariate Normality Test (Mardia) ──
#>
#> H0: x and y are bivariate normal
#> N = 32
#> Mardia skewness: 11.9364, df = 4, p = 0.01783
#> Mardia kurtosis: 0.5583, p = 0.5766
#> Conclusion: reject. Evidence against H0 at alpha = 0.05 (Mardia skewness).
test_bivariate_normality(mtcars$mpg, mtcars$vs)
#>
#> ── Bivariate Normality Test (Jarque-Bera and LR chi-square) ──
#>
#> H0: the within-group residuals of the continuous variable are normal, and the
#> table arises from a bivariate normal threshold model
#> N = 32
#> Table dimensions: 10 x 2; polychoric rho: 0.8554
#> Jarque-Bera: 1.6924, df = 2, p = 0.429
#> LR chi-square: 10.5534, df = 8, p = 0.2283
#> Conclusion: fail_to_reject. No evidence against H0 at alpha = 0.05.
test_bivariate_normality(mtcars$gear, mtcars$carb)
#>
#> ── Bivariate Normality Test (LR chi-square) ──
#>
#> H0: the table arises from a bivariate normal threshold model
#> N = 32
#> Table dimensions: 3 x 6; polychoric rho: 0.2381
#> LR chi-square: 14.0043, df = 9, p = 0.1222
#> Conclusion: fail_to_reject. No evidence against H0 at alpha = 0.05.For continuous/count pairs, test_linearity() regresses
each variable on the other and orders the residuals by the regressor. A
direction rejects linearity when a residual structure test (Box-Pierce,
Ljung-Box, or BDS) and RESET both reject, and Spearman is selected only
when both directions reject.
diagnostics = test_linearity(mtcars$mpg, mtcars$wt)
diagnostics
#>
#> ── Linearity Test (sorted residual diagnostics) ──
#>
#> H0: the relationship between x and y is linear
#> N = 32, alpha = 0.05
#> y on x p = 0.05617, x on y p = 0.01263
#> Components (y on x):
#> Box-Pierce p = 0.5509, Ljung-Box p = 0.527, BDS p (min) = 0.006207
#> Breusch-Pagan p = 0.2699, RESET p = 0.05617
#> Combined p = 0.05617 (fail_to_reject)
#> Decision: Pearson (robust/bootstrap SE)
#> No evidence of nonlinearity in at least one direction at alpha = 0.05; this
#> does not establish linearity.The result does not depend on argument order. Inconclusive diagnostics retain Pearson and report the reason. Non-rejection does not prove linearity.
Every test uses the significance level alpha, which
defaults to 0.05 and is an argument of smart_cor(),
smart_cormat(), smart_cor_df(), and
compare_methods(). The tests that ran are stored in the
diagnostics element of the result; entries that did not run
are NULL.
res = smart_cor(mtcars$mpg, mtcars$vs, alpha = 0.01, verbose = FALSE)
names(res$diagnostics)
#> [1] "outliers" "normality" "latent" "linearity"
res$diagnostics$normality
#>
#> ── Bivariate Normality Test (Jarque-Bera) ──
#>
#> H0: the within-group residuals are normal
#> N = 32
#> Jarque-Bera: 1.6924, df = 2, p = 0.429
#> Conclusion: fail_to_reject. No evidence against H0 at alpha = 0.01.An explicit method skips all diagnostics.
When both variables have a small number of unique values, smartcor detects them as ordinal:
# numeric variables with 10 or fewer unique values are treated as ordinal;
# assume_latent_normal controls whether polychoric or Kendall is used
smart_cor(mtcars$gear, mtcars$carb, assume_latent_normal = FALSE)
#>
#> ── Variable types ──
#>
#> ℹ Detected as ordinal (numeric with 3 unique values, at or below threshold of 10).
#> mtcars$gear: "ordinal"
#> ℹ Detected as ordinal (numeric with 6 unique values, at or below threshold of 10).
#> mtcars$carb: "ordinal"
#>
#> ── Method selection ──
#>
#> ✔ Selected: Kendall's Tau-b
#>
#> ── Smart Correlation ───────────────────────────────────────────────────────────
#> Estimate: 0.0980
#> Method: Kendall's Tau-b
#> Variables: mtcars$gear ("ordinal") x mtcars$carb ("ordinal")
#> N: 32
#> p-value: 0.5302
#> Test: asymptotic normal (cor.test, exact = FALSE)
#> H0: tau = 0
#> Small p-values (e.g., p < 0.05) indicate evidence against H0.
#> 95% CI: [-0.1455, 0.3303]
#> Source: Fieller, Hartley, and Pearson (1957)
#>
#> Both variables are ordinal. Latent normality not assumed (assume_latent_normal
#> = FALSE); Kendall's Tau-b selected (the smaller variable has 3 categories).
#>
#> ℹ Alternatives: "polychoric", "spearman", and "gamma" (pass `method = "..."` to use)# some categorical data
set.seed(1)
colour = factor(sample(c("red", "blue", "green"), 100, replace = TRUE))
shape = factor(sample(c("circle", "square", "triangle"), 100, replace = TRUE))
smart_cor(colour, shape)
#>
#> ── Variable types ──
#>
#> ℹ Detected as categorical (factor/character with 3 levels).
#> colour: "categorical"
#> ℹ Detected as categorical (factor/character with 3 levels).
#> shape: "categorical"
#>
#> ── Method selection ──
#>
#> ✔ Selected: Cramer's V
#>
#> ── Smart Correlation ───────────────────────────────────────────────────────────
#> Estimate: 0.0945
#> Method: Cramer's V
#> Variables: colour ("categorical") x shape ("categorical")
#> N: 100
#> p-value: 0.7752
#> Test: Pearson chi-square test of independence
#> H0: variables are independent
#> Small p-values indicate evidence that the variables are not independent.
#> 95% CI: [0.0000, 0.1715]
#> Source: Noncentral chi-square pivot
#>
#> Both variables are categorical. Cramer's V selected as the standard symmetric
#> measure of nominal association.
#>
#> ℹ Alternatives: "theils_u" and "tschuprows_t" (pass `method = "..."` to use)A continuous variable paired with a categorical one gets the correlation ratio \(\eta^2\), or its Kruskal-Wallis version when the within-group residuals have outliers or are far from normal:
smart_cor(mtcars$mpg, factor(mtcars$cyl))
#>
#> ── Variable types ──
#>
#> ℹ Detected as continuous (numeric with 25 unique values).
#> mtcars$mpg: "continuous"
#> ℹ Detected as categorical (factor/character with 3 levels).
#> factor(mtcars$cyl): "categorical"
#>
#> ── Method selection ──
#>
#> ✔ Selected: Correlation Ratio (eta-squared)
#>
#> ── Smart Correlation ───────────────────────────────────────────────────────────
#> Estimate: 0.7325
#> Method: Correlation Ratio (eta-squared)
#> Variables: mtcars$mpg ("continuous") x factor(mtcars$cyl) ("categorical")
#> N: 32
#> p-value: < 0.001
#> Test: One-way ANOVA F test
#> H0: eta^2 = 0 (equal group means)
#> Small p-values (e.g., p < 0.05) indicate evidence against H0.
#> 95% CI: [0.5362, 0.8286]
#> Source: Noncentral F pivot
#>
#> Diagnostics
#> Outliers (outer fence): none
#> Normality (Jarque-Bera): p = 0.859 (not rejected)
#>
#> No outliers and no evidence against normality of the within-group residuals;
#> Correlation Ratio (eta-squared) selected.
#>
#> ℹ Alternatives: "kw_eta_squared" and "cramers_v" (pass `method = "..."` to use)If the detected type is wrong, you can set it yourself:
# force cyl to be treated as continuous
smart_cor(mtcars$mpg, mtcars$cyl, y_type = "continuous")
#>
#> ── Variable types ──
#>
#> ℹ Detected as continuous (numeric with 25 unique values).
#> mtcars$mpg: "continuous"
#> mtcars$cyl: "continuous" (user-specified)
#>
#> ── Method selection ──
#>
#> ✔ Selected: Pearson Correlation
#>
#> ── Smart Correlation ───────────────────────────────────────────────────────────
#> Estimate: -0.8522
#> Method: Pearson Correlation
#> Variables: mtcars$mpg ("continuous") x mtcars$cyl ("continuous")
#> N: 32
#> p-value: < 0.001
#> Test: t-test on r (cor.test)
#> H0: rho = 0
#> Small p-values (e.g., p < 0.05) indicate evidence against H0.
#> 95% CI: [-0.9258, -0.7163]
#> Source: Fisher z (cor.test)
#>
#> Diagnostics
#> Outliers (outer fence): none
#> Normality (Mardia): p = 0.216 (not rejected)
#> Linearity: p = NA (inconclusive)
#>
#> Linearity test inconclusive: nonlinearity was found with one variable as
#> regressor, and the other direction could not be tested.
#>
#> ℹ Alternatives: "spearman" and "kendall" (pass `method = "..."` to use)smart_cor(mtcars$mpg, mtcars$wt, method = "spearman")
#>
#> ── Variable types ──
#>
#> ℹ Detected as continuous (numeric with 25 unique values).
#> mtcars$mpg: "continuous"
#> ℹ Detected as continuous (numeric with 29 unique values).
#> mtcars$wt: "continuous"
#>
#> ── Method ──
#>
#> ℹ User specified method: "spearman"
#>
#> ── Smart Correlation ───────────────────────────────────────────────────────────
#> Estimate: -0.8864
#> Method: Spearman Rank Correlation
#> Variables: mtcars$mpg ("continuous") x mtcars$wt ("continuous")
#> N: 32
#> p-value: < 0.001
#> Test: t approximation (cor.test, exact = FALSE)
#> H0: rho_S = 0
#> Small p-values (e.g., p < 0.05) indicate evidence against H0.
#> 95% CI: [-0.9503, -0.7511]
#> Source: Bonett and Wright (2000)
#>
#> Method selected by user: spearmansmart_cormat() computes all pairwise correlations for a
data frame, selecting a method for each pair:
mat = smart_cormat(
mtcars[, c("mpg", "cyl", "vs", "gear")],
assume_latent_normal = FALSE
)
#>
#> ── Smart Correlation Matrix (4 variables) ──────────────────────────────────────
#>
#> ── Detected variable types ──
#>
#> mpg: "continuous"
#> cyl: "ordinal"
#> vs: "binary"
#> gear: "ordinal"
#>
#> ── Computing correlations ──
#>
#> ℹ Processing 6 variable pairs...
#> ✔ Done. 6 pairs computed.
mat
#>
#> ── Smart Correlation Matrix ────────────────────────────────────────────────────
#>
#> ── Variable types ──
#>
#> mpg: "continuous"
#> cyl: "ordinal"
#> vs: "binary"
#> gear: "ordinal"
#>
#> ── Correlations ──
#>
#> mpg cyl vs gear
#> mpg 1.00 -0.80 0.66 0.43
#> cyl -0.80 1.00 -0.88 -0.51
#> vs 0.66 -0.88 1.00 0.30
#> gear 0.43 -0.51 0.30 1.00
#> ── Methods used ──
#> mpg cyl vs gear
#> mpg Pr Kn PB Kn
#> cyl Kn Kn RB Kn
#> vs PB RB YQ RB
#> gear Kn Kn RB Kn
#>
#> ℹ Legend:
#> "Pr" = Pearson Correlation
#> "Kn" = Kendall's Tau-b
#> "PB" = Point-Biserial Correlation (= Pearson)
#> "RB" = Rank-Biserial Correlation
#> "YQ" = Yule's QBoth smart_cor() and smart_cormat() support
tidy(), which returns a data frame that works with
tidyverse tools. The last four columns summarise the diagnostics
(outlier_flag, normality_p,
latent_p, linearity_p) and are NA
when a test did not run.
res = smart_cor(mtcars$mpg, mtcars$wt, verbose = FALSE)
tidy(res)
#> # A tibble: 1 × 24
#> var_x var_y x_type y_type estimate method method_label rationale statistic
#> <chr> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <dbl>
#> 1 mtcars$m… mtca… conti… conti… -0.868 pears… Pearson Cor… Sorted r… -9.56
#> # ℹ 15 more variables: p.value <dbl>, p_method <chr>, null_hypothesis <chr>,
#> # p_interpretation <chr>, ci_lower <dbl>, ci_upper <dbl>, conf_level <dbl>,
#> # ci_method <chr>, ci_source <chr>, n <int>, n_complete <int>,
#> # outlier_flag <lgl>, normality_p <dbl>, latent_p <dbl>, linearity_p <dbl>tidy(mat)
#> # A tibble: 6 × 24
#> var_x var_y x_type y_type estimate method method_label rationale statistic
#> <chr> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <dbl>
#> 1 mpg cyl continuous ordin… -0.795 kenda… Kendall's T… No outli… -5.59
#> 2 mpg vs continuous binary 0.664 point… Point-Biser… No outli… 4.86
#> 3 mpg gear continuous ordin… 0.433 kenda… Kendall's T… No outli… 3.02
#> 4 cyl vs ordinal binary -0.881 rank_… Rank-Biseri… One vari… -3.91
#> 5 cyl gear ordinal ordin… -0.513 kenda… Kendall's T… Both var… -3.16
#> 6 vs gear binary ordin… 0.302 rank_… Rank-Biseri… One vari… 1.43
#> # ℹ 15 more variables: p.value <dbl>, p_method <chr>, null_hypothesis <chr>,
#> # p_interpretation <chr>, ci_lower <dbl>, ci_upper <dbl>, conf_level <dbl>,
#> # ci_method <chr>, ci_source <chr>, n <int>, n_complete <int>,
#> # outlier_flag <lgl>, normality_p <dbl>, latent_p <dbl>, linearity_p <dbl>Use available_methods() to see what methods are
available for a given variable-type pair:
available_methods("ordinal", "ordinal")
#> # A tibble: 4 × 2
#> method label
#> <chr> <chr>
#> 1 polychoric Polychoric Correlation
#> 2 kendall Kendall's Tau-b
#> 3 spearman Spearman Rank Correlation
#> 4 gamma Goodman-Kruskal's Gamma
available_methods("continuous", "binary")
#> # A tibble: 3 × 2
#> method label
#> <chr> <chr>
#> 1 point_biserial Point-Biserial Correlation (= Pearson)
#> 2 biserial Biserial Correlation
#> 3 rank_biserial Rank-Biserial Correlation
available_methods("categorical", "categorical")
#> # A tibble: 3 × 2
#> method label
#> <chr> <chr>
#> 1 cramers_v Cramer's V
#> 2 theils_u Theil's U (Uncertainty Coefficient)
#> 3 tschuprows_t Tschuprow's TThe method selection follows the decision matrix of Harshvardhan and Ranjan (2026), drawn from simulation studies comparing correlation methods:
| 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 numeric continuous variables, so each count
combination follows the corresponding continuous row. The
vignette("theory") vignette explains each rule, and the
decision framework article on the package website shows one example per
pair.
The simulation studies in the paper found that:
For ordinal and binary variables, and for a continuous variable
paired with one, the default assume_latent_normal = "auto"
runs the likelihood-ratio \(G^2\) test
of test_bivariate_normality() and follows its result:
# runs the LR test and picks polychoric or Kendall accordingly
smart_cor(mtcars$gear, mtcars$carb, verbose = FALSE)$rationale
#> [1] "Both variables are ordinal. The LR test does not reject latent bivariate normality (p = 0.122); Polychoric Correlation selected."A binary-binary pair is the exception: its 2 x 2 table is saturated
(df = 0), the test cannot run, so smartcor uses phi.
assume_latent_normal = TRUE gives tetrachoric and
FALSE gives Yule’s Q. For the other pairs,
TRUE or FALSE skips the \(G^2\) test and fixes the answer, but the
outlier and Jarque-Bera checks still run. Set the argument to
NULL to be asked interactively (non-interactive sessions
then default to TRUE).
The latent-variable methods rely on polycor
(polychoric, polyserial, biserial, and tetrachoric estimates, with their
standard errors). Mardia’s test comes from MVN and the
Jarque-Bera test from moments. All three are listed
under Imports, so they install with smartcor.
A latent-variable estimator can fail numerically on some data, for
example with a singular Hessian or a sparse contingency table. smartcor
does not stop with an error when this happens. For the polyserial
correlation fitted with polyserial_ml = TRUE it first
retries the fit from a safe starting value, since polycor’s default
start for the correlation can land outside \((-1, 1)\) on strongly correlated pairs and
abort the optimiser before it takes a step; this fixes most failures.
Only if no fit can be obtained does smartcor fall back to Kendall’s tau
(for polychoric), Spearman (for polyserial), or point-biserial (for
biserial); the returned method, method_label,
and rationale then name the method actually used, so you
can always see which method was used.
nearPD = FALSE is the default. Enable it explicitly when
a downstream procedure requires a positive-definite matrix:
adjusted = smart_cormat(mtcars[, c("mpg", "wt", "disp")], nearPD = TRUE, verbose = FALSE)
adjusted$nearPD
#> $normF
#> [1] 1.481222e-15
#>
#> $iterations
#> [1] 1
#>
#> $converged
#> [1] TRUEadjusted$correlations contains the adjusted matrix;
adjusted$original_correlations retains the original.
Methods, details, and tidy() describe the original pairwise
estimates and inference. Resolve missing or infinite entries before
adjustment. smart_cor_df(..., nearPD = TRUE) returns the
adjusted correlation table and original method table, without the extra
diagnostics.