The idea behind the test

Two variables \(Y\) and \(X\) are conditionally independent given \(Z\), written \(Y \perp\!\!\!\perp X \mid Z\), if \(X\) carries no information about \(Y\) once \(Z\) is known. The CCI test (Thorjussen et al. 2024, 2026) turns this into a question about prediction: does \(X\) improve the out-of-sample prediction of \(Y\) when \(Z\) is already used?

If \(Y \perp\!\!\!\perp X \mid Z\), then \(E(Y \mid X, Z) = E(Y \mid Z)\). The same holds if \(X\) is replaced by a random permutation \(X^*\), so a model for \(Y\) should predict equally well with \(X\) as with \(X^*\). With \(M\) a performance metric (e.g. RMSE) of a fitted model \(f_Y\) on new data, the hypotheses are

\[ H_0: Y \perp\!\!\!\perp X \mid Z \;\Rightarrow\; M(f_Y(X^*, Z)) = M(f_Y(X, Z)), \qquad H_1: Y \not\perp\!\!\!\perp X \mid Z \;\Rightarrow\; M(f_Y(X^*, Z)) \neq M(f_Y(X, Z)). \]

The test is computed with Monte Carlo cross-validation:

1. Null distribution. Repeat nperm times: permute \(X\) into \(X^*\), split the data randomly into a training part (share p) and a test part, train a model for \(Y\) from \(X^*\) and \(Z\) on the training part, and record its performance on the test part.
2. Test statistic. Train the same model with the real \(X\) on a random training split, and record its performance on the test part.
3. P-value. Compare the test statistic with the null distribution. If the model with the real \(X\) predicts clearly better than the models with \(X^*\), \(X\) carries information about \(Y\) beyond \(Z\), and \(H_0\) is rejected.

Permuting \(X\) also breaks its relationship with \(Z\). Exact validity of a permutation test would require that the permutation keeps the joint distribution of \((X, Z)\), which is only possible when \(Z\) is categorical (then CCI permutes \(X\) within the groups of \(Z\)). CCI does not try to enforce this. Instead, it relies on the fact that under \(H_0\) both models converge to the same prediction \(E(Y \mid Z)\), and on repeated out-of-sample evaluation. Simulations show good type I error control, and the test is on the conservative side (Thorjussen et al. 2026). The price is computation time, and some loss of power in small samples.

A first test

We simulate data where \(Y\) and \(X\) both depend on \(Z_1\) and \(Z_2\), but not on each other. Then \(Y \perp\!\!\!\perp X \mid Z_1, Z_2\) is true, while \(Y \perp\!\!\!\perp X \mid Z_1\) is false, since \(Z_2\) creates a dependence between \(Y\) and \(X\) when it is left out.

normal_data <- function(n) {
  Z1 <- rnorm(n)
  Z2 <- rnorm(n)
  X <- Z1 + Z2 + rnorm(n)
  Y <- Z1 + Z2 + rnorm(n)
  data.frame(Z1, Z2, X, Y)
}
set.seed(123)
dat <- normal_data(500)

The hypothesis \(Y \perp\!\!\!\perp X \mid Z_1, Z_2\) is written Y ~ X | Z1 + Z2. The only required arguments of CCI.test() are formula and data. Setting seed makes the result reproducible.

res_true <- CCI.test(Y ~ X | Z1 + Z2, data = dat, seed = 1, progress = FALSE)
summary(res_true)
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using rf 
#> Formula:   Y ~ X | Z1 + Z2 
#> Permutations:  160 
#> Metric:    RMSE 
#> Tail:      left 
#> Statistic: 1.143 
#> P-value:   0.7019 
#> 
#> MC sample:   1

The p-value is 0.702, so we do not reject conditional independence. Leaving out \(Z_2\) should give a rejection:

res_false <- CCI.test(Y ~ X | Z1, data = dat, seed = 1, progress = FALSE)
summary(res_false)
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using rf 
#> Formula:   Y ~ X | Z1 
#> Permutations:  160 
#> Metric:    RMSE 
#> Tail:      left 
#> Statistic: 1.492 
#> P-value:   0.006211 
#> 
#> MC sample:   1

The summary shows:

• Method: the learner used, here random forest (rf, the default).
• Formula: the hypothesis tested.
• Permutations: the number of Monte Carlo samples in the null distribution (nperm, default 160).
• Metric: the performance metric, RMSE for a continuous outcome.
• Tail: which side of the null distribution counts as evidence against \(H_0\). For RMSE, lower is better, so a test statistic in the left tail means the model with the real \(X\) predicts better than the models with \(X^*\).
• Statistic and P-value: the performance with the real \(X\), and the empirical p-value. The smallest possible empirical p-value is \(1/(\text{nperm} + 1)\), here 0.0062.
• MC sample: the share of the data used in each Monte Carlo sample (MC_sample) (see Large data sets).

Printing the result gives a shorter overview:

res_false
#> 
#> Computational Conditional Independence Test
#> 
#> Formula:   Y ~ X | Z1 
#> Method:    rf with metric RMSE and 160 permutations
#> Statistic: 1.492   P-value: 0.006211

The null distribution and the test statistic (dashed line) can be plotted:

plot(res_false)

The test statistic lies far to the left of the null distribution: the model predicts \(Y\) much better with the real \(X\) than with a permuted one.

The formula

The formula Y ~ X | Z1 + Z2 means “test \(Y \perp\!\!\!\perp X \mid Z_1, Z_2\)”. The first variable on the right-hand side is \(X\), and everything after | is the conditioning set \(Z\). The following formulas are equivalent ways of writing the same hypothesis:

Y ~ X | Z1 + Z2
Y ~ X + Z1 + Z2   # the first variable after ~ is X, the rest is Z
X ~ Y | Z1 + Z2   # conditional independence is symmetric in X and Y

The test itself is not symmetric: the variable on the left is the one that is predicted. Which direction gives the most power is discussed in vignette("diagnostics-and-tuning", package = "CCI").

To test unconditional independence \(Y \perp\!\!\!\perp X\), write the conditioning set as 1:

set.seed(2)
uncond <- data.frame(X = rnorm(300))
uncond$Y <- sin(2 * uncond$X) + rnorm(300, sd = 0.5)   # dependent, but almost uncorrelated
cor(uncond$X, uncond$Y)
#> [1] 0.2046094
summary(CCI.test(Y ~ X | 1, data = uncond, seed = 1, progress = FALSE))
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using rf 
#> Formula:   Y ~ X | 1 
#> Permutations:  160 
#> Metric:    RMSE 
#> Tail:      left 
#> Statistic: 0.5464 
#> P-value:   0.006211 
#> 
#> MC sample:   1

Y ~ X + 1 gives the same test, while Y ~ X gives an error, to avoid testing unconditional independence by mistake. The example also shows that CCI finds non-linear dependence that a correlation misses.

Choosing the learner

The method argument selects the machine learning model:

Tree-based models (rf, xgboost) are recommended in most cases, as they predict robustly out of sample even when fitted closely to the training data. Model parameters are passed directly to CCI.test():

summary(CCI.test(Y ~ X | Z1, data = dat, method = "xgboost",
                 nrounds = 200, eta = 0.1, max_depth = 3,
                 nperm = 100, seed = 1, progress = FALSE))
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using xgboost 
#> Formula:   Y ~ X | Z1 
#> Permutations:  100 
#> Metric:    RMSE 
#> Tail:      left 
#> Statistic: 1.494 
#> P-value:   0.009901 
#> 
#> MC sample:   1

nrounds is the number of trees for both rf and xgboost (default 600). For rf, mtry sets the number of variables tried in each split, for KNN, k sets the number of neighbours. Other arguments are passed on to the underlying model function (ranger::ranger(), xgboost::xgb.train(), e1071::svm() or kknn::kknn()).

Different types of data

CCI chooses the performance metric from the type of \(Y\) (metric = "Auto"):

• a numeric \(Y\) gives regression, evaluated with RMSE (lower is better);
• a factor, character or logical \(Y\) gives classification, evaluated with Cohen’s Kappa (higher is better).

For classification, metric = "LogLoss" evaluates the predicted class probabilities instead of the predicted classes, and metric can also be set to "Kappa" for a numeric \(Y\) with a few values.

three_classes <- function(n) {
  Z1 <- rnorm(n)
  Z2 <- rnorm(n)
  X <- exp(Z1) + Z2 + rnorm(n, sd = 0.2)
  score <- log(abs(Z1) + 1) + Z2
  Y <- ifelse(score > 0.5, "high", ifelse(score > 0, "medium", "low"))
  data.frame(Z1, Z2, X, Y)
}
set.seed(3)
cat_data <- three_classes(500)
summary(CCI.test(Y ~ X | Z1 + Z2, data = cat_data, seed = 1, progress = FALSE))
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using rf 
#> Formula:   Y ~ X | Z1 + Z2 
#> Permutations:  160 
#> Metric:    Kappa 
#> Tail:      right 
#> Statistic: 0.8978 
#> P-value:   0.4658 
#> 
#> MC sample:   1
summary(CCI.test(Y ~ X | Z1, data = cat_data, metric = "LogLoss", seed = 1, progress = FALSE))
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using rf 
#> Formula:   Y ~ X | Z1 
#> Permutations:  160 
#> Metric:    LogLoss 
#> Tail:      left 
#> Statistic: 0.475 
#> P-value:   0.006211 
#> 
#> MC sample:   1

\(Y\) is here a character vector, which is treated as categorical. For Kappa the tail is right (higher Kappa means better prediction), for RMSE and LogLoss it is left.

Categorical conditioning variables. When \(Z\) contains factor, character or logical variables, and robust = TRUE (default), \(X\) is permuted within the groups defined by the categorical variables in \(Z\). This keeps the relationship between \(X\) and the categorical part of \(Z\), which improves type I error control.

Polynomial and interaction terms. By default, CCI.test() adds polynomial terms up to degree 3 (poly = TRUE, degree = 3) and pairwise interactions (interaction = TRUE) of the numeric conditioning variables to the model. This helps the learner capture smooth relationships between \(Z\) and \(Y\). The terms are named Z1_d_2, Z1_int_Z2 and so on, and are only added to \(Z\), never to \(X\). Set poly = FALSE and interaction = FALSE to use the variables as they are.

Large data sets

The runtime grows with the sample size and the number of Monte Carlo samples. Three arguments help:

• MC_sample = "Auto" (default) uses only a share of the data in each Monte Carlo sample when \(n > 900\), namely \((900/n)^{0.75}\), e.g. about 30% for \(n = 5000\). A different random sample is drawn each time. Use MC_sample = "No" to always use all data, or MC_sample = "Yes" with MC_sample_set to set the share yourself. (Before version 0.3.7 these arguments were called subsample and subsample_set; the old names still work, with a warning.)
• method = "KNN" is by far the fastest learner.
• nthread sets the number of threads for rf and xgboost.

hard_case <- function(n) {
  Z1 <- runif(n, -2, 2)
  Z2 <- runif(n, -2, 2)
  hZ <- sin(Z1) * cos(Z2)
  X <- hZ + 0.2 * rnorm(n)
  Y <- hZ^2 + 0.2 * rnorm(n)
  data.frame(X, Y, Z1, Z2)
}
set.seed(4)
big <- hard_case(5000)
start <- Sys.time()
res_big <- CCI.test(Y ~ X | Z1 + Z2, data = big, method = "KNN", seed = 1, progress = FALSE)
Sys.time() - start
#> Time difference of 2.521336 secs
summary(res_big)
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using KNN 
#> Formula:   Y ~ X | Z1 + Z2 
#> Permutations:  160 
#> Metric:    RMSE 
#> Tail:      left 
#> Statistic: 0.2176 
#> P-value:   0.2671 
#> 
#> MC sample:   0.28

The “MC sample” line of the summary shows that a share of 0.28 of the data was used in each Monte Carlo sample, i.e. \((900/5000)^{0.75}\). For guidance on runtimes with the other learners, see Thorjussen et al. (2026).