Test Bivariate Normality for a Pair of Variables — test_bivariate

Tests the normality assumption behind Pearson's inference and behind the latent-variable correlations (polychoric, polyserial, biserial). The test depends on the variable types:

continuous + continuous: Mardia's multivariate skewness and kurtosis tests, via MVN::mardia().
continuous + ordinal: Jarque-Bera on the continuous variable (moments::jarque.test()), plus the LR chi-square threshold-model test after binning the continuous variable.
continuous + binary: Jarque-Bera on the residuals of the continuous variable around its group means, plus the LR chi-square test after binning (df = m1 - 2).
ordinal or binary pairs: The LR chi-square test. A binary-binary table is saturated (df = 0) and cannot be tested.

Count variables are treated as continuous. Categorical variables are not supported.

Usage

test_bivariate_normality(
  x,
  y = NULL,
  x_type = NULL,
  y_type = NULL,
  alpha = 0.05,
  max_levels = 10L
)

Arguments

x, y: Vectors of the same length. Alternatively, a pre-computed contingency table (matrix or table object) can be passed as x with y = NULL, which runs the LR chi-square test on that table.
x_type, y_type: Optional variable types ("continuous", "count", "binary", or "ordinal"). Detected with detect_type() when NULL.
alpha: Significance level for the test. Default: 0.05.
max_levels: Integer. Cap on the number of categories per variable when building the contingency table from numeric vectors. A numeric variable with more than max_levels distinct values (e.g. a continuous or count variable) is binned into max_levels quantile groups before the table is formed. This keeps the maximum-likelihood threshold model tractable: an unbinned continuous variable would induce one threshold per distinct value, making the ML fit prohibitively slow and the chi-square reference distribution meaningless. Genuine ordinal/binary inputs (at most max_levels distinct values) are never altered. Default: 10.

Value

An object of class "smartcor_normtest", a list with elements:

statistic: The LR chi-square statistic.
df: Degrees of freedom of the LR chi-square test.
p.value: The overall p-value (smallest component p-value).
rho: The two-step estimate of the polychoric correlation.
row_thresholds: Estimated thresholds for the first variable.
col_thresholds: Estimated thresholds for the second variable.
n: Total number of observations.
m1: Number of categories for the first variable.
m2: Number of categories for the second variable.
alpha: The significance level used.
conclusion: One of "reject", "fail_to_reject", or "saturated".
test: Name of the test, e.g. "Mardia", "Jarque-Bera and LR chi-square", or "LR chi-square".
components: Data frame with one row per component test and columns test, statistic, df, and p.value.
null_hypothesis: The null hypothesis in words.
interpretation: One plain sentence stating the result.

statistic, df, rho, the thresholds, m1, and m2 describe the LR chi-square component and are NA (or empty) when it was not run.

Details

The LR chi-square test follows Joreskog (2005). The bivariate normal threshold model is fitted with the two-step estimator (Olsson 1979) via polycor::polychor(): the thresholds come from the marginal proportions and the correlation is then estimated by one-dimensional maximum likelihood given those thresholds. The likelihood-ratio statistic \(G^2 = 2 \sum_{ij} O_{ij} \log(O_{ij} / E_{ij})\) compares the observed counts \(O_{ij}\) with the expected counts \(E_{ij}\) under the fitted model (cells with zero observed count are skipped). The degrees of freedom equal \(m_1 \times m_2 - m_1 - m_2\), where \(m_1\) and \(m_2\) are the numbers of categories of the two variables. For a \(2 \times 2\) table the model is saturated (df = 0) and no test is possible.

When several component tests are run, the overall test rejects if any component rejects at alpha, and the overall p.value is the smallest component p-value. The Python pysmartcor package computes the same statistics, so the two implementations agree.

Examples

csv = system.file("extdata", "gss_2024_casestudy.csv", package = "smartcor")
gss = read.csv(csv)

# Check the latent-normality assumption for two ordinal variables
test_bivariate_normality(gss$degree, gss$happy)
#> 
#> ── Bivariate Normality Test (LR chi-square) ──
#> 
#> H0: the table arises from a bivariate normal threshold model
#> N = 3000
#> Table dimensions: 5 x 3; polychoric rho: -0.1032
#> LR chi-square: 23.9731, df = 7, p = 0.001152
#> Conclusion: reject. Evidence against H0 at alpha = 0.05 (LR chi-square).

# Mardia's test for two continuous variables
test_bivariate_normality(gss$age, gss$coninc, x_type = "continuous",
                         y_type = "continuous")
#> 
#> ── Bivariate Normality Test (Mardia) ──
#> 
#> H0: x and y are bivariate normal
#> N = 3000
#> Mardia skewness: 666.3715, df = 4, p = < 2.2e-16
#> Mardia kurtosis: -6.6499, p = 2.934e-11
#> Conclusion: reject. Evidence against H0 at alpha = 0.05 (Mardia skewness,
#> Mardia kurtosis).