---
title: "A benchmark harness and recommender for causal discovery"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{A benchmark harness and recommender for causal discovery}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- Render time: ~2.5 s under rmarkdown::render(); macOS arm64, R 4.5.2,
     one core; pcalg, bnlearn, ggplot2 and proxymix installed, no Python
     discovery back-end installed on this build. -->

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  dpi = 150
)
```

```{r library}
library(cdzoo)
has_pcalg <- requireNamespace("pcalg", quietly = TRUE)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
has_proxymix <- requireNamespace("proxymix", quietly = TRUE)
set.seed(20260902)
```

## Why

An analyst sitting on a data set of several observational variables, who
needs a defensible adjustment set for one variable's effect on another, asks
a blunt question: *I do not trust running one causal-discovery algorithm and
taking its graph as the answer -- which method should I run on data that
looks like mine, how sure should I be, and what would the alternatives have
said?* `cdzoo` is built to answer exactly that, for a problem with no
run-time ground truth to check any single algorithm's answer against: it
reads the data's own signature, routes to the method a controlled benchmark
shows is suited to that signature, runs it, and hands back the graph
together with its reasoning, a confidence label, and the alternatives it did
not choose -- or, when even that is not certainty enough, the
stability-weighted agreement of several methods run side by side.

## What

`list_methods()` lists every registered discovery method, whether it needs
an R package (`pcalg`, `bnlearn`) or a Python module reached through
`reticulate`, and whether that back-end is installed on this machine.
`cd_discover()` runs one named method directly on a data matrix and returns
a `cd_result`: the recovered graph as an adjacency matrix in the `pcalg`
`amat.cpdag` convention, where `graph[i, j] == 1` and `graph[j, i] == 0`
encodes a directed edge from `i` to `j`, and
`graph[i, j] == graph[j, i] == 1` encodes an edge whose direction was not
resolved. `cd_edges()` turns that matrix into a tidy edge list.

Because discovery has no run-time ground truth, `cd_diagnose()` first
describes the data's own characteristics -- its non-Gaussianity, a
curvature heuristic for nonlinear mechanisms, the condition number of its
correlation matrix, the spread of its per-node residual variances -- and
`cd_auto()` routes on those characteristics to the method the benchmark
evidence
(`cd_field_benchmark`) shows is suited to that regime, runs it, and reports
why in its `rationale`, how sure in its `confidence`, and what it did not
choose in its `alternatives`. `cd_explain()` prints the fuller, auditable
account.
`mode = "consensus"` runs several benchmark-suited methods at once and
returns their stability-weighted union, with the share of methods that
agreed on each edge. `cd_score()` checks a recovered graph against a known
truth on the built-in scenario grid (`cd_scenario()`), and `cd_benchmark()`
runs that grid for a chosen set of methods and scenarios.

## Do

### The registry

```{r registry}
reg <- list_methods()
```

```{r registry-kable, echo = FALSE}
knitr::kable(
  reg[, c("id", "name", "family", "language", "installed")],
  caption = paste0(
    "Every registered discovery method (n = ", nrow(reg), "), and whether ",
    "its back-end is installed on the machine that rendered this ",
    "vignette (", sum(reg$installed), " are)."
  )
)
```

Of the `r nrow(reg)` registered methods, `r sum(reg$installed)` have an
installed back-end here. Nothing below needs every row of this table to be
`TRUE`: a method whose back-end is missing is simply not a candidate the
recommender can pick, and it says so when that changes its answer.

### Running one method directly

```{r discover, eval = has_pcalg}
sim <- cd_scenario("linear_gaussian", seed = 1L)
res <- cd_discover(sim$data, "pc_pcalg")
res
head(cd_edges(res))
```

```{r discover-skip, eval = !has_pcalg, echo = FALSE, results = "asis"}
cat("pcalg is not installed on this build, so the direct PC run above is",
    "skipped.\n")
```

`cd_discover()` is the plain dispatcher: name a method, get its graph, with
no routing and no opinion attached. It is what `cd_auto()` calls once it has
decided which method that should be.

### The recommender explains itself

A nonlinear additive-noise chain (`nonlinear_anm`) is the regime the
registry's kernel conditional-independence method was built for.
`cd_diagnose()` shows what the recommender reads before it decides
anything:

```{r nonlinear-diagnose}
nl <- cd_scenario("nonlinear_anm", seed = 1L)
knitr::kable(
  cd_diagnose(nl$data),
  caption = "The transparent characteristics cd_auto() reads for this data."
)
```

Asking for the recommendation alone, with no method actually run, costs
nothing and needs no back-end at all:

```{r nonlinear-recommend}
rec_only <- cd_auto(nl$data, report = FALSE)
cd_explain(rec_only)
```

Running the recommendation produces the graph on top of the same rationale
and confidence:

```{r nonlinear-run}
rec <- cd_auto(nl$data)
res_nl <- cd_edges(rec)
res_nl
```

`r rec@method_name` is not the route's first choice; the rationale above
says which method was preferred and why it was not used. This is the
package's identity in one call: it never decrees a method, it names the one
it ran, the evidence that pointed at it, how confident that evidence leaves
it, and the alternatives a differently-equipped machine would have reached
for instead.

### When unsure, ask several: consensus

```{r consensus, eval = has_pcalg}
lap <- cd_scenario("linear_laplace", seed = 1L)
con <- cd_auto(lap$data, mode = "consensus", top_k = 3L)
con
edges_con <- cd_edges(con, agreement = TRUE)
```

```{r consensus-kable, echo = FALSE, eval = has_pcalg}
# pandoc's smart-dash rendering turns the raw "-->"/"---" edge codes into
# unicode dashes inside an HTML table; relabel for display only.
edges_con_display <- edges_con
edges_con_display$type <- ifelse(
  edges_con_display$type == "-->", "directed", "undirected"
)
knitr::kable(
  edges_con_display, digits = 2,
  caption = paste0(
    "Edges in the consensus graph, with the share of the ",
    length(con@provenance$runs), " combined methods that recovered each ",
    "one's adjacency (orientation may still differ between them)."
  )
)
```

```{r consensus-skip, eval = !has_pcalg, echo = FALSE, results = "asis"}
cat("pcalg is not installed on this build, so the consensus example above",
    "is skipped.\n")
```

The full pairwise picture -- including the variable pairs no method
proposed an edge for -- is a heat map, which a table of found edges alone
cannot show:

```{r fig-consensus-heatmap, eval = has_pcalg && has_ggplot2, fig.height = 4.5, fig.cap = "Pairwise edge agreement among the consensus-mode methods: the share of the combined methods that recovered an edge (in either direction) between each pair of variables, ignoring orientation. Darker tiles mark pairs every method agreed carry an edge."}
agr <- con@provenance$agreement
idx <- which(upper.tri(agr), arr.ind = TRUE)
heat_df <- data.frame(
  from = rownames(agr)[idx[, 1L]],
  to = colnames(agr)[idx[, 2L]],
  agreement = agr[idx],
  stringsAsFactors = FALSE
)
ggplot2::ggplot(
  heat_df, ggplot2::aes(x = from, y = to, fill = agreement)
) +
  ggplot2::geom_tile(colour = "white") +
  ggplot2::scale_fill_viridis_c(
    name = "share agreeing", limits = c(0, 1), option = "D"
  ) +
  ggplot2::labs(
    x = NULL, y = NULL,
    title = "Consensus edge agreement, linear_laplace scenario"
  ) +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::coord_fixed()
```

```{r fig-consensus-heatmap-skip, eval = !(has_pcalg && has_ggplot2), echo = FALSE, results = "asis"}
cat("The consensus heat map needs pcalg and ggplot2 both installed; it is",
    "skipped when either is absent.\n")
```

```{r consensus-text, eval = has_pcalg, echo = FALSE, results = "asis"}
n_methods <- length(con@provenance$runs)
n_edges <- nrow(edges_con)
n_full <- sum(edges_con$agreement >= 0.999)
n_majority <- n_edges - n_full
cat(sprintf(
  paste0(
    "Of the %d edges that cleared the consensus's majority-agreement ",
    "threshold out of the %d combined methods, %d were found by every ",
    "one of the %d methods and %d by a bare majority alone; the heat ",
    "map's palest tiles are the pairs no majority agreed on at all, ",
    "which the edge table above does not show because it only lists ",
    "edges the consensus actually carried.\n"
  ),
  n_edges, n_methods, n_full, n_methods, n_majority
))
```

```{r consensus-text-skip, eval = !has_pcalg, echo = FALSE, results = "asis"}
cat("pcalg is not installed on this build, so the consensus agreement",
    "reading above is skipped.\n")
```

### The benchmark evidence behind the routing

`cd_field_benchmark` is shipped with the package, so this section needs no
back-end. It is a summary of `r length(unique(cd_field_benchmark$method))`
established methods run across the
`r nrow(cd_scenarios())` structural-causal-model scenarios of
`cd_scenarios()`, each scored against its known DAG and averaged over
random-graph seeds -- the evidence `cd_auto()` routes on.

```{r ranking}
agg <- aggregate(
  cbind(shd, skeleton_f1, runtime_s) ~ method_name + family,
  data = cd_field_benchmark, FUN = mean
)
agg <- agg[order(agg$shd), ]
```

```{r ranking-kable, echo = FALSE}
knitr::kable(
  agg, row.names = FALSE, digits = 3L,
  caption = paste0(
    "Mean structural Hamming distance (lower is better), skeleton F1, and ",
    "runtime for each method, averaged over the ",
    nrow(cd_scenarios()), " scenarios."
  )
)
```

```{r fig-accuracy, eval = has_ggplot2, fig.height = 4.5, fig.cap = "Mean structural Hamming distance (lower is better) for each registered method, averaged over the nine scenarios and coloured by algorithm family."}
ord <- agg[order(-agg$shd), "method_name"]
plot_df <- agg
plot_df$method_name <- factor(plot_df$method_name, levels = ord)
ggplot2::ggplot(plot_df, ggplot2::aes(method_name, shd, fill = family)) +
  ggplot2::geom_col() +
  ggplot2::coord_flip() +
  ggplot2::scale_fill_viridis_d(option = "D", end = 0.9) +
  ggplot2::labs(
    x = NULL, y = "mean structural Hamming distance (lower is better)",
    fill = "family", title = "Discovery accuracy across nine scenarios"
  ) +
  ggplot2::theme_minimal(base_size = 11)
```

```{r fig-accuracy-skip, eval = !has_ggplot2, echo = FALSE, results = "asis"}
cat("ggplot2 is not installed on this build, so the accuracy figure above",
    "is skipped.\n")
```

```{r fig-speed, eval = has_ggplot2, fig.height = 4.5, fig.cap = "The speed-accuracy trade-off: mean runtime (log scale) against mean structural Hamming distance for each method, coloured by algorithm family."}
ggplot2::ggplot(
  agg, ggplot2::aes(runtime_s, shd, colour = family, label = method_name)
) +
  ggplot2::geom_point(size = 2.4) +
  ggplot2::scale_x_log10() +
  ggplot2::scale_colour_viridis_d(option = "D", end = 0.9) +
  ggplot2::labs(
    x = "mean runtime (seconds, log scale)",
    y = "mean structural Hamming distance", colour = "family",
    title = "The speed-accuracy trade-off"
  ) +
  ggplot2::theme_minimal(base_size = 11)
```

```{r fig-speed-skip, eval = !has_ggplot2, echo = FALSE, results = "asis"}
cat("ggplot2 is not installed on this build, so the speed-accuracy figure",
    "above is skipped.\n")
```

```{r speed-text, echo = FALSE, results = "asis"}
kci_rt <- agg$runtime_s[agg$method_name == "PC + KCI (kernel CI)"]
lingam_rt <- min(agg$runtime_s[agg$family == "functional"])
ratio <- round(kci_rt / lingam_rt)
cat(sprintf(paste0(
  "The continuous-optimisation and LiNGAM families cluster at the ",
  "low-error end of the accuracy figure; the kernel- and constraint-based ",
  "search methods sit at the high-error end on these mostly-linear ",
  "scenarios, which is why the recommender does not route to them by ",
  "default outside a nonlinear-mechanism signature. The kernel ",
  "conditional-independence method runs roughly %d times slower than the ",
  "fastest LiNGAM entrant while landing at a worse mean structural ",
  "Hamming distance, which is why the route above confined it to the ",
  "nonlinear signature rather than using it as a general-purpose ",
  "default.\n"
), ratio))
```

### Scoring and benchmarking your own runs

```{r score, eval = has_pcalg}
sim_b <- cd_scenario("linear_bimodal", seed = 2L)
score_one <- cd_score(cd_discover(sim_b$data, "icalingam_pcalg"), sim_b$truth)
```

```{r score-kable, echo = FALSE, eval = has_pcalg}
knitr::kable(
  score_one, digits = 3,
  caption = "cd_score() of ICA-LiNGAM against the known truth graph."
)
```

`cd_benchmark()` runs a grid of methods over a grid of scenarios and seeds,
isolating each cell so one back-end's failure records a status rather than
stopping the sweep:

```{r benchmark, eval = has_pcalg}
bench <- cd_benchmark(
  methods = c("pc_pcalg", "ges_pcalg"),
  scenarios = c("linear_gaussian", "linear_bimodal"),
  seeds = 1:2, verbose = FALSE
)
```

```{r benchmark-kable, echo = FALSE, eval = has_pcalg}
knitr::kable(
  bench, digits = 3, row.names = FALSE,
  caption = "cd_benchmark() over two methods, two scenarios, two seeds."
)
```

```{r score-benchmark-skip, eval = !has_pcalg, echo = FALSE, results = "asis"}
cat("pcalg is not installed on this build, so the scoring and benchmark",
    "examples above are skipped.\n")
```

### Beyond a sample: an evaluable density

The sampling-based methods need a sample. When a density can be evaluated
but not sampled -- given as a mixture rather than draws -- the
conditional-independence (Markov) structure can be read off it in closed
form instead. Passing a `proxymix` mixture to `cd_auto()` takes this route,
which no sampling-based method in the registry can serve:

```{r density, eval = has_proxymix}
# a three-variable mixture whose precision makes x1 and x3 conditionally
# independent given x2 -- a chain x1 -- x2 -- x3
sigma_g <- solve(matrix(
  c(1, -0.5, 0, -0.5, 1.25, -0.5, 0, -0.5, 1), 3L, 3L
))
g <- proxymix::gmm(
  weights = 1, means = list(c(0, 0, 0)), covariances = list(sigma_g)
)
res_density <- cd_auto(g)
res_density@rationale
cd_edges(res_density)
```

```{r density-skip, eval = !has_proxymix, echo = FALSE, results = "asis"}
cat("proxymix is not installed on this build, so the evaluable-density",
    "example above is skipped.\n")
```

## Read

`r rec@method_name` answers the opening question in miniature. On the
nonlinear scenario, `cd_auto()` did not decree an answer: its rationale
named the method the benchmark prefers once a curvature heuristic clears
`r rec_only@provenance$thresholds$nonlinearity`, said plainly that method is
not installed here, named the substitute it ran instead, and downgraded its
own confidence from `moderate` to `r rec@confidence` because the substitute
is not the route's preferred method, exactly as the rationale states. An
analyst reading only the confidence label, without the rationale attached to
it, would not know whether `low` means the data itself is ambiguous or the
right tool for it simply was not on the machine; `cd_auto()` never lets that
distinction go unstated.

Consensus mode answers the same trust question a second way: instead of one
run's self-reported confidence, it shows how much several benchmark-suited
methods agree once they are actually run side by side. The heat map above
carries information the edge table cannot -- every pair the combined
methods considered and rejected, not only the ones they proposed an edge
for -- which is exactly the fuller picture an analyst who does not trust a
single algorithm's graph is asking for.

The benchmark figures ground the routing itself: the continuous-optimisation
and LiNGAM families cluster at the accurate, cheap end of the trade-off,
and the kernel conditional-independence method is markedly slower and less
accurate on these mostly-linear scenarios, which is exactly why `cd_auto()`
reserves it for a signature that specifically flags nonlinearity rather than
using it as a default.

## Limits

This vignette never checks a recommendation against a ground truth outside
the built-in scenario grid. On genuinely unlabelled data, `cd_score()` has
nothing to score against, and a confidence label describes how well the
data's signature matches a regime the benchmark was run on, not a promise
about this particular data set. Consensus mode's agreement is stability
across methods run on the same data, not independent corroboration from a
second data source: methods that share an assumption can agree confidently
and still be wrong together. `cd_field_benchmark`'s nine scenarios are
simulated and, however varied their noise families and graph shapes, cannot
cover every mechanism a real data set might hide. This build has no Python
discovery back-end installed, so every route that would have preferred
DirectLiNGAM, RCD, NOTEARS, GOLEM, DAGMA or PC+KCI ran its next-best R
alternative instead and said so in its rationale, exactly as shown above; a
build with those Python modules installed through `reticulate` would run
each route's first choice and report higher confidence on the regimes those
methods target. Installing a different set of back-ends, adding scenarios
covering mechanisms outside this grid, or asking for a different `top_k` in
consensus mode would each change which method is preferred, which is
substituted, and which edges the consensus carries.

## What to read next

*Discovery from sufficient statistics* covers the case this vignette does
not: no row-level data at all, only an aggregated correlation matrix or
pairwise co-occurrence counts, and what the recommender can and cannot say
about a data set it can only see through its second moments.

## Reproduce

`set.seed(20260902)` guards the top-level session; every simulated data set
in this vignette is additionally seeded explicitly through `cd_scenario()`'s
own `seed` argument (`1L` or `2L` as shown in each call), so the same call
reproduces the same data regardless of what ran before it in the session.

```{r session-info}
sessionInfo()
```