---
title: "Causal Discovery with lingamr"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Causal Discovery with lingamr}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 5
)
set.seed(0)
```

This vignette walks through a complete causal discovery workflow with
`lingamr`, step by step, using sample data.

```{r setup}
library(lingamr)
```


## Sample Data

`lingamr` provides five sample data generators. Each returns a list
containing `data` (a data frame) and `true_adjacency` (the true adjacency matrix).

| Function | Variables | Default n | Characteristics |
|------|:------:|:-----------:|------|
| `generate_lingam_sample_6()` | 6 | 1,000 | Standard fixed structure. The main example in this vignette |
| `generate_lingam_sample_10()` | 10 | 1,000 | An extension of the 6-variable case (used in [A Larger Dataset](#a-larger-dataset-10-variables)) |
| `generate_lingam_hard_sample()` | 9 | 200 | A difficult setting with strong multicollinearity |
| `generate_lingam_large_sample()` | variable | 1,000 | A random sparse DAG with an arbitrary number of variables (used in [The Scalability Wall](#when-there-are-many-variables-the-scalability-wall)) |
| `generate_lingam_paradox_data()` | 4 | 2,000 | The measurement error paradox (used in [The Paradox Example](#a-case-where-directlingam-struggles-the-measurement-error-paradox)) |

### generate_lingam_sample_6()

`generate_lingam_sample_6()` returns artificial data following a 6-variable
LiNGAM model, together with its true adjacency matrix.
The data is stored in `data` and the adjacency matrix in `true_adjacency`.

```{r sample1_matrix6_1}
x1k <- generate_lingam_sample_6(n = 1000)

x1k$data |>
  head()
```

```{r sample1_matrix6_2}
x1k$true_adjacency
```

`plot_adjacency()` draws a causal graph based on the adjacency matrix.

```{r sample1_matrix6_3}
x1k$true_adjacency |>
  plot_adjacency(
    labels  = colnames(x1k$data),
    title   = "True causal structure",
    rankdir = "TB",
    shape   = "circle"
  )
```


## Causal Discovery

`lingam_direct()` runs Direct LiNGAM. By default, independence is assessed
using mutual information, and path coefficients are computed with adaptive
LASSO regression.

```{r lingam_direct}
model <- x1k$data |>
  lingam_direct()
```

To use HSIC for assessing independence, set the `measure` argument to
"kernel". HSIC is computationally expensive; for `n > 1000`,
`lingam_direct()` automatically switches to a low-rank approximation.


### Causal Order

The estimated causal order is stored in `causal_order` as index numbers.
```{r causal_order}
# index number
model$causal_order

# variable name
colnames(x1k$data)[model$causal_order]
```


### Estimated Adjacency Matrix

We inspect the estimated effect magnitudes. By default, the regression
coefficients from adaptive LASSO regression are used.

```{r adjacency_matrix}
model$adjacency_matrix |>
  round(3)
```


### Drawing the Causal Graph

We draw the causal graph based on the adjacency matrix estimated by Direct LiNGAM.

```{r plot_adjacency}
model$adjacency_matrix |>
  plot_adjacency(
    labels    = colnames(model$adjacency_matrix),
    title     = "Estimated Causal Structure (Direct LiNGAM)",
    rankdir   = "TB",
    shape     = "ellipse",
    fillcolor = "lightgreen"
  )
```


### Comparing the Estimated and True Structures

When the true structure is known, as with sample data, you can pass the true
adjacency matrix to the `true_B` argument of `plot_adjacency()` to color-code
the estimated edges by comparing them against the true structure. This lets you
assess estimation accuracy at a glance, which is useful for validating methods
or for educational purposes.

- **Green (solid)**: correctly detected edges (estimated and true)
- **Red (solid)**: falsely detected edges (estimated but not true)
- **Orange (dashed)**: missed edges (true but not estimated; the true coefficient is shown)

```{r plot_adjacency_compare}
model$adjacency_matrix |>
  plot_adjacency(
    labels  = colnames(model$adjacency_matrix),
    true_B  = x1k$true_adjacency,
    title   = "Estimated vs. True Structure",
    rankdir = "TB",
    shape   = "ellipse"
  )
```


### Static Plotting with ggplot2

While `plot_adjacency()` returns an interactive HTML figure via DiagrammeR,
`autoplot()` draws the same causal structure as a static, ggplot2-based figure.
This is stable for image and PDF output in R Markdown / Quarto, and you can
layer ggplot2 functions on top to set themes or titles afterward. Node positions
are computed using `igraph`'s hierarchical layout, so the causal flow generally
runs from top to bottom.

`autoplot()` is a ggplot2 generic, so call it as `ggplot2::autoplot()` or load
it beforehand with `library(ggplot2)` (plotting requires `ggplot2` and `igraph`).

```{r autoplot, eval = requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("igraph", quietly = TRUE)}
ggplot2::autoplot(model)
```


## Total Causal Effect

The **total causal effect** is the overall impact of changing one variable by
one unit, combining the direct path and all indirect paths (paths through
mediating variables).

```{r total_effects}
total_effects <- x1k$data |>
  estimate_all_total_effects(model)

round(total_effects, 3)
```


### Comparison with Multiple Regression Coefficients

Multiple regression coefficients and total causal effects do not agree when
mediating variables are present.

In the true causal structure of `generate_lingam_sample_6()`, there are two
paths from x3 to x1 (there is no **direct** edge from x3 to x1).

- x3 -> x0 -> x1 (indirect effect: 3.0 x 3.0 = **9.0**)
- x3 -> x2 -> x1 (indirect effect: 6.0 x 2.0 = **12.0**)
- **Total causal effect of x3 on x1 = 9.0 + 12.0 = 21.0**

We compare the coefficients from an OLS regression that includes all variables
to predict x1 against the results of `estimate_all_total_effects()`.

```{r total_vs_regression}
# Multiple regression: include all variables to predict x1
lm_coefs <- coef(lm(x1 ~ ., data = x1k$data))

# Comparison (variables causally related to x1: x0, x2, x3)
data.frame(
  variable           = c("x0", "x2", "x3"),
  OLS_coefficient    = round(lm_coefs[c("x0", "x2", "x3")], 3),
  total_causal_effect = round(total_effects["x1", c("x0", "x2", "x3")], 3)
)
```

The OLS coefficient for x3 is nearly **0**. This is because including x0 and x2
(the mediating variables) in the model causes x3's "effect through mediation"
to be absorbed into the coefficients of x0 and x2.

In contrast, the value of x3 from `estimate_all_total_effects()` is **~21**,
which correctly represents how much x1 ultimately changes when x3 is moved by
one unit.

| Question | Metric to use |
|------|------------|
| "How does x1 change if I move x3 while holding x0 and x2 fixed?" | OLS multiple regression coefficient |
| "How does x1 change if I move x3, through all paths?" | Total causal effect |

When you want to know "the ultimate impact of intervening on a variable," use
the total causal effect rather than the multiple regression coefficient.


## Inference with Prior Knowledge

With `make_prior_knowledge()`, you can incorporate domain knowledge about the
causal relationships among variables into Direct LiNGAM. This narrows the search
space and stabilizes estimation.

### Format of the Prior Knowledge Matrix

`make_prior_knowledge()` returns a $p \times p$ integer matrix. It uses the
indexing convention **row = effect variable (to), column = cause variable (from)**,
the same convention as the adjacency matrix.

| Value | Meaning |
|----|------|
| `-1` | Unknown (default; Direct LiNGAM searches freely) |
| `0` | This edge does not exist |
| `1` | This edge definitely exists |

The following shows how each argument affects the matrix.

| Argument | Value set | Meaning |
|------|------------|------|
| `exogenous_variables` | the entire **row** of the specified variable -> `0` | Receives no influence from any variable (root variable) |
| `sink_variables` | the entire **column** of the specified variable -> `0` | Exerts no influence on any variable (sink variable) |
| `paths` | `pk[to, from] = 1` | Specifies that this edge exists |
| `no_paths` | `pk[to, from] = 0` | Specifies that this edge does not exist |

Variables can be specified either by **1-based index** or by **variable name**
(which requires the `labels` argument).

### Usage Example

We supply domain knowledge about the true structure of `generate_lingam_sample_6()`.

- **x3** (index 4) is exogenous -- it receives no influence from any other variable
- **x1, x4, x5** (indices 2, 5, 6) are sink variables -- they exert no influence on other variables
- **Between x0 and x2** there is no path (in either direction)

#### Specifying by Index

```{r make_prior_knowledge1_index}
pk1 <- make_prior_knowledge(
  n_variables         = 6,
  exogenous_variables = 4,          # x3
  sink_variables      = c(2, 5, 6), # x1, x4, x5
  no_paths            = list(c(3, 1), c(1, 3)) # no x2<->x0
)

pk1
```

How to read the matrix: if `pk1["x1", "x3"]` is `-1`, then "x3->x1 is unknown
(LiNGAM searches for it)"; if `0`, then "x3->x1 does not exist".

#### Specifying by Variable Name

Passing `labels` lets you specify by variable name. This improves readability
and is robust to adding or reordering columns.

```{r make_prior_knowledge1_name}
pk1_named <- make_prior_knowledge(
  n_variables         = 6,
  exogenous_variables = "x3",
  sink_variables      = c("x1", "x4", "x5"),
  no_paths            = list(c("x2", "x0"), c("x0", "x2")),
  labels              = colnames(x1k$data)
)

# Equivalent in content to pk1
identical(pk1, pk1_named)
```

### Running Direct LiNGAM with Prior Knowledge

Simply pass it to the `prior_knowledge` argument and it is reflected in the search.

```{r LiNGAM_with_pk_1}
model_pk1 <- x1k$data |>
  lingam_direct(prior_knowledge = pk1, lambda = "BIC")

cat("Causal Order: ", colnames(x1k$data)[model_pk1$causal_order], "\n")
```

```{r LiNGAM_with_pk_2}
model_pk1$adjacency_matrix |>
  round(3)

model_pk1$adjacency_matrix |>
  plot_adjacency(
    labels    = colnames(model_pk1$adjacency_matrix),
    title     = "Estimated (with Prior Knowledge)",
    rankdir   = "TB",
    shape     = "circle",
    fillcolor = "lightgreen"
  )
```


## Choosing a Regression Method (reg_method)

In Direct LiNGAM, the adjacency matrix is estimated by regression after the
causal order is determined. The `reg_method` argument selects that regression method.

| `reg_method` | `glmnet` | Sparsification | Characteristics |
|---|---|---|---|
| `"ols"` | Not required | None | Estimates all edges. For sanity checks or environments without the package |
| `"lasso"` | Required | Yes | Shrinks weak edges to 0 |
| `"adaptive_lasso"` | Required | Yes (strong) | **Default**. Has the oracle property -- reliably sets truly zero edges to 0 |
| `"ridge"` | Required | None | Stabilizes coefficients with $\ell_2$ regularization. Robust to multicollinearity. Does not sparsify |

The oracle property is the theoretical guarantee that "the true structure can be
reliably recovered as the sample size grows," so `"adaptive_lasso"` is usually
recommended.

### Comparison of the Four Methods

```{r reg_method_compare}
fit_ols    <- lingam_direct(x1k$data, reg_method = "ols")
fit_lasso  <- lingam_direct(x1k$data, reg_method = "lasso",          lambda = "BIC")
fit_alasso <- lingam_direct(x1k$data, reg_method = "adaptive_lasso", lambda = "BIC")
fit_ridge  <- lingam_direct(x1k$data, reg_method = "ridge",          lambda = "BIC")

# Compare the adjacency matrices side by side
round(fit_ols$adjacency_matrix,    3)
round(fit_lasso$adjacency_matrix,  3)
round(fit_alasso$adjacency_matrix, 3)
round(fit_ridge$adjacency_matrix,  3)
```

OLS and Ridge tend to leave nonzero coefficients on all edges, whereas LASSO and
Adaptive LASSO shrink superfluous edges to 0. Ridge reduces the **magnitude** of
coefficients but does not set them to zero.

### Choosing lambda (common to LASSO / Adaptive LASSO)

The choice of penalty strength $\lambda$ directly determines the sparsity of the estimate.

| `lambda` | Method | Sparsity | Use |
|---|---|---|---|
| `"BIC"` | Information criterion | Highest | **Default**. Stable even with small samples |
| `"AIC"` | Information criterion | High | Leaves slightly more edges than BIC |
| `"lambda.min"` | CV (minimum prediction error) | Low | Prioritizes predictive accuracy. More edges |
| `"lambda.1se"` | CV (1SE rule) | Medium to high | Robust CV variant |
| `"oracle"` | Analytic formula (adaptive_lasso only) | - | $\lambda = 5 / n^{1.75}$. Guarantees the theoretical oracle property |

```{r reg_method_lambda}
# Compare BIC (default, sparsest) and lambda.min (minimum prediction error)
fit_bic     <- lingam_direct(x1k$data, lambda = "BIC")
fit_lam_min <- lingam_direct(x1k$data, lambda = "lambda.min")

# Number of nonzero edges
sum(fit_bic$adjacency_matrix     != 0)
sum(fit_lam_min$adjacency_matrix != 0)
```


## Independence between Error Variables

LiNGAM assumes that the residuals are independent.
`get_error_independence_p_values()` returns the p-values from tests of
independence between the residuals.

```{r Independence_between_error_variables1}
result <- x1k$data |>
  lingam_direct()

p_vals <- x1k$data |>
  get_error_independence_p_values(result)

round(p_vals, 3)
```


## The Non-Gaussianity Assumption

The theoretical heart of LiNGAM is the assumption that **the error terms follow
a non-Gaussian distribution**. When the errors are Gaussian, the **direction** of
causation becomes fundamentally unidentifiable (a reverse-direction model that
explains the same distribution exists), and the estimates are unreliable.

We verify this difference in practice by switching the error distribution with
the `noise_dist` argument of `generate_lingam_sample_6()`. The true structure is
as follows (the root is x3).

```{r nongauss_true}
set.seed(0)
truth <- generate_lingam_sample_6(noise_dist = "uniform")

truth$true_adjacency |>
  round(1)
```

The causal graph of the true structure:

```{r nongauss_true_plot}
truth$true_adjacency |>
  plot_adjacency(
    labels = colnames(truth$data),
    title  = "True structure"
  )
```

### Non-Gaussian Errors (Uniform Distribution) -- When It Works

```{r nongauss_uniform}
fit_uniform <- lingam_direct(truth$data)

# Estimated causal order (the true root x3 comes first)
colnames(truth$data)[fit_uniform$causal_order]

# The estimated adjacency matrix recovers the true structure almost perfectly
fit_uniform$adjacency_matrix |>
  round(1)
```

The estimated graph matches the true structure. Edges are color-coded against
the truth: green = correct, red = false positive, orange dashed = missed.

```{r nongauss_uniform_plot}
fit_uniform$adjacency_matrix |>
  plot_adjacency(
    labels = colnames(truth$data),
    true_B = truth$true_adjacency,
    title  = "Estimated (uniform errors)"
  )
```

### Gaussian Errors -- When It Fails

With the same causal structure, the results break down when the errors are Gaussian.

```{r nongauss_gaussian}
gauss <- generate_lingam_sample_6(noise_dist = "gaussian")
fit_gauss <- lingam_direct(gauss$data)

# The causal order does not match the true structure (root x3 does not come first)
colnames(gauss$data)[fit_gauss$causal_order]

fit_gauss$adjacency_matrix |>
  round(1)
```

Compared with the true structure, many edges are wrong (red) or missed (orange
dashed) -- the same color coding as above:

```{r nongauss_gaussian_plot}
fit_gauss$adjacency_matrix |>
  plot_adjacency(
    labels = colnames(gauss$data),
    true_B = truth$true_adjacency,
    title  = "Estimated (Gaussian errors)"
  )
```

With non-Gaussian errors the true adjacency matrix is recovered as-is, whereas
with Gaussian errors both the causal order and the coefficients deviate greatly
from the true structure. This is why it is said that "LiNGAM exploits the
non-Gaussianity of the data to determine the direction of causation." When
applying it to real data, it is important to **test the normality of the
residuals**, as in the next section, to check whether this assumption holds.


## Testing the Normality of Residuals

We test the normality of the residuals. Because LiNGAM assumes non-Gaussianity,
having normality **rejected** (a small p-value) is consistent with the model's
assumptions.

```{r nortest}
# Shapiro-Wilk (default)
x1k$data |>
  test_residual_normality(result)
```

We also check the normality of the residuals with a QQ plot.

```{r qqplot}
x1k$data |>
  plot_residual_qq(result)
```


## Model Summary

`summary_lingam()` runs the residual independence test and the normality test
together, letting you review at a glance how well the two assumptions LiNGAM
relies on hold (that the residuals are mutually independent, and that the
residuals are non-Gaussian). Instead of calling
`get_error_independence_p_values()` and `test_residual_normality()` separately,
you can survey the diagnostics in one place.

```{r summary_lingam}
x1k$data |>
  summary_lingam(result)
```


## Bootstrap Direct LiNGAM
We assess the reliability of the model using the bootstrap method.

```{r lingam_direct_bootstrap}
bs_model <- x1k$data |>
  lingam_direct_bootstrap(n_sampling = 100L, seed = 42)

bs_model
```

When the number of iterations or variables is large, specifying
`parallel = TRUE` lets it run faster on multiple cores. The number of cores is
specified with `n_cores` (when unspecified, it is capped at 2 cores for safety).

```{r lingam_direct_bootstrap_parallel, eval = FALSE}
bs_model <- x1k$data |>
  lingam_direct_bootstrap(
    n_sampling = 100L,
    seed       = 42,
    parallel   = TRUE,
    n_cores    = 4L
  )
```

Note that parallel execution uses L'Ecuyer's parallel random number streams, so
results are reproducible given the same `seed` and the same `n_cores`, but they
will not numerically match the results of sequential execution
(`parallel = FALSE`).


### Inspecting the Bootstrap Results

From the bootstrap results, we compute the frequency of occurrence of each path
and the mean of the coefficients.

```{r get_causal_direction_counts}
bs_model |>
  get_causal_direction_counts(labels = names(x1k$data))
```


### Adjacency Matrix of Mean Causal Effects

We construct an adjacency matrix from the bootstrap results.

```{r get_adjacency_matrix_summary}
bs_adjacency_matrix <- bs_model |>
  get_adjacency_matrix_summary(stat = "median")

bs_adjacency_matrix |>
  round(3)
```

We visualize the estimated adjacency matrix.
```{r plot_bootstrap}
bs_adjacency_matrix |>
  plot_adjacency(
    labels    = colnames(x1k$data),
    title     = "Estimated (with Bootstrap)",
    rankdir   = "TB",
    shape     = "circle",
    fillcolor = "lightgreen"
  )
```


### Matrix of Path Occurrence Frequencies

We compute the matrix of occurrence frequencies for each path.

```{r get_probabilities}
bs_model |>
  get_probabilities()
```


### Mean Total Effects

We compute the mean total effect of each path.

```{r get_total_causal_effects}
bs_model |>
  get_total_causal_effects()
```

We turn the bootstrap results into a causal graph. By default, only paths that
occur in at least 50% of samples are shown.

```{r plot_bootstrap_probabilities}
bs_model |>
  plot_bootstrap_probabilities()
```


### Stability of the Causal Order

`get_causal_order_stability()` aggregates the causal orders estimated in each
bootstrap sample and quantifies how stable the order is. It returns the rank
distribution of each variable, the precedence probability for variable pairs
(`P[i, j]` = the fraction of samples in which variable i came upstream of j),
and an overall stability score (0 = random, 1 = identical across all samples).

```{r get_causal_order_stability}
bs_model |>
  get_causal_order_stability(labels = names(x1k$data))
```


## Integration with broom (tidy / glance)

Estimation results can be converted to a data.frame with the `broom`-compatible
`tidy()` / `glance()`, making integration with `ggplot2` and `dplyr` easy.
`tidy()` returns an edge list (`from`, `to`, `estimate`), and `glance()` returns
a one-row summary of the whole model. `tidy()` also works on bootstrap results,
in which case it returns the occurrence frequencies for each direction, etc.

```{r broom}
# Convert the estimated adjacency matrix to an edge list
tidy(model)

# One-row summary of the whole model
glance(model)

# Direction-wise summary of the bootstrap results (variable names via labels)
tidy(bs_model, labels = names(x1k$data))
```


## A Larger Dataset (10 Variables)

An example of a larger dataset with 10 variables and 10,000 rows.

```{r sample1_matrix10}
x10k <- generate_lingam_sample_10(n = 10000)

x10k$true_adjacency |>
  plot_adjacency(
    labels  = colnames(x10k$data),
    title   = "True causal structure",
    rankdir = "TB",
    shape   = "circle"
  )
```



## Comparing ICA-LiNGAM and Direct LiNGAM

`pcalg::lingam()` is the original LiNGAM algorithm, which estimates the mixing
matrix with FastICA and obtains the causal order and coefficients
(Shimizu et al. 2006). It solves the same problem while taking an approach
independent of `lingam_direct()`.

### Running Both Algorithms

We analyze the same 6-variable dataset ($n = 1000$) with both methods.

```{r ica_run, eval = requireNamespace("pcalg", quietly = TRUE)}
d_cmp <- generate_lingam_sample_6(n = 1000, seed = 42)

t_cmp_direct <- system.time(res_cmp_direct <- lingam_direct(d_cmp$data))
t_cmp_ica    <- system.time(res_cmp_ica    <- pcalg::lingam(as.matrix(d_cmp$data)))

cat(sprintf("Direct LiNGAM : %.2f sec\nICA-LiNGAM    : %.2f sec\n",
            t_cmp_direct["elapsed"], t_cmp_ica["elapsed"]))
```

### Comparing the Estimated Coefficients

`$Bpruned` uses the same convention as the lingamr adjacency matrix
(`B[i, j]` = coefficient of $x_j \to x_i$).

```{r ica_compare, eval = requireNamespace("pcalg", quietly = TRUE)}
B_ica <- res_cmp_ica$Bpruned
rownames(B_ica) <- colnames(B_ica) <- names(d_cmp$data)

idx_ica  <- which(abs(B_ica) > 0, arr.ind = TRUE)
tidy_ica <- data.frame(
  from  = colnames(B_ica)[idx_ica[, 2]],
  to    = rownames(B_ica)[idx_ica[, 1]],
  ica   = round(B_ica[idx_ica], 3)
)

tidy_dir <- tidy(res_cmp_direct)
tidy_dir <- data.frame(from = tidy_dir$from, to = tidy_dir$to,
                       direct = round(tidy_dir$estimate, 3))

merge(tidy_dir, tidy_ica, by = c("from", "to"), sort = TRUE)
```

### Comparing the DAG Structures

We compare the structures with a full outer join over all edges and check
consistency with the true DAG.

```{r ica_structure, eval = requireNamespace("pcalg", quietly = TRUE)}
B_true   <- d_cmp$true_adjacency
idx_true <- which(abs(B_true) > 0, arr.ind = TRUE)
true_key <- paste(colnames(B_true)[idx_true[, 2]],
                  rownames(B_true)[idx_true[, 1]], sep = "->")

cmp <- merge(tidy_dir, tidy_ica, by = c("from", "to"), all = TRUE, sort = TRUE)
cmp$truth <- paste(cmp$from, cmp$to, sep = "->") %in% true_key
cmp
```

When the `direct` or `ica` column is `NA`, it means that method did not detect
that edge. `truth = TRUE` indicates an edge that exists in the true DAG.

---

## When There Are Many Variables: The Scalability Wall

At each step, Direct LiNGAM performs independence tests on all remaining pairs
of variables. Since the number of steps is $p$ and the number of tests per step
is at most $p(p-1)$, the total number of independence tests is approximately

$$\sum_{k=1}^{p} k(k-1) \approx \frac{p^3}{3}$$

giving a computational cost of **$O(p^3)$**.
By contrast, the FastICA used by ICA-LiNGAM is $O(p^2 n)$ (with BLAS
optimization), so the gap widens as the number of variables grows.

`generate_lingam_large_sample()` generates random sparse DAG data with a freely
configurable number of variables `p`. Each variable $x_i$ ($i \ge 1$) randomly
has at most `max_parents` parents chosen from $x_0, \ldots, x_{i-1}$. Since the
causal order is guaranteed to follow the index order, the adjacency matrix is
always a **lower triangular matrix**.

### Generating the Data

```{r large_gen}
d20 <- generate_lingam_large_sample(p = 20, n = 1000, seed = 42)

dim(d20$data)                    # 1000 rows x 20 columns
sum(d20$true_adjacency != 0)     # number of true edges (sparse DAG)
d20$true_causal_order            # 0, 1, ..., 19
```

### Comparing Execution Times

When $p$ grows by a factor of 1.5 (10 -> 15), the number of independence tests
grows by a factor of $15^3 / 10^3 \approx 3.4$.

```{r large_timing}
d10 <- generate_lingam_large_sample(p = 10, n = 500, seed = 42)
d15 <- generate_lingam_large_sample(p = 15, n = 500, seed = 42)

t10 <- system.time({ r10 <- lingam_direct(d10$data) })
t15 <- system.time({ r15 <- lingam_direct(d15$data) })

cat(sprintf(
  "p = 10 : %.2f sec\np = 15 : %.2f sec\ntheoretical factor %.1fx vs. observed %.1fx\n",
  t10["elapsed"],
  t15["elapsed"],
  15^3 / 10^3,
  t15["elapsed"] / max(t10["elapsed"], 0.01)
))
```

We run ICA-LiNGAM on the same data to compare speed directly.

```{r large_timing_ica, eval = requireNamespace("pcalg", quietly = TRUE)}
t10_ica <- system.time({ pcalg::lingam(as.matrix(d10$data)) })
t15_ica <- system.time({ pcalg::lingam(as.matrix(d15$data)) })

cat(sprintf(
  "              p = 10   p = 15\nDirect LiNGAM : %5.2f sec  %5.2f sec\nICA-LiNGAM    : %5.2f sec  %5.2f sec\n",
  t10["elapsed"], t15["elapsed"],
  t10_ica["elapsed"], t15_ica["elapsed"]
))
```

The larger $p$ becomes, the more Direct LiNGAM's $O(p^3)$ cost dominates, and the
gap between the two widens. In large-scale settings such as $p = 30$ or $p = 50$,
this trend becomes even more pronounced.

### Checking Estimation Accuracy (p = 10)

Even with a sparse DAG, as long as there are **non-Gaussian errors** (default:
uniform distribution), Direct LiNGAM can recover the correct causal order.

```{r large_causal_order}
# Estimated causal order
r10$causal_order

# Whether it matches the true causal order 0, 1, ..., 9 exactly
all(r10$causal_order == d10$true_causal_order)
```

We convert to an edge list with `tidy()` and inspect the estimated coefficients.

```{r large_tidy}
tidy(r10) |>
  head(10)
```


## High-Dimensional Direct LiNGAM

The $O(p^3)$ independence-test cost shown above becomes a real bottleneck once
$p$ grows into the tens or hundreds, and breaks down entirely once $p > n$
(more variables than observations), where the usual regression-based
adjacency estimation is no longer well defined.

`lingam_high_dim()` implements HighDimDirectLiNGAM (Wang & Drton 2020), a
variant designed for this regime. Instead of pairwise independence tests, it
searches the causal order using moment statistics of non-Gaussianity, computed
from a cached Gram matrix. The algorithm is deterministic (no random
restarts), and it returns the same `LingamResult` object as `lingam_direct()`,
so `print()`, `tidy()`, `plot_adjacency()` and `estimate_total_effect()` all
work unchanged.

```{r high_dim_basic}
hd_sample <- generate_lingam_sample_6(n = 500, seed = 1)
hd_result <- lingam_high_dim(hd_sample$data)

hd_result$causal_order
round(hd_result$adjacency_matrix, 3)
```

When `n_samples <= n_features`, the usual BIC-based Adaptive LASSO cannot be
used to estimate the adjacency matrix, so `lingam_high_dim()` falls back to a
cross-validated LASSO (`glmnet::cv.glmnet`) and emits a warning:

```{r high_dim_wide}
wide_sample <- generate_lingam_large_sample(p = 30, n = 25, seed = 1)
wide_result <- lingam_high_dim(wide_sample$data)

wide_result$causal_order
```

## A Case Where DirectLiNGAM Struggles: The Measurement Error Paradox

Causal discovery methods have assumptions, and when those are violated they may
fail to recover the correct structure. `generate_lingam_paradox_data()` is a
dataset designed to deliberately create such a difficult case. As with the other
sample generators, it returns a list containing `data` and `true_adjacency`.

The true structure of this data is a simple serial chain
**x0 -> x1 -> x2 -> x3** (each coefficient 0.8). However, it has two notable features.

- **Heavy measurement error is added to the root variable x0.** This disrupts the
  independence assessment performed in DirectLiNGAM's first step, causing it to
  choose the root incorrectly and making error propagation more likely.
- All variables are **standardized** with `scale()` (no differences in scale).

```{r paradox_generate}
paradox <- generate_lingam_paradox_data(n = 2000L, seed = 42)

head(paradox$data)

# All variables are standardized (sd = 1)
sapply(paradox$data, sd)
```

We visualize the true causal graph. The coefficient 0.8 is the structural
coefficient on the latent scale before standardization.

```{r paradox_true}
paradox$true_adjacency |>
  plot_adjacency(
    labels  = colnames(paradox$true_adjacency),
    title   = "True causal chain (x0 -> x1 -> x2 -> x3)",
    rankdir = "LR",
    shape   = "circle"
  )
```

Now let us apply Direct LiNGAM.

```{r paradox_fit}
model_p <- lingam_direct(paradox$data)

# Estimated causal order
colnames(paradox$data)[model_p$causal_order]
```

Note that the **head of the estimated causal order is x1, not the true root x0**.
Because of the measurement error on the root, DirectLiNGAM fails to select x0 as
the first exogenous variable.

```{r paradox_adjacency}
model_p$adjacency_matrix |>
  round(3)

model_p$adjacency_matrix |>
  plot_adjacency(
    labels    = colnames(model_p$adjacency_matrix),
    title     = "Estimated structure (paradox data)",
    rankdir   = "LR",
    shape     = "circle",
    fillcolor = "lightpink"
  )
```

While the downstream **x1 -> x2 -> x3** is recovered correctly, the **direction
between x0 and x1 is reversed** (the truth is x0 -> x1, but the estimate is
x1 -> x0), and x0 ends up being treated almost like a sink.

We use the bootstrap to check whether this error occurred by chance or is systematic.

```{r paradox_bootstrap}
bs_paradox <- paradox$data |>
  lingam_direct_bootstrap(n_sampling = 100L, seed = 42)

# Occurrence probability of each direction (row = to, column = from)
bs_paradox |>
  get_probabilities() |>
  round(2)
```

The important point is that the incorrect direction **x1 -> x0** is reproduced
with nearly 100% probability. In other words, this error is not coincidental but
**systematic**, and it appears stably across bootstrap samples.

> **Lesson:** Bootstrap stability (high reproduction probability) does not
> guarantee the *correctness* of the estimate. When the model's assumptions (here,
> the assumption that "upstream variables have no measurement error") are violated,
> the method may **stably** recover an incorrect structure. It is important to
> evaluate results critically, together with tests of residual independence and
> normality and domain knowledge about the data-generating process.


## VAR-LiNGAM: Causal Discovery in Time Series

Direct LiNGAM assumes that observations are **independent and identically
distributed (i.i.d.)**, a requirement that time series data violate.
**VAR-LiNGAM** (Hyvärinen et al., 2010) handles stationary time series by
first fitting a Vector Autoregression (VAR) model to absorb temporal
autocorrelation, then applying Direct LiNGAM to the VAR residuals to recover
the **instantaneous** causal structure $B_0$.
The model is:

$$X_t = B_0\,X_t + \sum_{k=1}^{p} B_k\,X_{t-k} + e_t$$

where $B_0$ encodes contemporaneous causal effects (strictly acyclic),
$B_1, \ldots, B_p$ encode lagged effects, and $e_t$ are mutually independent
non-Gaussian disturbances.

### Sample Data

`generate_varlingam_sample()` produces a three-variable time series from a
VAR(1)-LiNGAM model.
The instantaneous structure is $x_0 \to x_1 \to x_2$ (coefficients 0.6 and
−0.5), and the only cross-variable lag-1 effect is $x_2(t-1) \to x_0(t)$
(coefficient 0.3).

```{r varlingam_sample}
s <- generate_varlingam_sample(n = 1000, seed = 42)

# True instantaneous coefficient matrix B0  (B0[i, j]: x_j -> x_i)
s$true_B0

# True lag-1 coefficient matrix  (M1[i, j]: x_j(t-1) -> x_i(t), structural)
s$true_M1
```

### Fitting VAR-LiNGAM

Pass the data matrix to `lingam_var()`.
Rows must be in chronological order (earliest first).

```{r lingam_var}
model <- lingam_var(s$data, lags = 1)
model
```

The result object contains `adjacency_matrices`, a three-dimensional array of
shape `[1 + lags, n_features, n_features]`:

- **`[1, , ]` (`"lag0"`):** instantaneous matrix $B_0$.
  `B0[i, j]` is the direct effect of $x_j$ on $x_i$ at the *same* time step.
- **`[k + 1, , ]` (`"lag`*k*`"`):** lagged matrix $B_k$.
  `Bk[i, j]` is the direct structural effect of $x_j(t-k)$ on $x_i(t)$.

Both can be extracted by their dimension label:

```{r varlingam_matrices}
B0 <- model$adjacency_matrices["lag0", , ]
B1 <- model$adjacency_matrices["lag1", , ]

round(B0, 2)  # compare with s$true_B0
round(B1, 2)  # compare with s$true_M1
```

### Lag Order Selection

By default, `lingam_var()` automatically selects the lag order among `1:lags`
using the Bayesian Information Criterion (`criterion = "bic"`).
The alternatives `"aic"`, `"hqic"`, and `"fpe"` are also supported.
To use a fixed lag order without automatic selection, set `criterion = NULL`:

```{r varlingam_lags, eval = FALSE}
# Fix lag order to 2 without IC-based selection
model_lag2 <- lingam_var(s$data, lags = 2, criterion = NULL)
```

### Stationarity Check

VAR-LiNGAM is defined for **stationary** processes.
`check_var_stationarity()` inspects the eigenvalues of the VAR companion matrix:
the process is stationary when all moduli are **strictly less than 1**.

```{r varlingam_stationarity}
check_var_stationarity(model)
```

A `max_modulus` at or above 1 indicates a unit-root or explosive process.
In that case, differencing the series before analysis is recommended.

### Residual Diagnostics

LiNGAM assumes that the error terms $e_t$ are **non-Gaussian**.
`test_varlingam_residual_normality()` tests whether the LiNGAM innovations
$e_t = (I - B_0)\,n_t$ (where $n_t$ are the stored VAR residuals) depart
from normality.
A small p-value (reject $H_0$: Gaussian) supports the model assumption.

```{r varlingam_normality}
test_varlingam_residual_normality(model)
```

`test_varlingam_residual_normality_all()` runs several tests at once and
appends skewness and excess kurtosis columns for a quick overview:

```{r varlingam_normality_all}
test_varlingam_residual_normality_all(model, methods = c("shapiro", "jb"))
```

`plot_varlingam_residual_qq()` draws per-variable normal Q-Q plots.
Deviations from the straight reference line indicate non-Gaussianity.

```{r varlingam_qq, fig.width = 7, fig.height = 3}
plot_varlingam_residual_qq(model)
```

### Total Causal Effects

`estimate_var_total_effect()` estimates the **total** causal effect of one
variable on another, integrating over all direct and mediated paths.
The `from_lag` argument controls the time offset of the cause:
`from_lag = 0` (default) gives the contemporaneous total effect;
`from_lag = 1` gives the one-step-ahead effect of $x_j(t-1)$ on $x_i(t)$.

```{r varlingam_total_effect}
# Total effect x0 -> x2 (contemporaneous)
estimate_var_total_effect(s$data, model, from_index = 1, to_index = 3)

# Total effect x0(t-1) -> x2(t) (one-step-ahead)
estimate_var_total_effect(s$data, model, from_index = 1, to_index = 3, from_lag = 1)
```

Variable indices are **1-based** integers, or column names as character strings.

### Bootstrap

`lingam_var_bootstrap()` quantifies the uncertainty of the estimated structure
by re-running VAR-LiNGAM on **residual bootstrap** samples.
Unlike the Direct LiNGAM bootstrap (which resamples i.i.d. rows),
VAR-LiNGAM holds the fitted values fixed and resamples only the VAR residuals,
preserving the temporal structure of the series.

```{r varlingam_bootstrap}
bs_var <- lingam_var_bootstrap(
  s$data,
  n_sampling = 100L,
  seed       = 42,
  verbose    = FALSE
)
```

`get_var_probabilities()` returns the proportion of bootstrap samples in which
each directed edge was detected.
The column layout mirrors `adjacency_matrices`: the first `n_features` columns
correspond to the instantaneous structure (lag 0), the next `n_features` to
lag 1, and so on.

```{r varlingam_probabilities}
round(get_var_probabilities(bs_var), 2)
```

`get_var_paths()` enumerates all causal paths between two variables found
across bootstrap samples, together with each path's average total effect and
detection probability.

```{r varlingam_paths}
# Paths from x0 to x2 at the same time step (from_lag = 0)
get_var_paths(bs_var, from_index = 1, to_index = 3)
```

```{r varlingam_paths_lag}
# Paths from x0(t-1) to x2(t)  (from_lag = 1)
get_var_paths(bs_var, from_index = 1, to_index = 3, from_lag = 1)
```


## LiNGAM for Mixed Data (LiM)

Direct LiNGAM assumes every variable is continuous. `lingam_lim()` relaxes
this assumption and estimates a causal structure from data containing a
mixture of continuous and binary (0/1) discrete variables, following
Zeng et al. (2022). It combines a NOTEARS-style continuous optimization
(the "global" phase) with a combinatorial local search over edge
directions, pruning, and edge addition (the "local" phase).

`generate_lim_sample()` generates a small dataset with a known causal chain
of continuous and discrete variables: `x1` (continuous) -> `x2` (discrete)
-> `x3` (continuous).

```{r lim_sample}
set.seed(1)
lim_dat <- generate_lim_sample(n = 2000)
head(lim_dat$data)
lim_dat$is_continuous
```

`lingam_lim()` requires `is_continuous`, a logical vector marking which
columns are continuous (`TRUE`) versus binary discrete (`FALSE`).
Because the optimization starts from a random initial point,
reproducibility requires `set.seed()`.

```{r lim_fit}
lim_result <- lingam_lim(lim_dat$data, is_continuous = lim_dat$is_continuous)
print(lim_result)
```

As with `lingam_direct()`, `adjacency_matrix` follows the
`B[i, j]` = j -> i convention (row = to, column = from), and `causal_order`
lists the estimated topological order as 1-based indices.

```{r lim_order}
colnames(lim_dat$data)[lim_result$causal_order]
```

Only binary (0/1) discrete variables are supported; see `?lingam_lim` for
details on the local phase's edge-weight convention and its numeric
differences from the Python implementation.


## Multi-Group Direct LiNGAM

`lingam_direct()` fits a single dataset. When data comes from several
sources that plausibly share the same causal structure but not the same
strength of effect (e.g. the same study run at multiple sites, or the same
process observed in different time periods), `lingam_multi_group()` jointly
estimates a **common causal order** across all groups while still allowing
each group its own adjacency matrix (structural coefficients), following
Shimizu (2012).

`generate_multi_group_sample()` generates two datasets that share the causal
structure of `generate_lingam_sample_6()` but with slightly different
coefficients per group.

```{r multi_group_sample}
mg <- generate_multi_group_sample(n = c(1000, 1000), seed = 42)
lapply(mg$data_list, head, 3)
```

```{r multi_group_fit}
mg_result <- lingam_multi_group(mg$data_list, reg_method = "ols")
print(mg_result)
```

`causal_order` is shared by all groups; `adjacency_matrices` holds one
matrix per group, each following the usual `B[i, j]` = j -> i convention.

To analyze a single group with the rest of `lingamr`'s single-group
tooling (total causal effects, independence tests, plotting), extract it
as a plain `LingamResult` with `get_group_result()`:

```{r multi_group_extract}
g1 <- get_group_result(mg_result, "group1")
class(g1)

estimate_all_total_effects(mg$data_list$group1, g1, method = "ols")
```

`lingam_multi_group_bootstrap()` provides bootstrap stability estimates in
the same joint fashion: every iteration resamples each group independently,
then jointly re-estimates the causal order and per-group adjacency
matrices. It returns a named list of per-group `BootstrapResult` objects, so
the existing bootstrap query functions apply directly per group:

```{r multi_group_bootstrap}
mg_bs <- lingam_multi_group_bootstrap(mg$data_list,
  n_sampling = 20L, reg_method = "ols", seed = 1, verbose = FALSE
)
get_probabilities(mg_bs$group1)
```

Note that `lingam_multi_group_bootstrap()`'s total causal effects are
computed as path-coefficient products over each iteration's adjacency
matrix, not via regression; this matches the upstream Python
implementation but differs from `lingam_direct_bootstrap()`'s
regression-based `estimate_all_total_effects()`.


## Causal Discovery with Missing Data

All algorithms above assume a complete data matrix. When `X` contains missing
values (`NA`), `bootstrap_with_imputation()` combines bootstrap resampling
with multiple imputation: each resample is imputed into several complete
datasets, and a common causal structure is jointly estimated across them with
`lingam_multi_group()` (treating the imputed copies as "groups" that share one
causal order). This is an R port of the Python
`lingam.tools.bootstrap_with_imputation()`.

```{r bootstrap_with_imputation_setup, eval = requireNamespace("mice", quietly = TRUE)}
sample6_na <- generate_lingam_sample_6(n = 1000, seed = 1)
X_na <- sample6_na$data
set.seed(1)
X_na$x5[sample.int(nrow(X_na), size = round(0.1 * nrow(X_na)))] <- NA # MCAR 10% on x5
```

```{r bootstrap_with_imputation_fit, eval = requireNamespace("mice", quietly = TRUE)}
bwi <- bootstrap_with_imputation(X_na,
  n_sampling = 20L, n_repeats = 5L, seed = 42, verbose = FALSE
)
print(bwi)
```

The default imputer is `mice::mice(method = "norm")` (Bayesian linear
regression), the closest standard R equivalent of the upstream Python default
(`IterativeImputer(sample_posterior = TRUE)`); numeric results will not match
the Python implementation. Both the imputer and the causal-discovery step can
be swapped for a custom `function` via the `imputer` and `cd_fit` arguments.

Because each iteration produces `n_repeats` adjacency matrices (one per
imputed dataset), the result's shape differs from `lingam_direct_bootstrap()`.
`as_bootstrap_result()` collapses the `n_repeats` dimension (median or mean)
into a regular `BootstrapResult`, so the existing bootstrap query functions
apply as usual:

```{r bootstrap_with_imputation_convert, eval = requireNamespace("mice", quietly = TRUE)}
bs_na <- as_bootstrap_result(bwi, aggregate = "median")
get_probabilities(bs_na)
```

`get_total_causal_effects()` is not available on this `BootstrapResult`, since
`bootstrap_with_imputation()` never computes total effects.


## Latent Confounders: Bottom-Up ParceLiNGAM

Every algorithm above assumes there is no latent (unobserved) confounder: any
variable that causes two or more of the observed variables must itself be
observed. When that assumption fails, `lingam_direct()` will still return a
full causal order, but silently -- some part of it may be wrong, with no
indication of which part.

`lingam_parce()` (BottomUpParceLiNGAM, Tashiro et al. 2014) is designed for
this situation. It searches for the causal order from the sink (most
downstream) side, testing at each step whether a candidate variable's
residual is independent of the others. As soon as that test is rejected, the
search stops, and every variable it could not yet place is returned together
as a single **unresolved block** -- a signal that those variables likely
share a latent confounder, rather than a (possibly wrong) guess at their
order.

`generate_parce_sample()` generates a 7-variable model in which `x6` is an
unobserved common cause of `x2` and `x3`; only `x0`-`x5` are returned as data.

```{r parce_sample}
# HSIC is O(n^2), so a moderate n keeps this vignette fast to build
confounded <- generate_parce_sample(n = 500, seed = 1)
head(confounded$data)
confounded$confounded_pair
```

```{r parce_fit}
parce_result <- lingam_parce(confounded$data, reg_method = "ols")
print(parce_result)
```

The causal order's first element is the unresolved block, shown in
parentheses; here it correctly contains `x2` and `x3`. The corresponding
entries of the adjacency matrix are `NA`, while edges among the remaining,
fully-resolved variables are estimated as usual:

```{r parce_block}
parce_result$causal_order[[1]]
parce_result$adjacency_matrix[confounded$confounded_pair, confounded$confounded_pair]
```

Because a confounded variable's true parents cannot be identified,
`estimate_total_effect_parce()` warns and returns `NA` when asked for a total
effect *from* a variable in the unresolved block, but still computes normal
estimates for well-identified pairs:

```{r parce_total_effect}
# from a confounded variable: warns and returns NA
estimate_total_effect_parce(confounded$data, parce_result,
  from_index = confounded$confounded_pair[1], to_index = "x1"
)

# a well-identified pair: a normal numeric estimate
estimate_total_effect_parce(confounded$data, parce_result,
  from_index = "x0", to_index = "x5"
)
```

`lingam_parce_bootstrap()` provides bootstrap stability estimates in the same
style as `lingam_direct_bootstrap()`. `NA` (unresolved) edges are treated as
absent when aggregating, so `get_probabilities()` and the other
`BootstrapResult` query functions work as usual; `get_causal_order_stability()`
is the one exception, since ParceLiNGAM's blocked causal order does not fit
its fixed-length format.

```{r parce_bootstrap}
parce_bs <- lingam_parce_bootstrap(confounded$data,
  n_sampling = 10L, reg_method = "ols", seed = 1, verbose = FALSE
)
get_probabilities(parce_bs)
```


## Latent Confounders: RCD

`lingam_rcd()` (Repetitive Causal Discovery; Maeda and Shimizu 2020) tackles
the same latent-confounder problem as `lingam_parce()`, but from a different
angle: rather than searching for a causal order and giving up on an
**unresolved block** once a test is rejected, RCD directly estimates each
variable's **ancestor set** and then checks individual, parent-free pairs for
a shared latent confounder. This makes RCD's output pair-level (which specific
pairs are confounded) rather than block-level (which set of variables could
not be ordered).

`generate_rcd_sample()` generates a 7-variable model in which `x6` is an
unobserved common cause of `x2` and `x4`; only `x0`-`x5` are returned as data.

```{r rcd_sample}
# HSIC is O(n^2), so a moderate n keeps this vignette fast to build
rcd_confounded <- generate_rcd_sample(n = 300, seed = 1)
head(rcd_confounded$data)
rcd_confounded$confounded_pair
```

```{r rcd_fit}
rcd_result <- lingam_rcd(rcd_confounded$data)
print(rcd_result)
```

`ancestors_list` gives each variable's estimated ancestors (not a causal
order), and the confounded pair's adjacency-matrix entries are `NA`:

```{r rcd_ancestors}
rcd_result$ancestors_list
rcd_result$adjacency_matrix[rcd_confounded$confounded_pair, rcd_confounded$confounded_pair]
```

As with ParceLiNGAM, `estimate_total_effect_rcd()` warns and returns `NA`
when asked for a total effect *from* a confounded variable:

```{r rcd_total_effect}
# from a confounded variable: warns and returns NA
estimate_total_effect_rcd(rcd_confounded$data, rcd_result,
  from_index = rcd_confounded$confounded_pair[1], to_index = rcd_confounded$confounded_pair[2]
)

# a well-identified pair: a normal numeric estimate
estimate_total_effect_rcd(rcd_confounded$data, rcd_result,
  from_index = "x5", to_index = "x0"
)
```


## Evaluating Model Fit

`evaluate_model_fit()` treats an estimated adjacency matrix as a structural
equation model (SEM) and reports standard SEM fit measures (CFI, RMSEA,
AIC/BIC, etc.) via the `lavaan` package (an optional dependency; install it
with `install.packages("lavaan")`). This is useful for judging whether an
estimated causal graph is consistent with the data, independent of how it
was estimated.

```{r evaluate_model_fit, eval = requireNamespace("lavaan", quietly = TRUE)}
sample6 <- generate_lingam_sample_6()
fit_result <- lingam_direct(sample6$data, reg_method = "ols")

# fit measures for the estimated graph
evaluate_model_fit(fit_result, sample6$data)
```

Reversing the direction of every edge produces a mis-specified model, and
its fit measures are visibly worse (lower CFI, higher RMSEA):

```{r evaluate_model_fit_reversed, eval = requireNamespace("lavaan", quietly = TRUE)}
reversed_adjacency <- t(fit_result$adjacency_matrix)
evaluate_model_fit(reversed_adjacency, sample6$data)
```


## When LiNGAM Cannot Be Used

LiNGAM (and `lingamr`) requires several assumptions.
When these are not met, estimation either fails or systematically recovers an
incorrect structure.

| Assumption | When problems arise | Remedy / alternative |
|---|---|---|
| **Non-Gaussian errors** | When all errors follow a Gaussian distribution, the causal direction becomes unidentifiable | See the "The Non-Gaussianity Assumption" section of this vignette. ICA-LiNGAM and Direct LiNGAM fail equally |
| **Acyclic graph (DAG)** | When feedback loops (x -> y -> x) exist | Consider Cyclic LiNGAM (implemented in the Python version) |
| **No latent common causes** | When unobserved common causes (hidden confounders) exist | Consider LvLiNGAM (Latent variable LiNGAM) |
| **Linear causal relationships** | When the relationships among variables are nonlinear | Consider additive noise models (ANM) or nonlinear ICA |
| **No measurement error (upstream variables)** | When heavy measurement error is present on variables near the root, the direction is systematically reversed | See the "The Measurement Error Paradox" section of this vignette |
| **Independent and identically distributed (i.i.d.)** | When there is time-series data, hierarchical data, or cluster structure | Consider VAR-LiNGAM (time series), MultiBench (multi-domain), etc. |
| **Sufficient sample size** | When $n$ is extremely small relative to the number of variables $p$ (rule of thumb: $n < 10p$), estimation tends to be unstable | Reduce the number of variables; sparsify with `reg_method = "adaptive_lasso"` |

### A Checklist to Verify in Advance

Before starting an actual analysis, we recommend confirming the following.

1. **Acyclicity of the graph** -- Can feedback loops be ruled out from domain expertise?
2. **Absence of latent variables** -- Are the key observed variables all present?
3. **Non-Gaussianity of the errors** -- Can be checked with `test_residual_normality()`
   (though this is a post-estimation diagnostic). As a quick check beforehand,
   visually inspect each variable's histogram and skewness.
4. **Presence of measurement error** -- Is there measurement error on variables near
   the root? If so, interpret with care.
5. **Sample size** -- Aim for $n \geq 10p$. If it falls short, do not over-trust the results.

> **Summary:** LiNGAM is powerful when all five assumptions -- linear, acyclic,
> non-Gaussian, no latent variables, and i.i.d. -- hold.
> Verifying these with domain knowledge and residual diagnostics before analysis
> is the first step toward reliable causal inference.
