Skip to contents

Causal Discovery with lingamr

Source: vignettes/lingamr.Rmd
lingamr.Rmd

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

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)
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)
generate_lingam_paradox_data() 4 2,000 The measurement error paradox (used in The Paradox Example)

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.

x1k <- generate_lingam_sample_6(n = 1000)

x1k$data |>
  head()
#>         x0        x1       x2        x3        x4        x5
#> 1 2.814924 18.017120 4.543655 0.6333728 18.160090 12.236660
#> 2 1.889685 10.956005 2.188091 0.3175366 13.172754  7.932657
#> 3 1.008905  6.990652 1.953131 0.2409218  6.702107  4.797122
#> 4 1.965690 12.296763 2.847148 0.3784141 13.224002  8.685252
#> 5 1.698178  9.698147 2.145058 0.3521443 11.673495  7.366258
#> 6 1.412372  8.640107 1.929980 0.2977585 10.024075  6.340899
x1k$true_adjacency
#>    x0 x1 x2 x3 x4 x5
#> x0  0  0  0  3  0  0
#> x1  3  0  2  0  0  0
#> x2  0  0  0  6  0  0
#> x3  0  0  0  0  0  0
#> x4  8  0 -1  0  0  0
#> x5  4  0  0  0  0  0

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

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.

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.

# index number
model$causal_order
#> [1] 4 3 1 5 6 2

# variable name
colnames(x1k$data)[model$causal_order]
#> [1] "x3" "x2" "x0" "x4" "x5" "x1"

Estimated Adjacency Matrix

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

model$adjacency_matrix |>
  round(3)
#>       x0 x1     x2    x3 x4 x5
#> x0 0.000  0  0.000 3.033  0  0
#> x1 2.988  0  2.002 0.000  0  0
#> x2 0.000  0  0.000 5.993  0  0
#> x3 0.000  0  0.000 0.000  0  0
#> x4 8.000  0 -1.000 0.000  0  0
#> x5 4.015  0  0.000 0.000  0  0

Drawing the Causal Graph

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

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)
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).

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).

total_effects <- x1k$data |>
  estimate_all_total_effects(model)

round(total_effects, 3)
#>       x0 x1     x2     x3 x4 x5
#> x0 0.000  0  0.000  3.033  0  0
#> x1 2.872  0  1.937 21.059  0  0
#> x2 0.000  0  0.000  5.993  0  0
#> x3 0.000  0  0.000  0.000  0  0
#> x4 7.910  0 -1.129 18.276  0  0
#> x5 4.015  0  0.000 12.179  0  0

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().

# 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)
)
#>    variable OLS_coefficient total_causal_effect
#> x0       x0           3.237               2.872
#> x2       x2           1.965               1.937
#> x3       x3           0.014              21.059

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×pp \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 

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
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]   -1    0    0   -1    0    0
#> [2,]   -1   -1   -1   -1    0    0
#> [3,]    0    0   -1   -1    0    0
#> [4,]    0    0    0   -1    0    0
#> [5,]   -1    0   -1   -1   -1    0
#> [6,]   -1    0   -1   -1    0   -1

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.

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)
#> [1] FALSE

Running Direct LiNGAM with Prior Knowledge

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

model_pk1 <- x1k$data |>
  lingam_direct(prior_knowledge = pk1, lambda = "BIC")

cat("Causal Order: ", colnames(x1k$data)[model_pk1$causal_order], "\n")
#> Causal Order:  x3 x2 x0 x4 x5 x1
model_pk1$adjacency_matrix |>
  round(3)
#>       x0 x1     x2    x3 x4 x5
#> x0 0.000  0  0.000 3.033  0  0
#> x1 2.988  0  2.002 0.000  0  0
#> x2 0.000  0  0.000 5.993  0  0
#> x3 0.000  0  0.000 0.000  0  0
#> x4 8.000  0 -1.000 0.000  0  0
#> x5 4.015  0  0.000 0.000  0  0

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 ℓ2\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 

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)
#>       x0 x1     x2     x3     x4    x5
#> x0 0.000  0 -0.040  3.274  0.000 0.000
#> x1 3.237  0  1.965  0.014 -0.034 0.006
#> x2 0.000  0  0.000  5.993  0.000 0.000
#> x3 0.000  0  0.000  0.000  0.000 0.000
#> x4 7.992  0 -1.062  0.394  0.000 0.000
#> x5 3.873  0  0.069 -0.315  0.018 0.000
round(fit_lasso$adjacency_matrix,  3)
#>       x0 x1     x2    x3 x4 x5
#> x0 0.000  0  0.000 3.030  0  0
#> x1 2.939  0  1.965 0.185  0  0
#> x2 0.000  0  0.000 5.993  0  0
#> x3 0.000  0  0.000 0.000  0  0
#> x4 7.924  0 -0.960 0.000  0  0
#> x5 3.975  0  0.000 0.000  0  0
round(fit_alasso$adjacency_matrix, 3)
#>       x0 x1     x2    x3 x4 x5
#> x0 0.000  0  0.000 3.033  0  0
#> x1 2.988  0  2.002 0.000  0  0
#> x2 0.000  0  0.000 5.993  0  0
#> x3 0.000  0  0.000 0.000  0  0
#> x4 8.000  0 -1.000 0.000  0  0
#> x5 4.015  0  0.000 0.000  0  0
round(fit_ridge$adjacency_matrix,  3)
#>       x0 x1     x2     x3    x4    x5
#> x0 0.000  0 -0.017  3.132 0.000 0.000
#> x1 1.863  0  1.987  0.656 0.071 0.127
#> x2 0.000  0  0.000  5.993 0.000 0.000
#> x3 0.000  0  0.000  0.000 0.000 0.000
#> x4 7.927  0 -0.997  0.203 0.000 0.000
#> x5 2.407  0  0.254 -0.251 0.197 0.000

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) - λ=5/n1.75\lambda = 5 / n^{1.75}. Guarantees the theoretical oracle property 
# 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)
#> [1] 7
sum(fit_lam_min$adjacency_matrix != 0)
#> [1] 7

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.

result <- x1k$data |>
  lingam_direct()

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

round(p_vals, 3)
#>       x0    x1    x2    x3    x4    x5
#> x0    NA 0.988 0.214 0.976 0.876 0.952
#> x1 0.988    NA 0.986 0.991 0.328 0.882
#> x2 0.214 0.986    NA 0.919 0.051 0.124
#> x3 0.976 0.991 0.919    NA 0.934 0.978
#> x4 0.876 0.328 0.051 0.934    NA 0.650
#> x5 0.952 0.882 0.124 0.978 0.650    NA

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).

set.seed(0)
truth <- generate_lingam_sample_6(noise_dist = "uniform")

truth$true_adjacency |>
  round(1)
#>    x0 x1 x2 x3 x4 x5
#> x0  0  0  0  3  0  0
#> x1  3  0  2  0  0  0
#> x2  0  0  0  6  0  0
#> x3  0  0  0  0  0  0
#> x4  8  0 -1  0  0  0
#> x5  4  0  0  0  0  0

The causal graph of the true structure:

truth$true_adjacency |>
  plot_adjacency(
    labels = colnames(truth$data),
    title  = "True structure"
  )

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

fit_uniform <- lingam_direct(truth$data)

# Estimated causal order (the true root x3 comes first)
colnames(truth$data)[fit_uniform$causal_order]
#> [1] "x3" "x2" "x0" "x4" "x5" "x1"

# The estimated adjacency matrix recovers the true structure almost perfectly
fit_uniform$adjacency_matrix |>
  round(1)
#>    x0 x1 x2 x3 x4 x5
#> x0  0  0  0  3  0  0
#> x1  3  0  2  0  0  0
#> x2  0  0  0  6  0  0
#> x3  0  0  0  0  0  0
#> x4  8  0 -1  0  0  0
#> x5  4  0  0  0  0  0

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

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.

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]
#> [1] "x1" "x2" "x5" "x3" "x4" "x0"

fit_gauss$adjacency_matrix |>
  round(1)
#>    x0  x1   x2 x3  x4  x5
#> x0  0 0.1  0.0  0 0.1 0.0
#> x1  0 0.0  0.0  0 0.0 0.0
#> x2  0 0.3  0.0  0 0.0 0.0
#> x3  0 0.0  0.2  0 0.0 0.0
#> x4  0 0.9 -2.6  0 0.0 1.3
#> x5  0 1.2 -2.1  0 0.0 0.0

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

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.

# Shapiro-Wilk (default)
x1k$data |>
  test_residual_normality(result)
#> === Residual Normality Test ===
#> Method:         shapiro
#> Sample size:    1000
#> Significance:   0.050
#> Non-Gaussian:   6 / 6 variables
#> 
#>  variable statistic   p_value is_non_gauss skewness kurtosis
#>        x0    0.9516 < 2.2e-16         TRUE    0.061   -1.215
#>        x1    0.9521 < 2.2e-16         TRUE    0.026   -1.213
#>        x2    0.9557 < 2.2e-16         TRUE    0.083   -1.170
#>        x3    0.9578  2.25e-16         TRUE    0.025   -1.163
#>        x4    0.9544 < 2.2e-16         TRUE   -0.003   -1.206
#>        x5    0.9536 < 2.2e-16         TRUE   -0.052   -1.206
#> 
#> Interpretation:
#>   is_non_gauss = TRUE  -> rejects normality (supports LiNGAM assumption)
#>   is_non_gauss = FALSE -> cannot reject normality (LiNGAM may not fit)
#> 
#> All residuals are non-Gaussian. LiNGAM assumption is supported.

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

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.

x1k$data |>
  summary_lingam(result)
#> === Direct LiNGAM Model Summary ===
#> Variables:    6
#> Observations: 1000
#> Edges:        7
#> Causal order: x3 -> x2 -> x0 -> x4 -> x5 -> x1
#> 
#> --- Assumption 1: Independence of residuals ---
#> Method:           spearman
#> Dependent pairs:  0 / 15  (p < 0.050)
#> Min p-value:      0.0510
#> => Residuals appear mutually independent (assumption supported).
#> 
#> --- Assumption 2: Non-Gaussianity of residuals ---
#> Method:           shapiro
#> Non-Gaussian:     6 / 6  (p <= 0.050)
#> => All residuals are non-Gaussian (assumption supported).

Bootstrap Direct LiNGAM

We assess the reliability of the model using the bootstrap method.

bs_model <- x1k$data |>
  lingam_direct_bootstrap(n_sampling = 100L, seed = 42)
#> Bootstrap: 100 iterations, method=adaptive_lasso (sequential)
#>   iteration 1 / 100
#>   iteration 10 / 100
#>   iteration 20 / 100
#>   iteration 30 / 100
#>   iteration 40 / 100
#>   iteration 50 / 100
#>   iteration 60 / 100
#>   iteration 70 / 100
#>   iteration 80 / 100
#>   iteration 90 / 100
#>   iteration 100 / 100
#> Completed in 3.6 seconds.

bs_model
#> BootstrapResult: 100 samplings, 6 features

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).

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.

bs_model |>
  get_causal_direction_counts(labels = names(x1k$data))
#>    from to count proportion mean_effect median_effect  sd_effect    ci_lower
#> 1     1  6   100       1.00  4.01532920    4.01513886 0.01126767  3.99550980
#> 2     1  2    99       0.99  2.98181621    2.97864538 0.02849338  2.92980702
#> 3     1  5    99       0.99  8.00994011    8.00748238 0.02951185  7.95680521
#> 4     3  2    99       0.99  2.00498455    2.00660933 0.01479861  1.97675886
#> 5     3  5    99       0.99 -1.00529230   -1.00485827 0.01523485 -1.03801290
#> 6     4  1    99       0.99  3.03521019    3.03586526 0.03001961  2.97855949
#> 7     4  3    99       0.99  5.99644109    5.99745219 0.03186571  5.94050363
#> 8     2  1     1       0.01  0.05304916    0.05304916 0.00000000  0.05304916
#> 9     2  3     1       0.01  0.40196452    0.40196452 0.00000000  0.40196452
#> 10    2  5     1       0.01  0.90679690    0.90679690 0.00000000  0.90679690
#> 11    3  4     1       0.01  0.16166764    0.16166764 0.00000000  0.16166764
#> 12    5  1     1       0.01  0.10453910    0.10453910 0.00000000  0.10453910
#> 13    5  3     1       0.01 -0.13636255   -0.13636255 0.00000000 -0.13636255
#>       ci_upper from_name to_name
#> 1   4.03698551        x0      x5
#> 2   3.03860193        x0      x1
#> 3   8.07414013        x0      x4
#> 4   2.03193816        x2      x1
#> 5  -0.97488336        x2      x4
#> 6   3.09304642        x3      x0
#> 7   6.06134091        x3      x2
#> 8   0.05304916        x1      x0
#> 9   0.40196452        x1      x2
#> 10  0.90679690        x1      x4
#> 11  0.16166764        x2      x3
#> 12  0.10453910        x4      x0
#> 13 -0.13636255        x4      x2

Adjacency Matrix of Mean Causal Effects

We construct an adjacency matrix from the bootstrap results.

bs_adjacency_matrix <- bs_model |>
  get_adjacency_matrix_summary(stat = "median")

bs_adjacency_matrix |>
  round(3)
#>       [,1]  [,2]   [,3]  [,4]   [,5] [,6]
#> [1,] 0.000 0.053  0.000 3.036  0.105    0
#> [2,] 2.979 0.000  2.007 0.000  0.000    0
#> [3,] 0.000 0.402  0.000 5.997 -0.136    0
#> [4,] 0.000 0.000  0.162 0.000  0.000    0
#> [5,] 8.007 0.907 -1.005 0.000  0.000    0
#> [6,] 4.015 0.000  0.000 0.000  0.000    0

We visualize the estimated adjacency matrix.

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.

bs_model |>
  get_probabilities()
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.00 0.01 0.00 0.99 0.01    0
#> [2,] 0.99 0.00 0.99 0.00 0.00    0
#> [3,] 0.00 0.01 0.00 0.99 0.01    0
#> [4,] 0.00 0.00 0.01 0.00 0.00    0
#> [5,] 0.99 0.01 0.99 0.00 0.00    0
#> [6,] 1.00 0.00 0.00 0.00 0.00    0

Mean Total Effects

We compute the mean total effect of each path.

bs_model |>
  get_total_causal_effects()
#>    from to      effect probability
#> 1     1  6  4.01520158        1.00
#> 2     1  2  2.87431611        0.99
#> 3     1  5  7.90813117        0.99
#> 4     3  2  1.95874622        0.99
#> 5     3  5 -1.06193484        0.99
#> 6     4  1  3.03586526        0.99
#> 7     4  2 21.07027271        0.99
#> 8     4  3  5.99805118        0.99
#> 9     4  5 18.28272145        0.99
#> 10    4  6 12.18719857        0.99
#> 11    3  6 -0.24574320        0.04
#> 12    2  1  0.14794503        0.01
#> 13    2  3  0.27850920        0.01
#> 14    2  4  0.04611007        0.01
#> 15    2  5  0.90679690        0.01
#> 16    2  6  0.59359217        0.01
#> 17    3  4  0.16192779        0.01
#> 18    5  1  0.10498716        0.01
#> 19    5  3 -0.13625059        0.01
#> 20    5  6  0.42156703        0.01
#> 21    6  2  0.24518629        0.01

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

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).

bs_model |>
  get_causal_order_stability(labels = names(x1k$data))
#> === Causal Order Stability ===
#> Bootstrap samples:       100
#> Overall stability score: 0.736  (0 = random, 1 = fully stable)
#> 
#> Rank summary (sorted by mean rank; 1 = most upstream):
#>  variable mean_rank sd_rank median_rank mode_rank
#>        x3      1.05    0.50           1         1
#>        x0      2.62    0.51           3         3
#>        x2      2.75    0.95           2         2
#>        x5      4.41    1.23           4         3
#>        x4      4.92    0.77           5         5
#>        x1      5.25    0.88           5         6
#> 
#> Precedence probability P[i, j] = P(variable i precedes j):
#>      x0   x1   x2   x3   x4   x5
#> x0 0.00 0.99 0.39 0.01 0.99 1.00
#> x1 0.01 0.00 0.01 0.01 0.38 0.34
#> x2 0.61 0.99 0.00 0.01 0.99 0.65
#> x3 0.99 0.99 0.99 0.00 0.99 0.99
#> x4 0.01 0.62 0.01 0.01 0.00 0.43
#> x5 0.00 0.66 0.35 0.01 0.57 0.00

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.

# Convert the estimated adjacency matrix to an edge list
tidy(model)
#>   from to  estimate
#> 1   x0 x1  2.987705
#> 2   x0 x4  8.000096
#> 3   x0 x5  4.014962
#> 4   x2 x1  2.001708
#> 5   x2 x4 -1.000306
#> 6   x3 x0  3.032952
#> 7   x3 x2  5.992677

# One-row summary of the whole model
glance(model)
#>   n_variables n_edges                     causal_order
#> 1           6       7 x3 -> x2 -> x0 -> x4 -> x5 -> x1

# Direction-wise summary of the bootstrap results (variable names via labels)
tidy(bs_model, labels = names(x1k$data))
#>    from to count proportion mean_effect median_effect  sd_effect    ci_lower
#> 1     1  6   100       1.00  4.01532920    4.01513886 0.01126767  3.99550980
#> 2     1  2    99       0.99  2.98181621    2.97864538 0.02849338  2.92980702
#> 3     1  5    99       0.99  8.00994011    8.00748238 0.02951185  7.95680521
#> 4     3  2    99       0.99  2.00498455    2.00660933 0.01479861  1.97675886
#> 5     3  5    99       0.99 -1.00529230   -1.00485827 0.01523485 -1.03801290
#> 6     4  1    99       0.99  3.03521019    3.03586526 0.03001961  2.97855949
#> 7     4  3    99       0.99  5.99644109    5.99745219 0.03186571  5.94050363
#> 8     2  1     1       0.01  0.05304916    0.05304916 0.00000000  0.05304916
#> 9     2  3     1       0.01  0.40196452    0.40196452 0.00000000  0.40196452
#> 10    2  5     1       0.01  0.90679690    0.90679690 0.00000000  0.90679690
#> 11    3  4     1       0.01  0.16166764    0.16166764 0.00000000  0.16166764
#> 12    5  1     1       0.01  0.10453910    0.10453910 0.00000000  0.10453910
#> 13    5  3     1       0.01 -0.13636255   -0.13636255 0.00000000 -0.13636255
#>       ci_upper from_name to_name
#> 1   4.03698551        x0      x5
#> 2   3.03860193        x0      x1
#> 3   8.07414013        x0      x4
#> 4   2.03193816        x2      x1
#> 5  -0.97488336        x2      x4
#> 6   3.09304642        x3      x0
#> 7   6.06134091        x3      x2
#> 8   0.05304916        x1      x0
#> 9   0.40196452        x1      x2
#> 10  0.90679690        x1      x4
#> 11  0.16166764        x2      x3
#> 12  0.10453910        x4      x0
#> 13 -0.13636255        x4      x2

A Larger Dataset (10 Variables)

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

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=1000n = 1000) with both methods.

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"]))
#> Direct LiNGAM : 0.01 sec
#> ICA-LiNGAM    : 0.02 sec

Comparing the Estimated Coefficients

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

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)
#>   from to direct    ica
#> 1   x0 x1  2.988  3.245
#> 2   x0 x4  8.000  7.999
#> 3   x0 x5  4.015  3.876
#> 4   x2 x1  2.002  1.973
#> 5   x2 x4 -1.000 -1.060
#> 6   x3 x0  3.033  3.027
#> 7   x3 x2  5.993  6.101

Comparing the DAG Structures

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

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
#>   from to direct    ica truth
#> 1   x0 x1  2.988  3.245  TRUE
#> 2   x0 x4  8.000  7.999  TRUE
#> 3   x0 x5  4.015  3.876  TRUE
#> 4   x2 x1  2.002  1.973  TRUE
#> 5   x2 x4 -1.000 -1.060  TRUE
#> 6   x3 x0  3.033  3.027  TRUE
#> 7   x3 x2  5.993  6.101  TRUE

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 pp and the number of tests per step is at most p(p−1)p(p-1), the total number of independence tests is approximately

∑k=1pk(k−1)≈p33\sum_{k=1}^{p} k(k-1) \approx \frac{p^3}{3}

giving a computational cost of O(p3)O(p^3). By contrast, the FastICA used by ICA-LiNGAM is O(p2n)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 xix_i (i≥1i \ge 1) randomly has at most max_parents parents chosen from x0,…,xi−1x_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 

d20 <- generate_lingam_large_sample(p = 20, n = 1000, seed = 42)

dim(d20$data)                    # 1000 rows x 20 columns
#> [1] 1000   20
sum(d20$true_adjacency != 0)     # number of true edges (sparse DAG)
#> [1] 32
d20$true_causal_order            # 0, 1, ..., 19
#>  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19

Comparing Execution Times

When pp grows by a factor of 1.5 (10 -> 15), the number of independence tests grows by a factor of 153/103≈3.415^3 / 10^3 \approx 3.4.

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)
))
#> p = 10 : 0.03 sec
#> p = 15 : 0.06 sec
#> theoretical factor 3.4x vs. observed 2.1x

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

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"]
))
#>               p = 10   p = 15
#> Direct LiNGAM :  0.03 sec   0.06 sec
#> ICA-LiNGAM    :  0.02 sec   0.03 sec

The larger pp becomes, the more Direct LiNGAM’s O(p3)O(p^3) cost dominates, and the gap between the two widens. In large-scale settings such as p=30p = 30 or p=50p = 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.

# Estimated causal order
r10$causal_order
#>  [1]  1  2  3  7  4  5  9  8  6 10

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

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

tidy(r10) |>
  head(10)
#>    from to   estimate
#> 1    x0 x1 -1.3787175
#> 2    x0 x2  1.1069608
#> 3    x0 x3  0.9365537
#> 4    x0 x5  1.2879537
#> 5    x1 x2  0.9099343
#> 6    x1 x3  1.4225647
#> 7    x1 x5 -1.2930266
#> 8    x1 x6  1.4634025
#> 9    x1 x9  1.2511988
#> 10   x2 x3 -1.4992033

High-Dimensional Direct LiNGAM

The O(p3)O(p^3) independence-test cost shown above becomes a real bottleneck once pp grows into the tens or hundreds, and breaks down entirely once p>np > 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.

hd_sample <- generate_lingam_sample_6(n = 500, seed = 1)
hd_result <- lingam_high_dim(hd_sample$data)

hd_result$causal_order
#> [1] 4 3 1 5 2 6
round(hd_result$adjacency_matrix, 3)
#>       x0 x1     x2    x3 x4 x5
#> x0 0.000  0  0.000 2.968  0  0
#> x1 2.970  0  2.013 0.000  0  0
#> x2 0.000  0  0.000 6.010  0  0
#> x3 0.000  0  0.000 0.000  0  0
#> x4 8.023  0 -1.000 0.000  0  0
#> x5 4.013  0  0.000 0.000  0  0

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:

wide_sample <- generate_lingam_large_sample(p = 30, n = 25, seed = 1)
wide_result <- lingam_high_dim(wide_sample$data)
#> Warning: Since n_samples <= n_features, the adjacency matrix is estimated with
#> cross-validated lasso (cv.glmnet) instead of BIC-based lambda selection.
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold
#> Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
#> fold

wide_result$causal_order
#>  [1]  2  1  3  8 20 13 27  4  5 12 26  6  7 16  9 15 14 23 10 18 21 11 22 24 29
#> [26] 30 25 19 17 28

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).
paradox <- generate_lingam_paradox_data(n = 2000L, seed = 42)

head(paradox$data)
#>             x0         x1          x2         x3
#> 1  0.780627610  2.0872183  1.95046049  1.1209218
#> 2  0.529343129  1.1562639  1.86870201  1.6129261
#> 3 -1.193165251 -0.2515850 -0.43614264 -0.9056694
#> 4 -0.056001104  1.6615506  2.07542227  0.7890187
#> 5  0.004312424  1.0175487 -0.02532253 -0.3155891
#> 6  0.658064158  0.4833892  0.25385608  0.0167021

# All variables are standardized (sd = 1)
sapply(paradox$data, sd)
#> x0 x1 x2 x3 
#>  1  1  1  1

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

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.

model_p <- lingam_direct(paradox$data)

# Estimated causal order
colnames(paradox$data)[model_p$causal_order]
#> [1] "x1" "x2" "x0" "x3"

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.

model_p$adjacency_matrix |>
  round(3)
#>    x0    x1    x2 x3
#> x0  0 0.558 0.000  0
#> x1  0 0.000 0.000  0
#> x2  0 0.833 0.000  0
#> x3  0 0.000 0.822  0

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.

bs_paradox <- paradox$data |>
  lingam_direct_bootstrap(n_sampling = 100L, seed = 42)
#> Bootstrap: 100 iterations, method=adaptive_lasso (sequential)
#>   iteration 1 / 100
#>   iteration 10 / 100
#>   iteration 20 / 100
#>   iteration 30 / 100
#>   iteration 40 / 100
#>   iteration 50 / 100
#>   iteration 60 / 100
#>   iteration 70 / 100
#>   iteration 80 / 100
#>   iteration 90 / 100
#>   iteration 100 / 100
#> Completed in 1.5 seconds.

# Occurrence probability of each direction (row = to, column = from)
bs_paradox |>
  get_probabilities() |>
  round(2)
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    1 0.05 0.01
#> [2,]    0    0 0.00 0.00
#> [3,]    0    1 0.00 0.00
#> [4,]    0    0 1.00 0.00

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 B0B_0. The model is:

Xt=B0Xt+∑k=1pBkXt−k+etX_t = B_0\,X_t + \sum_{k=1}^{p} B_k\,X_{t-k} + e_t

where B0B_0 encodes contemporaneous causal effects (strictly acyclic), B1,…,BpB_1, \ldots, B_p encode lagged effects, and ete_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 x0→x1→x2x_0 \to x_1 \to x_2 (coefficients 0.6 and −0.5), and the only cross-variable lag-1 effect is x2(t−1)→x0(t)x_2(t-1) \to x_0(t) (coefficient 0.3).

s <- generate_varlingam_sample(n = 1000, seed = 42)

# True instantaneous coefficient matrix B0  (B0[i, j]: x_j -> x_i)
s$true_B0
#>      [,1] [,2] [,3]
#> [1,]  0.0  0.0    0
#> [2,]  0.6  0.0    0
#> [3,]  0.0 -0.5    0

# True lag-1 coefficient matrix  (M1[i, j]: x_j(t-1) -> x_i(t), structural)
s$true_M1
#>      [,1] [,2] [,3]
#> [1,]  0.4  0.0  0.3
#> [2,]  0.0  0.3  0.0
#> [3,]  0.0  0.0  0.5

Fitting VAR-LiNGAM

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

model <- lingam_var(s$data, lags = 1)
model
#> VAR-LiNGAM Result
#>   Variables : 3
#>   Lag order : 1
#>   Causal order (instantaneous): x0 -> x1 -> x2
#> 
#> Instantaneous adjacency matrix B0 (row = to, col = from):
#>       x0     x1 x2
#> x0 0.000  0.000  0
#> x1 0.576  0.000  0
#> x2 0.000 -0.491  0
#> 
#> Lagged adjacency matrix B1 (row = to, col = from):
#>     x0    x1    x2
#> x0 0.4 0.000 0.309
#> x1 0.0 0.225 0.000
#> x2 0.0 0.000 0.495

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

  • [1, , ] ("lag0"): instantaneous matrix B0B_0. B0[i, j] is the direct effect of xjx_j on xix_i at the same time step.
  • [k + 1, , ] ("lagk"): lagged matrix BkB_k. Bk[i, j] is the direct structural effect of xj(t−k)x_j(t-k) on xi(t)x_i(t).

Both can be extracted by their dimension label:

B0 <- model$adjacency_matrices["lag0", , ]
B1 <- model$adjacency_matrices["lag1", , ]

round(B0, 2)  # compare with s$true_B0
#>      x0    x1 x2
#> x0 0.00  0.00  0
#> x1 0.58  0.00  0
#> x2 0.00 -0.49  0
round(B1, 2)  # compare with s$true_M1
#>     x0   x1   x2
#> x0 0.4 0.00 0.31
#> x1 0.0 0.23 0.00
#> x2 0.0 0.00 0.50

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:

# 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.

check_var_stationarity(model)
#> === VAR Stationarity Check ===
#> Lag order:         1
#> Max |eigenvalue|:  0.4942  (threshold 1.00)
#> Stationary:        YES

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 ete_t are non-Gaussian. test_varlingam_residual_normality() tests whether the LiNGAM innovations et=(I−B0)nte_t = (I - B_0)\,n_t (where ntn_t are the stored VAR residuals) depart from normality. A small p-value (reject H0H_0: Gaussian) supports the model assumption.

test_varlingam_residual_normality(model)
#> === Residual Normality Test ===
#> Method:         shapiro
#> Sample size:    999
#> Significance:   0.050
#> Non-Gaussian:   3 / 3 variables
#> 
#>  variable statistic   p_value is_non_gauss skewness kurtosis
#>        x0    0.9498 < 2.2e-16         TRUE    0.088   -1.220
#>        x1    0.9536 < 2.2e-16         TRUE   -0.007   -1.238
#>        x2    0.9518 < 2.2e-16         TRUE   -0.046   -1.221
#> 
#> Interpretation:
#>   is_non_gauss = TRUE  -> rejects normality (supports LiNGAM assumption)
#>   is_non_gauss = FALSE -> cannot reject normality (LiNGAM may not fit)
#> 
#> All residuals are non-Gaussian. LiNGAM assumption is supported.

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

test_varlingam_residual_normality_all(model, methods = c("shapiro", "jb"))
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
#>   variable     skewness  kurtosis    p_shapiro         p_jb all_non_gauss
#> 1       x0  0.088013433 -1.219504 6.041150e-18 1.898481e-14          TRUE
#> 2       x1 -0.007060832 -1.238431 3.114110e-17 1.365574e-14          TRUE
#> 3       x2 -0.046381574 -1.220525 1.452023e-17 2.864375e-14          TRUE

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

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 xj(t−1)x_j(t-1) on xi(t)x_i(t).

# Total effect x0 -> x2 (contemporaneous)
estimate_var_total_effect(s$data, model, from_index = 1, to_index = 3)
#> [1] -0.2582049

# 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)
#> [1] -0.2752551

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.

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.

round(get_var_probabilities(bs_var), 2)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.00    0    0 1.00    0 1.00
#> [2,] 1.00    0    0 0.01    1 0.01
#> [3,] 0.02    1    0 0.00    0 1.00

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.

# Paths from x0 to x2 at the same time step (from_lag = 0)
get_var_paths(bs_var, from_index = 1, to_index = 3)
#>      path      effect probability
#> 1 1, 2, 3 -0.28289316        1.00
#> 2    1, 3  0.08085689        0.02
# 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)
#>            path       effect probability
#> 1    4, 1, 2, 3 -0.112845948        1.00
#> 2    4, 5, 2, 3 -0.062338986        1.00
#> 3  4, 5, 6,....  0.024970830        1.00
#> 4    4, 5, 6, 3 -0.139242633        1.00
#> 5       4, 1, 3  0.034930592        0.02
#> 6  4, 5, 6,.... -0.007793044        0.02
#> 7  4, 6, 1,.... -0.007793044        0.02
#> 8    4, 6, 1, 3  0.002044126        0.02
#> 9       4, 6, 3  0.040411503        0.02
#> 10      4, 2, 3 -0.068176510        0.01
#> 11 4, 5, 6,.... -0.018273144        0.01

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).

set.seed(1)
lim_dat <- generate_lim_sample(n = 2000)
head(lim_dat$data)
#>           x1 x2         x3
#> 1  0.1182559  0 -0.8936636
#> 2 -1.8695490  0 -1.2651618
#> 3 -3.3867259  0  2.1530815
#> 4 -0.4395899  1  0.6618645
#> 5  0.3215812  0 -0.4652433
#> 6  1.6721779  0  0.8429619
lim_dat$is_continuous
#> [1]  TRUE FALSE  TRUE

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().

lim_result <- lingam_lim(lim_dat$data, is_continuous = lim_dat$is_continuous)
print(lim_result)
#> LiM Result
#>   Variables : 3
#>   Variable types: continuous, discrete, continuous
#>   Causal order: x1 -> x2 -> x3
#> 
#> Adjacency matrix (row = to, col = from):
#>    x1    x2 x3
#> x1  0 0.000  0
#> x2  1 0.000  0
#> x3  0 1.657  0

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.

colnames(lim_dat$data)[lim_result$causal_order]
#> [1] "x1" "x2" "x3"

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.

mg <- generate_multi_group_sample(n = c(1000, 1000), seed = 42)
lapply(mg$data_list, head, 3)
#> $group1
#>         x0        x1       x2        x3        x4        x5
#> 1 2.814924 18.017120 4.543655 0.6333728 18.160090 12.236660
#> 2 1.889685 10.956005 2.188091 0.3175366 13.172754  7.932657
#> 3 1.008905  6.990652 1.953131 0.2409218  6.702107  4.797122
#> 
#> $group2
#>          x0        x1        x2         x3        x4        x5
#> 1 0.7259014  5.482225 0.9301592 0.01259095  5.275903  3.321061
#> 2 2.4321051 17.252303 3.3989705 0.41696287 16.459975 11.368701
#> 3 1.5550457 10.342355 1.8713591 0.24518297 10.520165  7.859683
mg_result <- lingam_multi_group(mg$data_list, reg_method = "ols")
print(mg_result)
#> Multi-Group Direct LiNGAM Result
#>   Groups      : 2 (group1, group2)
#>   Variables   : 6
#>   Causal order (common): x3 -> x0 -> x5 -> x2 -> x4 -> x1
#> 
#> [group1] Adjacency matrix (row = to, col = from):
#>        x0 x1     x2     x3     x4    x5
#> x0  0.000  0  0.000  3.033  0.000 0.000
#> x1  3.237  0  1.965  0.014 -0.034 0.006
#> x2 -0.236  0  0.000  6.112  0.000 0.049
#> x3  0.000  0  0.000  0.000  0.000 0.000
#> x4  7.921  0 -1.063  0.399  0.000 0.018
#> x5  4.016  0  0.000 -0.003  0.000 0.000
#> 
#> [group2] Adjacency matrix (row = to, col = from):
#>       x0 x1     x2     x3    x4     x5
#> x0 0.000  0  0.000  3.504 0.000  0.000
#> x1 2.732  0  2.568  0.083 0.034  0.093
#> x2 0.154  0  0.000  6.322 0.000 -0.024
#> x3 0.000  0  0.000  0.000 0.000  0.000
#> x4 8.483  0 -1.487 -0.110 0.000  0.006
#> x5 4.515  0  0.000 -0.045 0.000  0.000

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():

g1 <- get_group_result(mg_result, "group1")
class(g1)
#> [1] "LingamResult"

estimate_all_total_effects(mg$data_list$group1, g1, method = "ols")
#>             x0 x1        x2        x3          x4          x5
#> x0  0.00000000  0  0.000000  3.033460  0.00000000  0.00000000
#> x1  2.90911952  0  2.001580 21.058733 -0.03397056  0.10299386
#> x2 -0.03933572  0  0.000000  5.992677  0.00000000  0.04894766
#> x3  0.00000000  0  0.000000  0.000000  0.00000000  0.00000000
#> x4  8.03407606  0 -1.062516 18.276121  0.00000000 -0.03416285
#> x5  4.01586857  0  0.000000 12.179395  0.00000000  0.00000000

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:

mg_bs <- lingam_multi_group_bootstrap(mg$data_list,
  n_sampling = 20L, reg_method = "ols", seed = 1, verbose = FALSE
)
get_probabilities(mg_bs$group1)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]  0.0  0.0 0.30    1 0.00 0.00
#> [2,]  1.0  0.0 1.00    1 0.70 0.80
#> [3,]  0.7  0.0 0.00    1 0.00 0.45
#> [4,]  0.0  0.0 0.00    0 0.00 0.00
#> [5,]  1.0  0.3 1.00    1 0.00 0.55
#> [6,]  1.0  0.2 0.55    1 0.45 0.00

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().

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
bwi <- bootstrap_with_imputation(X_na,
  n_sampling = 20L, n_repeats = 5L, seed = 42, verbose = FALSE
)
print(bwi)
#> ImputationBootstrapResult: 20 samplings x 5 repeats, 6 features, 100 missing cells (original data)

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:

bs_na <- as_bootstrap_result(bwi, aggregate = "median")
get_probabilities(bs_na)
#>    x0 x1 x2 x3 x4 x5
#> x0  0  0  0  1  0  0
#> x1  1  0  1  0  0  0
#> x2  0  0  0  1  0  0
#> x3  0  0  0  0  0  0
#> x4  1  0  1  0  0  0
#> x5  1  0  0  0  0  0

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.

# 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)
#>          x0        x1        x2        x3          x4        x5
#> 1 0.6154746 1.7104554 1.0618261 1.0851944  0.57917376 0.8337563
#> 2 1.5905703 2.4770365 1.4291087 1.4325230  0.56022597 0.8686206
#> 3 1.1007549 2.1817888 1.5289901 1.8037643 -0.03671602 1.4001192
#> 4 1.7744689 2.7106515 2.7714036 2.4797583 -0.10133399 1.3102925
#> 5 0.5433612 0.7244786 0.5217204 0.8755981  0.82504738 1.2597767
#> 6 1.7671488 1.8838085 1.8358794 2.7663075  0.90699176 1.3624485
confounded$confounded_pair
#> [1] 3 4
parce_result <- lingam_parce(confounded$data, reg_method = "ols")
print(parce_result)
#> Bottom-Up ParceLiNGAM Result
#>   Variables : 6
#>   Independence measure: hsic
#>   Causal order: (x2, x3) -> x0 -> x4 -> x5 -> x1
#>   (NA entries in the adjacency matrix = unresolved order / suspected latent confounding)
#> 
#> Adjacency matrix (row = to, col = from):
#>       x0 x1     x2     x3    x4     x5
#> x0 0.000  0 -0.010  0.516 0.000  0.000
#> x1 0.479  0  0.447  0.060 0.025 -0.049
#> x2 0.000  0  0.000     NA 0.000  0.000
#> x3 0.000  0     NA  0.000 0.000  0.000
#> x4 0.497  0 -0.490 -0.001 0.000  0.000
#> x5 0.436  0  0.068  0.023 0.050  0.000

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:

parce_result$causal_order[[1]]
#> [1] 3 4
parce_result$adjacency_matrix[confounded$confounded_pair, confounded$confounded_pair]
#>    x2 x3
#> x2  0 NA
#> x3 NA  0

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:

# 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"
)
#> Warning in estimate_total_effect_parce(confounded$data, parce_result,
#> from_index = confounded$confounded_pair[1], : x2 is part of an unresolved
#> causal order (suspected latent confounding); total effect cannot be estimated.
#> [1] NA

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

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.

parce_bs <- lingam_parce_bootstrap(confounded$data,
  n_sampling = 10L, reg_method = "ols", seed = 1, verbose = FALSE
)
get_probabilities(parce_bs)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]  0.0  0.2  0.5  0.4  0.0  0.0
#> [2,]  0.5  0.0  0.5  0.4  0.2  0.3
#> [3,]  0.0  0.0  0.0  0.0  0.0  0.0
#> [4,]  0.1  0.2  0.2  0.0  0.1  0.0
#> [5,]  0.7  0.5  0.7  0.6  0.0  0.3
#> [6,]  0.6  0.4  0.6  0.6  0.4  0.0

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.

# 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)
#>            x0           x1         x2           x3         x4            x5
#> 1 -0.96839700 -0.023398020 -0.7514703 -0.473147026 -0.9674898 -0.0307310344
#> 2  1.31480945  0.424388510  1.0389688  0.001304296  1.4111680  0.0007741684
#> 3 -0.85973904 -0.025330472 -1.1731736 -0.034157640 -1.3389435 -0.0729373397
#> 4 -0.76463976  0.324413350 -0.7116352  0.078719765 -0.3333176  0.5074829194
#> 5  0.08152383  0.002364107 -0.2546673  0.036715357 -0.2040856  0.0044720536
#> 6 -0.30389184 -0.225029302 -0.5032249 -0.172267536  0.5221642 -0.0690391705
rcd_confounded$confounded_pair
#> [1] 3 5
rcd_result <- lingam_rcd(rcd_confounded$data)
print(rcd_result)
#> RCD Result
#>   Variables : 6
#> 
#> Ancestor sets:
#>   M(x0) = {x1, x3, x5}
#>   M(x1) = {x5}
#>   M(x2) = {x0, x1, x3, x5}
#>   M(x3) = {x5}
#>   M(x4) = {x0, x1, x3, x5}
#>   M(x5) = {}
#> 
#>   (NA entries in the adjacency matrix = suspected shared latent confounder)
#> 
#> Adjacency matrix (row = to, col = from):
#>       x0    x1 x2    x3 x4    x5
#> x0 0.000 1.116  0 0.989  0 0.000
#> x1 0.000 0.000  0 0.000  0 0.588
#> x2 0.810 0.000  0 0.000 NA 0.000
#> x3 0.000 0.000  0 0.000  0 0.449
#> x4 1.015 0.000 NA 0.000  0 0.000
#> x5 0.000 0.000  0 0.000  0 0.000

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

rcd_result$ancestors_list
#> $x0
#> [1] 2 4 6
#> 
#> $x1
#> [1] 6
#> 
#> $x2
#> [1] 1 2 4 6
#> 
#> $x3
#> [1] 6
#> 
#> $x4
#> [1] 1 2 4 6
#> 
#> $x5
#> integer(0)
rcd_result$adjacency_matrix[rcd_confounded$confounded_pair, rcd_confounded$confounded_pair]
#>    x2 x4
#> x2  0 NA
#> x4 NA  0

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

# 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]
)
#> Warning in estimate_total_effect_rcd(rcd_confounded$data, rcd_result,
#> from_index = rcd_confounded$confounded_pair[1], : x2 is part of a suspected
#> latent confounder pair; total effect cannot be estimated.
#> [1] NA

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

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.

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)
#>   DoF DoF Baseline chi2 chi2 p-value chi2 Baseline CFI GFI AGFI NFI TLI RMSEA
#> 1   0           15    0           NA       23023.7   1   1    1   1   1     0
#>        AIC      BIC    LogLik
#> 1 1860.598 1958.753 -910.2991

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

reversed_adjacency <- t(fit_result$adjacency_matrix)
evaluate_model_fit(reversed_adjacency, sample6$data)
#>   DoF DoF Baseline         chi2 chi2 p-value chi2 Baseline CFI GFI AGFI NFI TLI
#> 1   0           15 3.197442e-11           NA       23023.7   1   1    1   1   1
#>   RMSEA       AIC       BIC   LogLik
#> 1     0 -4264.864 -4166.708 2152.432

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 nn is extremely small relative to the number of variables pp (rule of thumb: n<10pn < 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≥10pn \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.