---
title: "User guide to the citcdf R package"
author: "Boris Hejblum"
date: today
format:
  html:
    toc: true                      
    number-sections: true         
    html-math-method: mathjax
bibliography: citcdf.bib
link-citations: true
vignette: >
  %\VignetteIndexEntry{User guide to the citcdf R package}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r knitrsetup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  dpi = 96
)
# Datasets borrowed from other packages (Suggests, not hard dependencies).
has_bnlearn <- requireNamespace("bnlearn", quietly = TRUE)
has_seurat <- requireNamespace("SeuratObject", quietly = TRUE) &&
  utils::packageVersion("SeuratObject") >= "5.0.0"
has_sessioninfo <- requireNamespace("sessioninfo", quietly = TRUE)
has_reactable <- requireNamespace("reactable", quietly = TRUE)
```

![](../man/figures/logo.svg){width=139px}

# Overview of `citcdf`

`citcdf` performs **c**onditional **i**ndependence **t**esting through
conditional **c**umulative **d**istribution **f**unctions
[@gauthier2021]. It addresses the following null hypothesis: 
$$H_0: Y \perp\!\!\!\perp X \mid Z$$ 

by testing whether $F_{Y \mid X, Z}(y) = F_{Y \mid Z}(y)$. The associated test statistic is computed across a grid of thresholds $\omega_1 < \dots < \omega_p$ spanning $Y$, where the indicator $\mathbb{1}_{Y_i \le \omega_j}$ is
regressed on both $X$ and $Z$ at each threshold $j$. Under $H_0$, the coefficients $\beta_j$ carried by $X$ are all null, and thus `citcdf` test statistic is the sum of squares whose asymptotic distribution is then a weighted mixture of $\chi^2_1$.

No distributional assumption is made on $Y$, and in that sense `citcdf` is *distribution-free*. It is therefore robust to zero-inflation, multi-modality and skewness that typical occur in single-cell RNA-seq data.


## Main user functions

Three main functions build form the `citcdf`package leverage this test statistic:

 - `cit_asymp()`: one hypothesis, asymptotic p-value
 - `cit_perm()`: one hypothesis, permutation p-value
 - `cit_multi()`: many outcomes at once (with a Benjamini-Hochberg adjustment), either asymptotic or permutation test

## Inputs

`Y` is a numeric vector. `X` and `Z` are data frames, one column per variable,
numeric or factor. Gene set analysis is available through `cit_gsa()` and is not
covered here.

```{r load}
library(citcdf)
set.seed(20260817)
```

# Testing a single hypothesis

## The data

We use the `marks`dataset, from the `bnlearn` package [@scutari2010]. It records the exam marks of 88 students in five mathematics topics [@mardia1979]. Of note, *mechanics* and *statistics* marks are **correlated**, but the association is mediated by the *algebra* mark (relevant to both disciplines).

```{r marks-data, eval = has_bnlearn}
data("marks", package = "bnlearn")
str(marks)
```

## Asymptotic test with `cit_asymp()`

$Y$ is the statistics mark, $X$ the mechanics mark. Without conditioning:

```{r asymp-uncond, eval = has_bnlearn}
Y <- marks$STAT
X <- data.frame(MECH = marks$MECH)

cit_asymp(Y, X)
```

Conditioning on the algebra mark asks whether mechanics adds anything beyond
algebra:

```{r asymp-cond, eval = has_bnlearn}
Z <- data.frame(ALG = marks$ALG)

cit_asymp(Y, X, Z)
```

The evidence disappears. 

**NB:** `cit_asymp()` assumed neither normality nor linearity.  
It remains valid for skewed, zero-inflated or multi-modal outcomes.

Of note, `bnlearn::ci.test()` reaches the same conclusion from a
Gaussian correlation statistic:
```{r bnlearn-check, eval = has_bnlearn}
bnlearn::ci.test("STAT", "MECH", data = marks, test = "cor")$p.value
bnlearn::ci.test("STAT", "MECH", "ALG", data = marks, test = "cor")$p.value
```

## Permutation test with `cit_perm()`

The asymptotic null distribution requires a reasonable sample size. For small
`n`, `cit_perm()` calibrates the same observed statistic against a permutation
null:

```{r perm-uncond, eval = has_bnlearn}
cit_perm(Y, X, n_perm = 1000)
```

`test_statistic` remains identical to the one returned by `cit_asymp()`, and `score` counts the permutations reaching the observed statistic: the (unbiased) permutation p-value is then computed as `(score + 1) / (n_perm + 1)` (so `n_perm` also fixes the smallest reachable p-value, ie. `1/(n_perm +1)`.

Permuting $X$ freely would destroy its association with $Z$ and test the wrong
null. `X_perm()` computes design permutations conditionally on $Z$:

 - $Z$ discrete: permutes within each strata of $Z$
 - $Z$ continuous: permutes by matching observations on the fitted value of $X$ given $Z$, following @berrett2020 (see `?perm_cont`)

```{r perm-cond, eval = has_bnlearn}
cit_perm(Y, X, Z = Z, n_perm = 1000)
```

Both tests agree on both hypotheses.

**NB:** `cit_perm()` only accepts a single covariate column when `Z` is not `NULL`. Multiple `X` with `Z` present require the asymptotic test.

## Visualization

`plot_compare_ccdf()` displays what the statistic compares: the CCDF of $Y$
given $X$ against its marginal counterpart. The further apart the 2 are, the larger $D$:

```{r ccdf-plot, eval = has_bnlearn, fig.height = 3.2}
plot_compare_ccdf(Y = data.frame(STAT = Y), X = X)
```

# Testing many outcomes with `cit_multi()`

`cit_multi()` loops the test over the columns of a matrix `M` of `n`
observations by `r` outcomes, and returns Benjamini-Hochberg adjusted p-values (alongside the raw ones).

## scRNA-seq data

This vignette uses `pbmc_small`, the PBMC scRNA-seq excerpt distributed with the `SeuratObject` package [@satija2023] (230 genes, 80 cells).

**NB:** `M` must be oriented *cells-by-genes*, the transpose of usual
genes-by-cells expression matrix. Also, **outcomes with no variability must be filtered-out beforehand** (otherwise an error is triggered).

```{r pbmc-data, eval = has_seurat}
data("pbmc_small", package = "SeuratObject")
expr <- as.matrix(SeuratObject::LayerData(pbmc_small, layer = "data"))
cell_md <- pbmc_small[[]]

expr <- expr[apply(expr, 1, function(g) length(unique(g)) > 1), ]
M <- as.data.frame(t(expr))
dim(M)
```

The variable of interest is the cell population (`letter.idents`, two populations). The two populations were not sequenced to comparable depth:

```{r pbmc-confounding, eval = has_seurat}
tapply(cell_md$nCount_RNA, cell_md$letter.idents, median)
```

Library size is therefore a candidate confounder. An unadjusted comparison will
identify significant genes that are only associated with sequencing depth.

## Asymptotic analysis

```{r pbmc-asymp, eval = has_seurat}
X_pop <- data.frame(pop = cell_md$letter.idents)
Z_lib <- data.frame(libsize = cell_md$nCount_RNA)

res_unadj <- cit_multi(M, X = X_pop, test = "asymptotic", parallel = FALSE)
res_adj <- cit_multi(M, X = X_pop, Z = Z_lib, test = "asymptotic",
                     parallel = FALSE)

c(unadjusted = sum(res_unadj$pvals$adj_pval < 0.05),
  adjusted   = sum(res_adj$pvals$adj_pval < 0.05))
```

Conditioning on library size withdraws a third of the hits. The strongest
signals — canonical monocyte markers — survive:

```{r pbmc-top, eval = has_seurat}
head(res_adj$pvals[order(res_adj$pvals$raw_pval), ], 5)
```

`res_adj` is a list containing:

 - `which_test` and `n_perm`, the test performed and the number of permutations
 (`NA` for the asymptotic test)
 - `pvals`, the gene-wise raw and BH-adjusted p-values, plus the test statistics

The `plot()` method sorts the raw p-values against the BH threshold and the
nominal level:

```{r pbmc-plot, eval = has_seurat}
plot(res_adj)
```

## Visualization

`plot_compare_ccdf()` details one gene. Panel A contrasts the CCDF given the
cell population with the marginal CDF. Panel B repeats the comparison within
library-size quartiles.

```{r pbmc-ccdf, eval = has_seurat, fig.height = 5}
top_gene <- rownames(res_adj$pvals)[which.min(res_adj$pvals$raw_pval)]
plot_compare_ccdf(Y = M[, top_gene, drop = FALSE], X = X_pop, Z = Z_lib,
                  space_y = TRUE, number_y = 20)
```

## Permutation analysis

With 80 cells, the asymptotic approximation should work (low end of its range).
`test = "permutation"` confirms a few genes. `adaptive = TRUE` (the default)
spends additional computation time to increase p-value resolution where it matters (especially for FDR-adjusted p-values): all outcomes start at `n_perm`, and only
those still significant proceed to the larger stages of `n_perm_adaptive`.

```{r pbmc-perm, eval = has_seurat}
res_perm <- cit_multi(M, X = X_pop, Z = Z_lib, test = "permutation",
                      n_perm = 100, parallel = FALSE)
res_perm$n_perm
```

```{r pbmc-perm-compare, eval = has_seurat && has_reactable, echo = has_seurat && has_reactable}
order_asymp <- rownames(res_adj$pvals[order(res_adj$pvals$adj_pval),])
reactable::reactable(data.frame(asymptotic = signif(res_adj$pvals[order_asymp, "raw_pval"], 3),
           permutation = signif(res_perm$pvals[order_asymp, "raw_pval"], 3),
           row.names = order_asymp),
           defaultPageSize=15
)
```

```{r pbmc-perm-compare-noreactable, eval = has_seurat && !has_reactable, echo = has_seurat && !has_reactable}
order_asymp <- rownames(res_adj$pvals[order(res_adj$pvals$adj_pval),])
knitr::kable(data.frame(asymptotic = signif(res_adj$pvals[order_asymp, "raw_pval"], 3),
           permutation = signif(res_perm$pvals[order_asymp, "raw_pval"], 3),
           row.names = order_asymp)[1:15,])
```

The orderings nearly coincide. Individual p-values, as expected, differ: the permutation ones carry Monte Carlo noise and cannot fall below
`1 / (total permutations + 1)`.

# Practical considerations

 - [ ] **Threshold grid:** `space_y = FALSE` places a threshold at every distinct value of `Y`, while `space_y = TRUE` uses a regular grid of `number_y` points instead, at a cost independent of `n`.  
Defaults differ: `cit_asymp()` and `cit_perm()` use the exhaustive grid, `cit_multi()` uses `number_y = 10`. Expect small numerical differences between a single call and a `cit_multi()` run.

 - [ ] **Parallelism.** `cit_multi()` parallelises over outcomes and defaults to
`parallel = interactive()`. Calls above set `parallel = FALSE` for reproducible
builds. Drop it in real analyses, and set `n_cpus`.

 - [ ] **Multiplicity.** `cit_multi()` returns Benjamini-Hochberg adjusted p-values in `adj_pval` for FDR control. `cit_asymp()` and `cit_perm()` are single-hypothesis functions and do not need such multiple-testing adjustment.

 - [ ] **Design.** The asymptotic test accepts several variables of interest (`ncol(X) > 1`) and several covariates (`ncol(Z) > 1`). The permutation test is less flexible accepting at most one covariate.

 - [ ] **Preprocessing.** Normalization remains the user's choice and responsibility.

# Session information

```{r session, eval = has_sessioninfo, echo=FALSE, collapse = FALSE, comment = ""}
sessioninfo::session_info()
```

```{r session-base, eval = !has_sessioninfo, echo = FALSE}
sessionInfo()
```

# References

<div id="refs"></div>
