Skip to contents

Weighted canonical correlation for categorical variables

weighted_cancor.Rd

Compute the first canonical correlation coefficient between two categorical variables via their indicator (dummy) sets. The coefficient equals the largest singular value of the standardized cross-covariance between the indicator sets and lies in [0, 1]. For a \(2 \times 2\) table, this reduces to the absolute \(\phi\) coefficient and Cramer's \(V\).

Usage

weighted_cancor(x, y, w = NULL)

weighted_cancor_sq(x, y, w = NULL)

Arguments

x

(factor)
The first categorical variable (with at least two observed levels). Observations with missing values are excluded pairwise. Each distinct level defines one indicator in the dummy set used for canonical correlation. Non-factor vectors will be coerced to factor.

y

(factor)
The second categorical variable (with at least two observed levels). Observations with missing values are excluded pairwise. Each distinct level defines one indicator in the dummy set used for canonical correlation. Non-factor vectors will be coerced to factor.

w

(numeric or NULL)
Weights used to compute the weighted canonical correlation. If numeric, values must be nonnegative and the same length as x and y. Zero weights are allowed. Observations with missing, negative, or non-finite values are excluded pairwise. Set to NULL (default) for an unweighted computation.

Value

  • weighted_cancor(): Scalar numeric in \([0, 1]\) or NA_real_.

  • weighted_cancor_sq(): Scalar numeric \([0, 1]\) or NA_real_.

Returns NA_real_ if all observations are removed during filtering or the total weight is nonpositive. Returns 0 if all mass lies in a single level of either variable or if numerical degeneracy prevents computation.

Details

Let \(X \in \{1,\ldots,r\}\) and \(Y \in \{1,\ldots,s\}\) denote categorical variables, and let \(w_k \ge 0\) be case weights (internally normalized to sum to one). Denote the weighted joint distribution by

$$P = \{p_{ij}\}_{i=1..r,\,j=1..s}$$

with

$$p_{ij} = \sum_{k} w_k \mathbb{1}(x_k=i, y_k=j)$$

and marginals

$$p_x = P \mathbf{1}_s$$ $$p_y = P^\top \mathbf{1}_r.$$

The covariance blocks for the indicator sets are

$$S_{xx} = \mathrm{diag}(p_x) - p_x p_x^\top$$ $$S_{yy} = \mathrm{diag}(p_y) - p_y p_y^\top$$ $$S_{xy} = P - p_x p_y^\top.$$

The first canonical correlation is then

$$\rho_1 = \sigma_{\max}(S_{xx}^{-1/2} \, S_{xy} \, S_{yy}^{-1/2}),$$

where \(\sigma_{\max}(\cdot)\) is the largest singular value and the inverses are taken on the support of the positive eigenvalues. This \(\rho_1\) equals the dominant singular value in correspondence analysis of the contingency table built from \(X\) and \(Y\), i.e., the square root of the largest principal inertia.

The helper weighted_cancor_sq() returns \(\rho^2_1\), the weighted squared first canonical correlation coefficient, which takes values in \([0,1]\).

Examples

#----------------------------------------------------------------------------
# weighted_cancor() examples
#----------------------------------------------------------------------------
library(vclust)

# 2x2 example: equals |phi| and Cramer's V
x <- factor(rep(c("A", "B"), times = c(30, 70)))
y <- factor(c(rep("U", 25), rep("V", 5),  # mostly A-U
              rep("U", 20), rep("V", 50))) # mostly B-V
weighted_cancor(x, y)
#> [1] 0.5044296

# Compare to |phi| computed from weighted proportions
tab <- prop.table(table(x, y))
px <- rowSums(tab)
py <- colSums(tab)
phi <- (tab[1,1]*tab[2,2] - tab[1,2]*tab[2,1]) / sqrt(px[1]*px[2]*py[1]*py[2])
phi <- as.numeric(phi)
c(cc = weighted_cancor(x, y), abs_phi = abs(phi))
#>        cc   abs_phi 
#> 0.5044296 0.5044296 

# Rows with negative weights and missing values filtered.
set.seed(1)
xw <- factor(sample(letters[1:3], 100, replace = TRUE))
yw <- factor(sample(LETTERS[1:4], 100, replace = TRUE))
w  <- runif(100)
w[c(5, 10)] <- NA_real_
w[20] <- -1  # dropped
xw[30] <- NA  # dropped
yw[40] <- NA  # dropped
weighted_cancor(xw, yw, w = w)
#> [1] 0.2706299

# Larger table where the value equals the dominant CA singular value
x3 <- factor(sample(paste0("x", 1:5), 500, TRUE))
y3 <- factor(sample(paste0("y", 1:6), 500, TRUE))
weighted_cancor(x3, y3)
#> [1] 0.155366

# Degenerate cases
# Single observed level in x after filtering
weighted_cancor(factor(rep("only", 10)), factor(sample(c("a","b"), 10, TRUE)))
#> [1] 0
# No usable data
weighted_cancor(factor(c(NA, NA)), factor(c(NA, NA)), w = c(1, 1))
#> [1] NA