smartcor in Python • smartcor

smartcor is available in R and in Python, and both versions have the same features. They detect the same variable types, run the same diagnostics (outer fence, Mardia, Jarque-Bera, the likelihood-ratio test of latent normality, and the linearity test), and choose the method with the same rules. Both return the estimate, confidence interval, p-value, and the reason for the choice. The numbers usually match, with a few exceptions listed at the end of this page. Use this page to translate any R example on this site into Python.

Install

The package is on PyPI:

pip install pysmartcor

Optional plotting support (matplotlib and seaborn):

pip install 'pysmartcor[viz]'

The package needs NumPy, SciPy, pandas, and statsmodels. It does not need R or any R package. The latent-variable estimators (polychoric, polyserial, biserial, tetrachoric) use our own corrected copy of the maximum-likelihood code from semopy. Mardia’s test uses the same formulas as MVN::mardia(), and the Jarque-Bera test uses scipy.stats.jarque_bera.

Function-by-function translation

R	Python
`smart_cor(x, y)`	`smart_cor(x, y)`
`smart_cormat(data)`	`smart_cormat(data)`
`smart_cor_df(data)`	`smart_cor_df(data)`
`compare_methods(x, y)`	`compare_methods(x, y)`
`detect_outliers(x)`	`detect_outliers(x)`
`test_bivariate_normality(x, y)`	`test_bivariate_normality(x, y)`
`test_linearity(x, y)`	`test_linearity(x, y)`
`detect_type(x)`	`detect_type(x)`
`available_methods(tx, ty)`	`available_methods(tx, ty)`
`tidy(result)`	`result.tidy()`
`print(result)`	`print(result)`
`plot(m)` / `ggcor_heatmap(m)`	`cor_heatmap(m)`
`ggcor_method_heatmap(m)`	`method_heatmap(m)`
`alpha = 0.01`	`alpha=0.01`
`seed = 42`	`rng = 42`
`result$p.value`	`result.p_value`
`result$diagnostics$outliers`	`result.diagnostics["outliers"]`
`lin$p.value`, `lin$null_hypothesis` (from `test_linearity()`)	`lin.p_value`, `lin.null_hypothesis`

The result classes match up as follows: smartcor_outliers is OutlierResult, smartcor_normtest is NormalityTestResult, and smartcor_lintest is LinearityTestResult. The fields have the same names, except that p.value becomes p_value (also in the components table of a normality test). detect_outliers() returns 1-based positions in index in R and 0-based positions in Python.

The tidy output has the same 24 columns in the same order in both languages. In Python they are listed in pysmartcor.TIDY_COLUMNS; the last four are outlier_flag, normality_p, latent_p, and linearity_p.

All 18 method codes are shared, including the four added in 1.2.0: biserial, eta_squared, kw_eta_squared, and freemans_theta. So are the two-letter codes used in the method heatmaps (Bs, E2, EH, Ft; tetrachoric is Tc and Tschuprow’s T is Ts).

A worked pair, side by side

library(smartcor)
gss = read.csv(system.file("extdata", "gss_2024_casestudy.csv", package = "smartcor"))
res = smart_cor(gss$coninc, gss$sex)
tidy(res)

Python:

from importlib.resources import as_file, files
import pandas as pd
from pysmartcor import smart_cor

with as_file(files("pysmartcor").joinpath("data/gss_2024_casestudy.csv")) as csv:
    gss = pd.read_csv(csv)
res = smart_cor(gss["coninc"], gss["sex"])
res.tidy()

Both detect the same types and run the same diagnostics. Income residuals within each sex fail the Jarque-Bera test, so both choose the rank-biserial correlation and report the same estimate, interval, and p-value.

Diagnostics in Python

from pysmartcor import detect_outliers, test_bivariate_normality, test_linearity

detect_outliers(gss["coninc"])
test_bivariate_normality(gss["coninc"], gss["degree"])
lin = test_linearity(gss["coninc"], gss["age"])
lin.p_value, lin.decision

res = smart_cor(gss["coninc"], gss["sex"], alpha=0.01, verbose=False)
res.diagnostics["normality"]
res.outlier_flag, res.normality_p, res.latent_p, res.linearity_p

res.diagnostics is a dictionary with the keys outliers, normality, latent, and linearity; a test that did not run is None. For two continuous variables the outliers entry is a dictionary with "x" and "y". smart_cormat() does not keep diagnostics for each pair, but it keeps the four summary fields, so to_long() still includes them.

Changes in 1.2.0

from pysmartcor import smart_cormat, smart_cor_df

result = smart_cor(gss["coninc"], gss["age"], verbose=False)
print(result.linearity_test)

# Adjustment is off unless explicitly enabled.
adjusted = smart_cormat(gss[["coninc", "age", "sex"]], nearPD=True, verbose=False)
original = adjusted.original_correlations
print(adjusted.nearPD)
tables = smart_cor_df(gss[["coninc", "age", "sex"]], nearPD=True)

In both languages the rules now follow Table 3 of the paper. The default method changes for eight of the ten type pairs, and the changelog lists them. Each rule starts with the Step 0 diagnostics: the outer-fence rule, Mardia’s test for two continuous variables, the Jarque-Bera test for one continuous variable, and the likelihood-ratio test of latent normality.

The linearity test runs in both directions. It orders the residuals by the regressor, averages them at tied values, and runs Box-Pierce, Ljung-Box, BDS, Breusch-Pagan, and RESET. For two continuous or count variables without outliers, Spearman is chosen only when both directions reject linearity; otherwise Pearson is kept. Each direction’s p-value is max(min(p_bp, p_lb, p_bds), p_reset), and the combined p_value is the larger of the two. The BDS p-value is the smallest over the distance thresholds. The result does not depend on argument order. If the tests are inconclusive, Pearson is kept and the rationale says why. Passing a method skips the rule.

When the rule recommends robust inference for Pearson, bootstrap="auto" gives a paired-bootstrap interval. The p-value is still the standard Pearson test. Use bootstrap=False to turn resampling off.

Both matrix functions still default to nearPD=False. With nearPD=True, correlations holds the adjusted matrix, while the methods, details, and to_long() describe the original estimates. Only smart_cormat() keeps the original matrix and the adjustment details. Remove missing or infinite entries before adjusting.

What matches and what differs

Mardia’s test uses the formulas of MVN::mardia(), including its small-sample correction. The Jarque-Bera test uses the same moment estimators as moments::jarque.test(). We test these and the linearity diagnostics against results from R, but only on our test inputs, so we cannot promise identical numbers for every function and dataset.

Each language uses its own random number generator. A seed makes the bootstrap and permutation results repeatable within a language, but R and Python will not give the same draws. This affects the intervals of Kruskal-Wallis $\eta^2_H$ , Freeman’s theta, and Theil’s U, and the permutation p-values of Freeman’s theta and Theil’s U.

The latent-variable estimates (polychoric, polyserial, biserial, tetrachoric) can differ in the last decimals when the correlation is close to $\pm 1$ . R fits rho and the thresholds together by maximum likelihood with polycor, while Python fits rho with the thresholds fixed at their marginal estimates.