--- title: "Getting started with CCI" output: rmarkdown::html_vignette bibliography: references.bib vignette: > %\VignetteIndexEntry{Getting started with CCI} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4 ) ``` ```{r setup} library(CCI) ``` This vignette shows how to test conditional independence with the **CCI** package: the idea behind the test, how to run it, how to read the output, and how it adapts to different kinds of data. The other vignettes cover: - `vignette("diagnostics-and-tuning", package = "CCI")`: checking that a result is reliable, tuning the learner and choosing the test direction. - `vignette("custom-models-and-metrics", package = "CCI")`: using your own machine learning model or performance metric. - `vignette("applied-examples", package = "CCI")`: testing a causal model (DAG) and time series. # 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 [@Thorjussen2024; @Thorjussen2026] 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 [@Thorjussen2026]. 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. ```{r} 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. ```{r} res_true <- CCI.test(Y ~ X | Z1 + Z2, data = dat, seed = 1, progress = FALSE) summary(res_true) ``` The p-value is `r round(res_true$p.value, 3)`, so we do not reject conditional independence. Leaving out $Z_2$ should give a rejection: ```{r} res_false <- CCI.test(Y ~ X | Z1, data = dat, seed = 1, progress = FALSE) summary(res_false) ``` 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 `r round(1 / 161, 4)`. - **MC sample**: the share of the data used in each Monte Carlo sample (`MC_sample`) (see [Large data sets](#large-data-sets)). Printing the result gives a shorter overview: ```{r} res_false ``` The null distribution and the test statistic (dashed line) can be plotted: ```{r} 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: ```{r, eval = FALSE} 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`: ```{r} 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) summary(CCI.test(Y ~ X | 1, data = uncond, seed = 1, progress = FALSE)) ``` `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: | `method` | Model | Comment | |---|---|---| | `"rf"` (default) | Random forest (ranger) | Robust, a good first choice | | `"xgboost"` | Gradient boosting (xgboost) | Often more power, slower, more tuning parameters | | `"svm"` | Support vector machine (e1071) | Smooth relationships | | `"KNN"` | k-nearest neighbours (kknn) | Very fast, useful for large data | 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()`: ```{r} summary(CCI.test(Y ~ X | Z1, data = dat, method = "xgboost", nrounds = 200, eta = 0.1, max_depth = 3, nperm = 100, seed = 1, progress = FALSE)) ``` `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. ```{r} 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)) summary(CCI.test(Y ~ X | Z1, data = cat_data, metric = "LogLoss", seed = 1, progress = FALSE)) ``` $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`. ```{r} 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 summary(res_big) ``` The "MC sample" line of the summary shows that a share of `r round(res_big$MC_sample, 2)` 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 @Thorjussen2026. # References