Package {lingamr}

lingamr: 'LiNGAM' Algorithms for Causal Discovery

Description

R implementation of 'LiNGAM' (Linear Non-Gaussian Acyclic Model) algorithms for causal discovery, following Shimizu et al. (2011) https://www.jmlr.org/papers/v12/shimizu11a.html. Based on the 'Python' implementation by Ikeuchi et al. (2023) https://github.com/cdt15/lingam. The 'VAR-LiNGAM' residual diagnostics are inspired by the 'VARLiNGAM' R code of Moneta et al. https://sites.google.com/site/dorisentner/publications/VARLiNGAM.

Details

tutorial: vignette("lingamr", package = "lingamr")

Author(s)

Maintainer: Osamu Morimoto galactic.supermarket@gmail.com [copyright holder]

Authors:

• Osamu Morimoto galactic.supermarket@gmail.com [copyright holder]

Other contributors:

• T. Ikeuchi [copyright holder]

• G. Haraoka [copyright holder]

• M. Ide [copyright holder]

• W. Kurebayashi [copyright holder]

• S. Shimizu [copyright holder]

References

Shimizu, S., et al. (2011). DirectLiNGAM: A direct method for learning a linear non-Gaussian structural equation model. Journal of Machine Learning Research, 12, 1225-1248.

Hyvärinen, A., Zhang, K., Shimizu, S., & Hoyer, P. O. (2010). Estimation of a structural vector autoregression model using non-Gaussianity. Journal of Machine Learning Research, 11, 1709-1731.

Moneta, A., Entner, D., Hoyer, P. O., & Coad, A. (2013). Causal inference by independent component analysis: Theory and applications. Oxford Bulletin of Economics and Statistics, 75(5), 705-730. VARLiNGAM R code: https://sites.google.com/site/dorisentner/publications/VARLiNGAM

Shimizu, S. (2012). Joint estimation of linear non-Gaussian acyclic models. Neurocomputing, 81, 104-107. (MultiGroup Direct LiNGAM)

Tashiro, T., Shimizu, S., Hyvärinen, A., & Washio, T. (2014). ParceLiNGAM: A causal ordering method robust against latent confounders. Neural Computation, 26(1), 57-83.

Maeda, T. N., & Shimizu, S. (2020). RCD: Repetitive causal discovery of linear non-Gaussian acyclic models with latent confounders. AISTATS 2020, PMLR 108, 735-745.

Wang, Y. S., & Drton, M. (2020). High-dimensional causal discovery under non-Gaussianity. Biometrika, 107(1), 41-59.

Zeng, Y., Shimizu, S., Matsui, H., & Sun, F. (2022). Causal discovery for linear mixed data. Proceedings of the First Conference on Causal Learning and Reasoning (CLeaR 2022), PMLR 177, 994-1009.

Python implementation (DirectLiNGAM, VAR-LiNGAM, MultiGroup, ParceLiNGAM, RCD, LiM, HighDim): cdt15/lingam, https://github.com/cdt15/lingam

See Also

Useful links:

• https://github.com/morimotoosamu/lingamr

• https://morimotoosamu.github.io/lingamr/

• Report bugs at https://github.com/morimotoosamu/lingamr/issues

Collapse an ImputationBootstrapResult into a BootstrapResult

Description

bootstrap_with_imputation()'s result carries an extra n_repeats dimension (one causal-structure estimate per imputed dataset per bootstrap iteration), which the existing bootstrap analysis functions (get_probabilities(), get_causal_direction_counts(), get_directed_acyclic_graph_counts(), get_causal_order_stability(), tidy()) do not expect. This collapses that dimension by aggregating the n_repeats adjacency matrices of each iteration into one, producing a regular lingam_direct_bootstrap()-style BootstrapResult.

Usage

as_bootstrap_result(x, aggregate = c("median", "mean"))

Arguments

Value

A BootstrapResult (see lingam_direct_bootstrap()) with total_effects = NULL (total effects are not computed by bootstrap_with_imputation()); calling get_total_causal_effects() on it raises the usual "not computed" error.

Examples

set.seed(1)
sample6 <- generate_lingam_sample_6(n = 300, seed = 1)
X <- sample6$data
X$x5[sample.int(nrow(X), size = 30)] <- NA

if (requireNamespace("mice", quietly = TRUE)) {
  res <- bootstrap_with_imputation(X,
    n_sampling = 5L, n_repeats = 3L, seed = 42, verbose = FALSE
  )
  bs <- as_bootstrap_result(res, aggregate = "median")
  get_probabilities(bs)

  # get_total_causal_effects() is not available: total effects were never computed
  tryCatch(get_total_causal_effects(bs), error = function(e) conditionMessage(e))
}

Plot the causal graph of a LiMResult with ggplot2

Description

Draws the estimated causal structure of a LiM model as a ggplot2-based directed graph, exactly like autoplot.LingamResult().

Usage

## S3 method for class 'LiMResult'
autoplot(
  object,
  threshold = 0,
  node_size = 16,
  node_color = "lightblue",
  label_edges = TRUE,
  label_pos = 0.35,
  ...
)

Arguments

Value

A ggplot object

Examples

if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("igraph", quietly = TRUE)) {
  library(ggplot2)
  set.seed(1)
  dat <- generate_lim_sample(n = 300)
  model <- lingam_lim(dat$data, is_continuous = dat$is_continuous)
  autoplot(model)
}

Plot the causal graph of a LingamResult with ggplot2

Description

Draws the estimated causal structure as a ggplot2-based directed graph. Node positions are computed with igraph's hierarchical layout (sugiyama), so the causal flow is generally arranged from top to bottom. Because the output is a static image, it is stable in RMarkdown / Quarto. If you need an interactive HTML figure, use plot_adjacency() (DiagrammeR-based).

Usage

## S3 method for class 'LingamResult'
autoplot(
  object,
  threshold = 0,
  node_size = 16,
  node_color = "lightblue",
  label_edges = TRUE,
  label_pos = 0.35,
  ...
)

Arguments

Details

autoplot() is a ggplot2 generic, so you must load ggplot2 with library(ggplot2) before using it. Plotting requires ggplot2 and igraph.

Value

A ggplot object

Examples

if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("igraph", quietly = TRUE)) {
  library(ggplot2)
  dat <- generate_lingam_sample_6()
  model <- lingam_direct(dat$data, reg_method = "ols")
  autoplot(model)
}

Plot one group of a MultiGroupLingamResult with ggplot2

Description

Extracts a single group's result with get_group_result() and draws it like autoplot.LingamResult(), with the group shown in the subtitle. The causal order is shared across groups, but the coefficients (and therefore the plotted edges) are group specific.

Usage

## S3 method for class 'MultiGroupLingamResult'
autoplot(
  object,
  group = 1,
  threshold = 0,
  node_size = 16,
  node_color = "lightblue",
  label_edges = TRUE,
  label_pos = 0.35,
  ...
)

Arguments

Value

A ggplot object

Examples

if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("igraph", quietly = TRUE)) {
  library(ggplot2)
  mg <- generate_multi_group_sample()
  model <- lingam_multi_group(mg$data_list, reg_method = "ols")
  autoplot(model, group = 2)
}

Plot the causal graph of a ParceLingamResult with ggplot2

Description

Draws the estimated causal structure as a ggplot2-based directed graph, like autoplot.LingamResult(). Variable pairs whose adjacency-matrix entries are NA (unresolved order / suspected latent confounding) are drawn as dashed, arrowless segments.

Usage

## S3 method for class 'ParceLingamResult'
autoplot(
  object,
  threshold = 0,
  node_size = 16,
  node_color = "lightblue",
  label_edges = TRUE,
  label_pos = 0.35,
  ...
)

Arguments

Value

A ggplot object

Examples

if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("igraph", quietly = TRUE)) {
  library(ggplot2)
  dat <- generate_parce_sample(n = 500, seed = 42)
  model <- lingam_parce(dat$data)
  autoplot(model)
}

Plot the causal graph of an RCDResult with ggplot2

Description

Draws the estimated causal structure as a ggplot2-based directed graph, like autoplot.LingamResult(). Variable pairs suspected to share a latent confounder (NA entries in the adjacency matrix) are drawn as dashed, arrowless segments.

Usage

## S3 method for class 'RCDResult'
autoplot(
  object,
  threshold = 0,
  node_size = 16,
  node_color = "lightblue",
  label_edges = TRUE,
  label_pos = 0.35,
  ...
)

Arguments

Value

A ggplot object

Examples

if (requireNamespace("ggplot2", quietly = TRUE) &&
    requireNamespace("igraph", quietly = TRUE)) {
  library(ggplot2)
  confounded <- generate_rcd_sample(n = 300, seed = 1)
  model <- lingam_rcd(confounded$data)
  autoplot(model)
}

Bootstrap with Multiple Imputation for Direct LiNGAM

Description

Causal discovery on data containing missing values (NA). Each bootstrap resample (drawn with replacement, missing values retained) is multiply imputed into n_repeats complete datasets, and a common causal structure is jointly estimated across those datasets with lingam_multi_group() (Shimizu 2012), treating the imputed copies as "groups" sharing one causal order. R port of the Python lingam.tools.bootstrap_with_imputation().

Usage

bootstrap_with_imputation(
  X,
  n_sampling,
  n_repeats = 10L,
  imputer = NULL,
  cd_fit = NULL,
  prior_knowledge = NULL,
  apply_prior_knowledge_softly = FALSE,
  seed = NULL,
  verbose = TRUE
)

Arguments

Details

Procedure: for each of n_sampling iterations, (1) resample X with replacement (missing values are retained), (2) impute the resample into n_repeats complete datasets, (3) jointly estimate one causal structure shared by all n_repeats datasets with lingam_multi_group(). This assumes the same causal structure underlies every imputed copy.

Default imputer. mice::mice(method = "norm") (Bayesian linear regression, multiple imputation by chained equations) is the closest standard R analogue of the upstream Python default (IterativeImputer(sample_posterior = TRUE)). The two do not produce numerically identical imputations.

Custom imputer / cd_fit. Supply your own imputation or causal-discovery routine by passing a function with the signature described above; the return value is validated and a descriptive error is raised on violation. This replaces the abstract base classes (BaseMultipleImputation, BaseMultiGroupCDModel) of the Python original.

Downstream analysis. The result's shape (an extra n_repeats dimension for adjacency_matrices and imputation_results) differs from lingam_direct_bootstrap()'s BootstrapResult, so it cannot be passed directly to get_probabilities() etc. Use as_bootstrap_result() to collapse the n_repeats dimension (aggregating by median or mean) into a BootstrapResult.

On iteration failures: each iteration is wrapped in tryCatch(); a failing iteration (e.g. mice fails to converge on a particular resample) is skipped with a warning, and only if every iteration fails is an error raised, mirroring lingam_direct_bootstrap().

Sequential execution only. Unlike lingam_direct_bootstrap(), this function does not support parallel = TRUE; the upstream Python implementation is sequential as well. If needed in the future, it can be parallelized following the parallel::makePSOCKcluster() pattern used by lingam_direct_bootstrap().

Value

An ImputationBootstrapResult (list) containing:

• causal_orders: n_sampling x p integer matrix (1-based causal order per iteration).

• adjacency_matrices: array(n_sampling, n_repeats, p, p); ⁠[, , i, j]⁠ follows the lingamr convention (B[i, j] = coefficient of j -> i).

• resampled_indices: n_sampling x n integer matrix of resampled row indices.

• imputation_results: array(n_sampling, n_repeats, n, p); non-NA only at positions that were missing in that iteration's bootstrap resample.

Examples

set.seed(1)
sample6 <- generate_lingam_sample_6(n = 300, seed = 1)
X <- sample6$data
X$x5[sample.int(nrow(X), size = round(0.1 * nrow(X)))] <- NA # MCAR 10% on x5


if (requireNamespace("mice", quietly = TRUE)) {
  res <- bootstrap_with_imputation(X,
    n_sampling = 5L, n_repeats = 3L, seed = 42, verbose = FALSE
  )
  print(res)

  # Collapse the n_repeats dimension to reuse the existing bootstrap tooling
  bs <- as_bootstrap_result(res, aggregate = "median")
  get_probabilities(bs)
}

Build the RCD adjacency matrix from parents and confounder pairs

Description

Faithful port of rcd.py's adjacency-matrix construction (original 368-408 lines). B[i, j] is the coefficient of j -> i, matching the lingamr convention; no transpose is needed.

Usage

build_adjacency_matrix_rcd(X, P, C)

Arguments

Value

adjacency matrix B (n_features x n_features), with NA entries for confounder pairs

Build the lagged design matrix for VAR models

Description

Constructs Z such that column block k (columns (k-1)*p+1 to k*p) contains ⁠X_{t-k}⁠ for ⁠t = lags+1, ..., n⁠.

Usage

build_lag_matrix(X, lags)

Arguments

Value

matrix of shape ⁠(n - lags, lags * p)⁠

Build a lavaan model string from an adjacency matrix

Description

Converts a lingamr-convention adjacency matrix (B[i, j] = causal coefficient from j to i) into a lavaan model syntax string. Non-zero elements become regression paths (xi ~ xj); NA elements (used by e.g. lingam_parce() to mark a latent confounder between two variables) become a residual covariance (xi ~~ xj) between the two machine-named variables, which is the standard lavaan idiom equivalent to a two-indicator latent common cause with one loading fixed (as used by the Python semopy-based original).

Usage

build_lavaan_model(B, var_names)

Arguments

Value

character scalar, lavaan model syntax (possibly empty string)

Build the true adjacency matrix from edge specifications

Description

Build the true adjacency matrix from edge specifications

Usage

build_true_adjacency(var_names, from, to, coef)

Arguments

Value

adjacency matrix (p x p). m[to, from] = coef (row = to, col = from)

Unconditional non-Gaussianity statistic (tau)

Description

Unconditional non-Gaussianity statistic (tau)

Usage

calc_tau(k, pa, ch)

Arguments

Value

a single non-negative numeric statistic

Conditional non-Gaussianity statistics (tau), minimized over all conditioning subsets

Description

Per Wang & Drton (2020), the pruning statistic for a candidate variable is the minimum of the conditional tau statistic over every size-appropriate conditioning subset. Upstream HighDimDirectLiNGAM (cdt15/lingam) has a return mis-indented inside its loop over conditioning sets, so it only ever evaluates the first subset; this R port intentionally does NOT replicate that bug and instead evaluates every subset in cond_sets, so results differ numerically from the Python package (see dev/high-dim-direct-lingam-implementation.md).

Usage

calc_taus(Y, yty, pa, ch, k, cond_sets, an_sets)

Arguments

Value

numeric vector of length ncol(Y), minimum statistic per variable

Compute the total causal effect between two variables from an adjacency matrix

Description

Returns the sum of all path effects enumerated by find_all_paths(). Returns 0 if no path exists.

Usage

calculate_total_effect(adjacency_matrix, from_index, to_index)

Arguments

Value

Total causal effect (scalar)

Check whether glmnet is available

Description

If it is not available, raise an error indicating which regression method required it.

Usage

check_glmnet_available(method)

Arguments

Check whether lavaan is available

Description

Check whether lavaan is available

Usage

check_lavaan_available()

Check that the 'mice' package is available

Description

Check that the 'mice' package is available

Usage

check_mice_available(context)

Arguments

Check the stationarity of a fitted VAR-LiNGAM model

Description

Recovers the reduced-form VAR coefficients ⁠M_k = (I - B0)^{-1} B_k⁠ from the structural matrices and inspects the eigenvalues of the VAR companion matrix. The process is stationary when every eigenvalue lies strictly inside the unit circle (all moduli < 1); a modulus on or outside it signals a (near-)unit root, under which the VAR-LiNGAM estimates are unreliable.

Usage

check_var_stationarity(result, tol = 1)

Arguments

Value

a var_stationarity object (list) with moduli (sorted descending), max_modulus, is_stationary (logical), lags, and tol.

References

Stationarity diagnostics in the spirit of the VARLiNGAM R code of Moneta, A., Entner, D., Hoyer, P. O., & Coad, A. (2013), Oxford Bulletin of Economics and Statistics, 75(5), 705-730. https://sites.google.com/site/dorisentner/publications/VARLiNGAM

Examples

s <- generate_varlingam_sample(n = 1000, seed = 42)
m <- lingam_var(s$data, lags = 1, reg_method = "ols", prune = FALSE)
check_var_stationarity(m)

Residual matrix to diagnose for a VAR-LiNGAM model

Description

Returns the series targeted by the residual diagnostics: either the LiNGAM innovations ⁠e_t = (I - B0) n_t⁠ (the independent errors) or the reduced-form VAR residuals n_t. Shared by the normality tests and the QQ plot.

Usage

compute_varlingam_residuals(result, on = c("innovations", "var"))

Arguments

Value

residual matrix (n_obs x n_features), column names preserved

Signal a hook contract violation (imputer / cd_fit return value)

Description

A plain error, but always called from outside the per-iteration tryCatch() in bootstrap_with_imputation() (see validate_imputer_output() / validate_cd_fit_output() call sites there), so that a programming error in a user-supplied hook aborts the whole call immediately instead of being swallowed as a per-iteration stochastic estimation failure.

Usage

contract_violation(msg)

Arguments

Create a BootstrapResult

Description

Create a BootstrapResult

Usage

create_bootstrap_result(
  adjacency_matrices,
  total_effects,
  resampled_indices = NULL,
  causal_orders = NULL
)

Arguments

Value

BootstrapResult (list with class attribute)

Create a VARBootstrapResult

Description

Create a VARBootstrapResult

Usage

create_var_bootstrap_result(
  adjacency_matrices,
  total_effects,
  lags,
  resampled_indices = NULL,
  causal_orders = NULL
)

Arguments

Value

a VARBootstrapResult (list with class attribute)

Default causal-discovery fit: joint estimation via lingam_multi_group()

Description

Treats the imputed datasets as "groups" sharing a common causal order.

Usage

default_cd_fit(X_list, prior_knowledge, apply_prior_knowledge_softly)

Arguments

Value

⁠list(causal_order = <integer vector>, adjacency_matrices = <list>)⁠

Default imputer: multiple imputation via mice::mice(method = "norm")

Description

The closest standard R equivalent of the upstream Python default (sklearn.impute.IterativeImputer(sample_posterior = TRUE), Bayesian linear regression chained equations with posterior sampling) is mice::mice(method = "norm") (Bayesian linear regression). Numeric values will not match the Python implementation; the multiple-imputation design is equivalent.

Usage

default_imputer(X_boot, n_repeats)

Arguments

Value

A list of n_repeats complete numeric matrices

Difference of mutual information

Description

Difference of mutual information

Usage

diff_mutual_info(xi_std, xj_std, ri_j, rj_i)

Arguments

Value

Difference of mutual information

Maximum-entropy approximation of entropy

Description

Maximum-entropy approximation of entropy

Usage

entropy_approx(u)

Arguments

Value

Approximate entropy value

Estimate the adjacency matrix from a causal order

Description

Estimate the adjacency matrix from a causal order

Usage

estimate_adjacency_matrix(
  X,
  causal_order,
  prior_knowledge = NULL,
  method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Value

adjacency matrix B (n_features x n_features)

Estimate the adjacency matrix by causal order using the n <= p route

Description

Estimate the adjacency matrix by causal order using the n <= p route

Usage

estimate_adjacency_matrix_high_dim_np(X, causal_order)

Arguments

Value

adjacency matrix B (n_features x n_features)

Estimate the adjacency matrix from a ParceLiNGAM causal order

Description

causal_order is a list whose first element may be an unresolved block (length > 1); all remaining elements are length-1. Block members are never regression targets (their parents cannot be identified), but they are valid predictors for downstream variables. Pairs within the block are set to NA.

Usage

estimate_adjacency_matrix_parce(
  X,
  causal_order,
  prior_knowledge,
  method,
  lambda,
  init_method
)

Arguments

Value

adjacency matrix B (n_features x n_features)

Estimate the total causal effects between all variables at once

Description

Estimate the total causal effects between all variables at once

Usage

estimate_all_total_effects(
  X,
  lingam_result,
  method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Value

Matrix of total causal effects (n_features x n_features). Convention: TE[i, j] is the total causal effect from variable j to variable i (j -> i). Same index convention as the adjacency matrix adjacency_matrix. The sum of direct and indirect effects.

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

model <- LiNGAM_sample_1000$data |>
  lingam_direct(reg_method = "ols")

LiNGAM_sample_1000$data |>
  estimate_all_total_effects(model)

Estimate the total causal effect between two specified variables

Description

Estimate the total causal effect between two specified variables

Usage

estimate_total_effect(
  X,
  lingam_result,
  from_index,
  to_index,
  method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Value

Estimated total causal effect

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

model <- LiNGAM_sample_1000$data |>
  lingam_direct(reg_method = "ols")

LiNGAM_sample_1000$data |>
  estimate_total_effect(model, 4, 1)

Estimate the total causal effect between two variables (ParceLiNGAM)

Description

Analogous to estimate_total_effect(), but for lingam_parce() results, which may contain NA entries in the adjacency matrix.

Usage

estimate_total_effect_parce(
  X,
  parce_result,
  from_index,
  to_index,
  method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Value

Estimated total causal effect, or NA (with a warning) if from_index is part of an unresolved block (its parents cannot be identified).

Examples

confounded <- generate_parce_sample(n = 500, seed = 1)
result <- lingam_parce(confounded$data, reg_method = "ols")

# A well-identified pair returns a numeric estimate
estimate_total_effect_parce(confounded$data, result, from_index = 1, to_index = 5)

Estimate the total causal effect between two variables (RCD)

Description

Analogous to estimate_total_effect(), but for lingam_rcd() results, which may contain NA entries in the adjacency matrix.

Usage

estimate_total_effect_rcd(
  X,
  rcd_result,
  from_index,
  to_index,
  method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Value

Estimated total causal effect, or NA (with a warning) if from_index is part of a suspected latent confounder pair (its parents cannot be identified). Also warns (without altering the estimate) if to_index is an ancestor of from_index according to ancestors_list, since that is inconsistent with a from -> to effect.

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)
result <- lingam_rcd(confounded$data)

# A well-identified pair returns a numeric estimate
estimate_total_effect_rcd(confounded$data, result, from_index = 6, to_index = 1)

Estimate a total causal effect in a VAR-LiNGAM model

Description

Estimates the total causal effect from from_index (optionally at lag from_lag) to to_index (at the current time) using the fitted VAR-LiNGAM model. Port of the Python reference estimate_total_effect: the destination variable is regressed on the source variable together with the source's parents (a back-door adjustment), and the source's coefficient is returned.

Usage

estimate_var_total_effect(X, result, from_index, to_index, from_lag = 0)

Arguments

Value

the estimated total effect (scalar)

Examples

sample <- generate_varlingam_sample(n = 1000, seed = 42)
model <- lingam_var(sample$data, lags = 1, reg_method = "ols", prune = FALSE)

# total effect of x0 (current) on x2 (current)
estimate_var_total_effect(sample$data, model, from_index = 1, to_index = 3)

Evaluate model fit of an estimated causal graph

Description

Fits the causal graph implied by adjacency_matrix as a structural equation model (SEM) via lavaan::sem() and returns standard SEM fit measures (CFI, RMSEA, AIC/BIC, etc.). This is an R port of the Python lingam.utils.evaluate_model_fit(), which delegates to the Python package semopy; this R version delegates to lavaan instead.

Usage

evaluate_model_fit(adjacency_matrix, X, is_ordinal = NULL)

Arguments

Details

• Optional dependency: this function requires the lavaan package (listed in Suggests, not Imports). Install it with install.packages("lavaan").

• Latent confounders: an NA element B[i, j] (as produced by e.g. lingam_parce() for a suspected latent confounder between variables i and j) is represented as a residual covariance xi ~~ xj in the lavaan model. This is algebraically equivalent to the two-indicator latent common cause (one loading fixed to 1) used by the Python semopy implementation, but is expressed with lavaan's standard residual-covariance idiom rather than an explicit latent variable.

• Numerical values will not match semopy exactly: lavaan and semopy use different default estimators/options. The fit measures returned are the same statistics, but exact numbers can differ.

• Convention: adjacency_matrix follows the lingamr convention B[i, j] = causal coefficient from variable j to variable i (j -> i).

Value

A one-row data.frame of fit measures: DoF, DoF Baseline, chi2, chi2 p-value, chi2 Baseline, CFI, GFI, AGFI, NFI, TLI, RMSEA, AIC, BIC, LogLik. When is_ordinal is used, AIC/BIC/LogLik and some other measures are not defined by the WLSMV estimator and are returned as NA.

References

Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1-36. doi:10.18637/jss.v048.i02

Examples

if (requireNamespace("lavaan", quietly = TRUE)) {
  dat <- generate_lingam_sample_6()
  result <- lingam_direct(dat$data, reg_method = "ols")
  evaluate_model_fit(result, dat$data)
}

Whether xi can be excluded from sink candidacy given the ancestor sets known so far

Description

Faithful port of the exclusion rule in rcd.py (original 165-174 lines): xi cannot be the sink of U if (a) xi is already known to be an ancestor of some other member of U, or (b) every other member of U is already known to be an ancestor of xi.

Usage

exists_ancestor_in_U(M, U, xi, xj_list)

Arguments

Value

TRUE if xi should be excluded from sink candidacy

Extract an adjacency matrix from a result object or matrix

Description

Extract an adjacency matrix from a result object or matrix

Usage

extract_adjacency_matrix(adjacency_matrix)

Extract the ancestor sets M(x_i) for every variable

Description

Faithful port of rcd.py's extract_ancestors() (original 254-316 lines). Repeatedly scans variable subsets of increasing size, growing each variable's ancestor set whenever a unique "sink" is found within a subset. See dev/rcd-implementation.md section 3.2 for the full algorithm description and the reasoning behind the cache-update position.

Usage

extract_ancestors(
  X,
  max_explanatory_num,
  cor_alpha,
  ind_alpha,
  shapiro_alpha,
  MLHSICR,
  independence,
  ind_corr
)

Arguments

Value

list of length ncol(X); each element is a sorted integer vector of ancestor indices (possibly empty)

Extract parents from the ancestor sets

Description

Faithful port of rcd.py's extract_parents() (original 318-343 lines).

Usage

extract_parents(X, M, cor_alpha)

Arguments

Value

list of length ncol(X); each element is a sorted integer vector of parent indices (subset of the corresponding ancestor set)

Extract partial orders from prior knowledge

Description

Extract partial orders from prior knowledge

Usage

extract_partial_orders(pk)

Arguments

Value

matrix (n x 2), each row is a (from, to) partial order

Detect variable pairs sharing an unobserved latent confounder

Description

Faithful port of rcd.py's extract_vars_sharing_confounders() (original 352-366 lines). Only pairs with no parent-child relationship in either direction are considered.

Usage

extract_vars_sharing_confounders(X, P, cor_alpha)

Arguments

Value

list of length ncol(X); each element is a sorted integer vector of indices sharing a latent confounder with that variable (symmetric)

F-correlation (kernel canonical correlation) between two variables

Description

Bach & Jordan (2002) kernel canonical correlation, as used by BottomUpParceLiNGAM's independence = "fcorr" option. Returns a value in (roughly) ⁠[0, 1]⁠; larger means more dependent.

Usage

f_correlation(x, y)

Arguments

Details

A constant x or y (zero variance) carries no dependence information and is treated as trivially independent (returns 0) rather than propagating a division-by-zero standardization into NaN/Inf.

Value

F-correlation value (scalar)

Low-rank SVD transform used inside f_correlation()

Description

Given the (already column-centered) incomplete-Cholesky factor G, returns the orthonormalized basis U and the shrinkage vector R used to assemble the block canonical-correlation matrix R_kappa.

Usage

fcorr_svd_transform(G, kappa, n)

Arguments

Value

list(U = n x d' matrix, R = length d' vector); R has length 0 if no eigenvalue of crossprod(G) clears the kappa threshold (i.e. G carries no retainable rank), which callers must handle explicitly.

Enumerate all paths in a DAG via depth-first search

Description

Takes an adjacency matrix where B[i, j] represents j -> i, and returns all paths from from_index to to_index together with each path effect (the product of the coefficients).

Usage

find_all_paths(adjacency_matrix, from_index, to_index, min_causal_effect = 0)

Arguments

Value

list(paths, effects)

Find the most sink-like candidate variable

Description

Faithful port of ⁠_search_causal_order.find_exo_vec()⁠ (original 227-276 lines). For each candidate j in Uc, regresses j on the other candidates (setdiff(Uc, j), not U) and evaluates how independent the residual is from those explanatory variables.

Usage

find_exo_vec(X, Uc, U, independence, Cov)

Arguments

Value

list(m = selected variable index, eval = its evaluation value)

Adaptive LASSO

Description

Adaptive LASSO

Usage

fit_adaptive_lasso(
  y,
  Xp,
  lambda = "BIC",
  gamma_weight = 1,
  init_method = "ols"
)

Arguments

Value

coefficient vector

LASSO regression (lambda selection by information criterion or CV)

Description

LASSO regression (lambda selection by information criterion or CV)

Usage

fit_lasso(y, Xp, lambda = "BIC")

Arguments

Value

coefficient vector

OLS regression

Description

OLS regression

Usage

fit_ols(y, Xp)

Arguments

Value

coefficient vector (excluding intercept)

OLS fit for a single predictor, pruned by information criterion

Description

glmnet requires at least two predictor columns, so the penalized methods fall back to OLS when only one predictor remains. Plain OLS never yields an exact zero, which would make single-predictor edges unprunable: the second variable in the causal order always has exactly one predictor, so a spurious edge would survive even for fully independent data. To preserve the sparse behavior of the penalized methods, the OLS coefficient is kept only when adding the predictor improves the information criterion over the intercept-only model; otherwise it is set to exactly zero.

Usage

fit_ols_ic_pruned(y, Xp, lambda)

Arguments

Details

The criterion is AIC for lambda = "AIC" and BIC otherwise (the CV / oracle lambdas have no single-predictor counterpart, so the sparsest criterion, BIC, is used for them as well).

Value

length-1 coefficient vector (0 when the predictor is pruned)

Penalized regression via glmnet (IC or CV lambda selection)

Description

Internal helper shared by fit_lasso() and fit_ridge_reg(). Both functions differ only in alpha and lambda_seq; this function encapsulates the duplicated IC / CV branches.

Usage

fit_penalized_regression(y, Xp_mat, alpha, lambda, lambda_seq)

Arguments

Value

coefficient vector (excluding intercept)

Ridge regression (lambda selection by information criterion or CV)

Description

Ridge regression (lambda selection by information criterion or CV)

Usage

fit_ridge_reg(y, Xp, lambda = "BIC")

Arguments

Value

coefficient vector

Fit a VAR(p) model by OLS (no intercept)

Description

Fit a VAR(p) model by OLS (no intercept)

Usage

fit_var_ols(X, lags)

Arguments

Value

list with coefs (array (lags, n_features, n_features); coefs[k, , ] is M_k such that ⁠X_t = sum_k M_k X_{t-k} + e_t⁠) and residuals (n_samples - lags, n_features)

Generate sample data for LiM (3 mixed variables)

Description

Generates a small dataset with a known causal chain of continuous and binary (0/1) discrete variables: x1 (continuous) -> x2 (discrete) -> x3 (continuous). Continuous variables use Laplace-distributed noise (non-Gaussian, as required by LiNGAM-family methods); the discrete variable is drawn from a Bernoulli distribution whose logit is a linear function of its parent.

Usage

generate_lim_sample(n = 1000L, seed = NULL)

Arguments

Value

A list with three elements:

• data: data frame with columns x1, x2, x3.

• adjacency_matrix: 3x3 true adjacency matrix, following the lingamr convention (m[to, from], i.e. adjacency_matrix["x2", "x1"] is the coefficient of the x1 -> x2 edge).

• is_continuous: logical vector c(TRUE, FALSE, TRUE) marking which columns of data are continuous.

Examples

dat <- generate_lim_sample(n = 500, seed = 1)
result <- lingam_lim(dat$data, is_continuous = dat$is_continuous)
result$adjacency_matrix

Generate a challenging sample data for Direct LiNGAM

Description

Generates a dataset with conditions that make causal estimation difficult:

1. High multicollinearity among predictors

2. Moderate sample size relative to variables

3. True coefficients of similar magnitude

Usage

generate_lingam_hard_sample(n = 200L, seed = 42L, collinearity = 0.95)

Arguments

Details

These conditions destabilize OLS initial estimates in Adaptive LASSO, making Ridge-initialized Adaptive LASSO preferable.

Value

list(data, true_adjacency)

Examples

result <- generate_lingam_hard_sample()
result$true_adjacency
head(result$data)

Generate large-scale sample data to benchmark Direct LiNGAM scalability

Description

Generates a dataset with many variables to demonstrate the computational scalability difference between Direct LiNGAM and ICA-LiNGAM.

Usage

generate_lingam_large_sample(
  p = 20L,
  n = 1000L,
  max_parents = 3L,
  coef_min = 0.5,
  coef_max = 1.5,
  seed = 42L,
  noise_dist = "uniform"
)

Arguments

Details

Why Direct LiNGAM slows down with large p

At each of its p steps, Direct LiNGAM evaluates an independence measure between every remaining candidate root and every other residual. The total number of evaluations is:

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

i.e., O(p^3). Each evaluation is itself O(n), giving O(p^3 n) overall. For p = 10 this is about 330 evaluations; for p = 20 about 2,660; for p = 40 about 21,320 — an 8x increase every time p doubles.

Why ICA-LiNGAM scales better

ICA-LiNGAM applies FastICA once to the whole p x n data matrix. Each FastICA iteration costs O(p^2 n), and the algorithm typically converges in far fewer than p iterations. Additionally, these matrix operations are fully vectorised (BLAS/LAPACK), whereas Direct LiNGAM iterates over pairs in an R loop.

Data-generating process

Variables are topologically ordered as x0, x1, ..., x(p-1). For each i >= 1, the number of parents is sampled uniformly from 1 to min(max_parents, i), and the parents are drawn without replacement from x0, ..., x(i-1). Edge coefficients are drawn uniformly from [-coef_max, -coef_min] U [coef_min, coef_max]. The resulting adjacency matrix is strictly lower-triangular.

Value

A list with three elements:

• data: data.frame with p columns (x0, x1, ..., x(p-1)).

• true_adjacency: p x p matrix. true_adjacency[i, j] is the structural coefficient of the edge xj -> xi (row = to, col = from). The matrix is strictly lower-triangular because variables are stored in causal order.

• true_causal_order: integer vector 0:(p-1). Variables are already in topological order, so the true causal order is always 0, 1, ..., p-1.

Examples

# Generate 20-variable data and check its sparsity
dataset <- generate_lingam_large_sample(p = 20, n = 500)
dim(dataset$data)                    # 500 x 20
sum(dataset$true_adjacency != 0)     # number of edges
dataset$true_causal_order            # 0, 1, ..., 19


# As the number of variables grows, Direct LiNGAM's run time increases sharply
t10 <- system.time(lingam_direct(generate_lingam_large_sample(p = 10)$data))
t20 <- system.time(lingam_direct(generate_lingam_large_sample(p = 20)$data))
cat(sprintf("p=10: %.1f sec,  p=20: %.1f sec\n", t10["elapsed"], t20["elapsed"]))

Generate Paradoxical Data Where DirectLiNGAM Struggles

Description

Generates a synthetic dataset designed to favor ICA-LiNGAM (due to standardized scales) while challenging DirectLiNGAM (due to heavy measurement noise on the root variable, which triggers error propagation). The true causal structure is a serial chain: x0 -> x1 -> x2 -> x3 (each coefficient 0.8).

Usage

generate_lingam_paradox_data(n = 2000L, seed = 42L)

Arguments

Details

This function intentionally injects strong measurement error into the root (causal upstream) variable x0. This noise corrupts the independence tests performed at the initial step of DirectLiNGAM, frequently causing it to misidentify the root variable and leading to a cascading failure (error propagation) throughout the causal ordering.

On the other hand, the output data is completely standardized using the scale() function. This eliminates any differences in scale among the variables, thereby neutralizing the major weakness of ICA-LiNGAM (scale-dependence) and allowing it to perform relatively better.

Because the data are standardized and the root carries measurement error, the coefficients estimated by lingam_direct() will not exactly match the 0.8 values stored in true_adjacency.

Value

list(data, true_adjacency)

• data: a data frame with 4 standardized variables (x0, x1, x2, x3); each column has a mean of 0 and a standard deviation of 1.

• true_adjacency: the 4x4 true adjacency matrix of the data-generating chain, following the m[row = to, col = from] convention and holding the structural coefficients (0.8) on the latent, pre-standardization scale.

Examples

# Generate the dataset
paradox <- generate_lingam_paradox_data(n = 1000, seed = 123)

# Verify the dataset
head(paradox$data)
sapply(paradox$data, sd)

# True data-generating structure
paradox$true_adjacency

Generate 10-variable sample data for Direct LiNGAM

Description

Generates a sample dataset with a known causal structure. The true causal structure is: x3 -> x0 (coef = 3.0) x3 -> x2 (coef = 6.0) x3 -> x9 (coef = 7.0) x0 -> x1 (coef = 3.0) x0 -> x5 (coef = 4.0) x0 -> x4 (coef = 8.0) x0 -> x7 (coef = 3.0) x2 -> x1 (coef = 2.0) x2 -> x4 (coef = -1.0) x2 -> x8 (coef = 0.5) x1 -> x6 (coef = 2.0) x5 -> x8 (coef = 2.0) x4 -> x7 (coef = 1.5) x6 -> x9 (coef = 1.0)

Usage

generate_lingam_sample_10(n = 1000L, seed = 42L, noise_dist = "uniform")

Arguments

Value

list(data, true_adjacency)

Examples

# Non-Gaussian (LiNGAM works well)
X_nongauss <- generate_lingam_sample_10(noise_dist = "uniform")
result <- lingam_direct(X_nongauss$data, reg_method = "ols")
result$causal_order

# Gaussian (LiNGAM may fail)
X_gauss <- generate_lingam_sample_10(noise_dist = "gaussian")
result <- lingam_direct(X_gauss$data, reg_method = "ols")
result$causal_order

Generate sample data for Direct LiNGAM (6 variables)

Description

Generates a sample dataset with a known causal structure. The true causal structure is: x3 -> x0 (coef = 3.0) x3 -> x2 (coef = 6.0) x0 -> x1 (coef = 3.0) x2 -> x1 (coef = 2.0) x0 -> x5 (coef = 4.0) x0 -> x4 (coef = 8.0) x2 -> x4 (coef = -1.0)

Usage

generate_lingam_sample_6(n = 1000L, seed = 42L, noise_dist = "uniform")

Arguments

Value

list(data, true_adjacency)

Examples

# Non-Gaussian (LiNGAM works well)
X_nongauss <- generate_lingam_sample_6(noise_dist = "uniform")
result <- lingam_direct(X_nongauss$data, reg_method = "ols")
result$causal_order

# Gaussian (LiNGAM may fail)
X_gauss <- generate_lingam_sample_6(noise_dist = "gaussian")
result <- lingam_direct(X_gauss$data, reg_method = "ols")
result$causal_order

Generate sample data for Multi-Group Direct LiNGAM (2 groups, 6 variables)

Description

Generates two datasets that share the same causal structure as generate_lingam_sample_6() (⁠x3 -> x0, x3 -> x2, x0 -> x1, x2 -> x1, x0 -> x4, x2 -> x4, x0 -> x5⁠) but with different structural coefficients per group, following the multi-dataset tutorial's setup.

Usage

generate_multi_group_sample(n = c(1000, 1000), seed = 42L)

Arguments

Value

A list with:

• data_list: named list (group1, group2) of data frames, each with columns x0..x5.

• adjacency_matrices: named list of the true adjacency matrix per group (m[to, from] = coef, same convention as lingam_direct()).

• causal_order: the true causal order shared by both groups (1-based indices into x0..x5).

Examples

mg <- generate_multi_group_sample()
lapply(mg$data_list, head)
mg$adjacency_matrices$group1

Generate a noise matrix with an independent seed per variable

Description

Column k is generated by noise_fn(n) immediately after set.seed(seed + k - 1). Fixing the seed per variable ensures that the same seed always yields the same noise columns.

Usage

generate_noise_matrix(n, n_vars, seed, noise_fn)

Arguments

Value

noise matrix (n x n_vars)

Generate sample data with a latent confounder (for BottomUpParceLiNGAM)

Description

Generates the 7-variable model used in the ParceLiNGAM tutorial, where x6 is an unobserved (latent) common cause of x2 and x3. Only x0-x5 are returned as observed data; x6 is not included.

Usage

generate_parce_sample(n = 1000L, seed = NULL)

Arguments

Details

The data-generating process (all error terms are Uniform(0, 1)):

x6 (latent) ~ Uniform(0, 1)
x3 = 2.0 * x6 + e
x2 = 2.0 * x6 + e
x0 = 0.5 * x3 + e
x1 = 0.5 * x0 + 0.5 * x2 + e
x5 = 0.5 * x0 + e
x4 = 0.5 * x0 - 0.5 * x2 + e

Value

list with three elements:

• data: data.frame of the 6 observed variables (x0-x5).

• adjacency_matrix: the true 6x6 adjacency matrix among the observed variables, following the m[to, from] convention. The x2-x3 entries (which share the latent confounder x6 and have no direct edge between them) are NA, matching the convention used by lingam_parce().

• confounded_pair: 1-based column positions of x2 and x3 (the pair sharing the latent confounder).

Examples

confounded <- generate_parce_sample(n = 500, seed = 42)
head(confounded$data)
confounded$adjacency_matrix
confounded$confounded_pair

Generate sample data with a latent confounder (for RCD)

Description

Generates the 7-variable model used in the RCD tutorial, where x6 is an unobserved (latent) common cause of x2 and x4. Only x0-x5 are returned as observed data; x6 is not included.

Usage

generate_rcd_sample(n = 300L, seed = NULL)

Arguments

Details

The data-generating process (all error terms e() are super-Gaussian, rnorm(n, 0, 0.5)^3):

x5 ~ e();  x6 (latent) ~ e()
x1 = 0.6 * x5 + e()
x3 = 0.5 * x5 + e()
x0 = 1.0 * x1 + 1.0 * x3 + e()
x2 = 0.8 * x0 - 0.6 * x6 + e()
x4 = 1.0 * x0 - 0.5 * x6 + e()

Value

list with four elements:

• data: data.frame of the 6 observed variables (x0-x5).

• adjacency_matrix: the true 6x6 adjacency matrix among the observed variables, following the m[to, from] convention. The x2-x4 entries (which share the latent confounder x6 and have no direct edge between them) are NA, matching the convention used by lingam_rcd().

• ancestors_list: the true ancestor sets (1-based column positions), usable as a test oracle for lingam_rcd()'s ancestors_list.

• confounded_pair: 1-based column positions of x2 and x4 (the pair sharing the latent confounder).

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)
head(confounded$data)
confounded$adjacency_matrix
confounded$ancestors_list
confounded$confounded_pair

Generate sample data from a VAR-LiNGAM model

Description

Generates a 3-variable time series following a VAR-LiNGAM model with a strictly acyclic instantaneous structure B0, a lag-1 coefficient matrix M1, and non-Gaussian (uniform) errors.

Usage

generate_varlingam_sample(n = 1000, seed = NULL)

Arguments

Value

list with data (data frame, n x 3), true_B0 (instantaneous matrix), and true_M1 (lag-1 coefficient matrix)

Examples

sample <- generate_varlingam_sample(n = 500, seed = 1)
head(sample$data)

Create an adjacency matrix of representative causal-effect values from bootstrap results

Description

Create an adjacency matrix of representative causal-effect values from bootstrap results

Usage

get_adjacency_matrix_summary(
  result,
  stat = "median",
  min_causal_effect = NULL,
  min_probability = NULL,
  labels = NULL
)

Arguments

Value

Adjacency matrix (n_features x n_features). Rule: B[i, j] is the causal coefficient from variable j to variable i (j -> i). Same rule as the adjacency_matrix of lingam_direct().

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L, reg_method = "ols", seed = 42
)
get_adjacency_matrix_summary(bs_model)

Get counts, proportions, and causal effects of causal directions

Description

Get counts, proportions, and causal effects of causal directions

Usage

get_causal_direction_counts(
  result,
  n_directions = NULL,
  min_causal_effect = NULL,
  split_by_causal_effect_sign = FALSE,
  labels = NULL
)

Arguments

Value

A data frame containing the following columns:

• from, to: 1-based indices of the causal (from) and effect (to) variables.

• count: Number of bootstrap samples in which this specific causal direction was identified.

• proportion: The frequency of the direction's occurrence (count / n_sampling), representing its bootstrap probability.

• mean_effect: The average value of the estimated causal effects across samples where this direction was identified.

• median_effect: The median value of the estimated causal effects, providing a robust estimate of the effect size.

• sd_effect: The standard deviation of the causal effect estimates, indicating the stability of the effect size.

• ci_lower, ci_upper: The lower (2.5%) and upper (97.5%) bounds of the bootstrap confidence interval for the causal effect.

• sign (optional): The sign of the causal effect (1 for positive, -1 for negative), included if split_by_causal_effect_sign = TRUE.

• from_name, to_name (optional): Character labels for the variables, included if labels were provided.

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L, reg_method = "ols", seed = 42
)

get_causal_direction_counts(bs_model, labels = names(LiNGAM_sample_1000$data))

Evaluate the stability of the causal order from bootstrap

Description

Aggregates the causal order (causal_order) estimated in each bootstrap sample and quantifies how stable the order is. Returns the rank distribution of each variable, the precedence probabilities for variable pairs, and an overall stability score.

Usage

get_causal_order_stability(result, labels = NULL)

Arguments

Value

A list of class causal_order_stability, containing:

• rank_summary: A summary of the rank of each variable (variable, mean_rank, sd_rank, median_rank, mode_rank). Sorted in ascending order of mean_rank (from upstream). A rank of 1 is the most upstream.

• precedence_matrix: A precedence probability matrix. P[i, j] is the proportion of bootstrap samples in which variable i was located upstream of (before) variable j.

• stability_score: An overall stability score, from 0 (random order) to 1 (order agrees across all samples). The closer the precedence probability of each variable pair is to 0/1, the higher the score.

• n_sampling: The number of bootstrap samples.

Examples

dat <- generate_lingam_sample_6()
bs <- lingam_direct_bootstrap(dat$data, n_sampling = 30L, reg_method = "ols", seed = 42)
get_causal_order_stability(bs, labels = names(dat$data))

Get DAG counts

Description

Get DAG counts

Usage

get_directed_acyclic_graph_counts(
  result,
  n_dags = NULL,
  min_causal_effect = NULL,
  split_by_causal_effect_sign = FALSE
)

Arguments

Value

list(dag = list of data.frames, count = integer vector)

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L, reg_method = "ols", seed = 42
)

get_directed_acyclic_graph_counts(bs_model)

Compute p-values for the independence test of the errors

Description

Compute p-values for the independence test of the errors

Usage

get_error_independence_p_values(X, lingam_result, method = "spearman")

Arguments

Value

matrix of p-values (n_features x n_features)

Examples

# Load the sample data
LiNGAM_sample_1000 <- generate_lingam_sample_6()

# Run Direct LiNGAM
result <- LiNGAM_sample_1000$data |>
  lingam_direct(reg_method = "ols")

# Compute p-values (default: Spearman)
p_vals <- get_error_independence_p_values(LiNGAM_sample_1000$data, result)
round(p_vals, 3)

# Compute with Kendall
p_vals_k <- get_error_independence_p_values(LiNGAM_sample_1000$data, result, method = "kendall")
round(p_vals_k, 3)

Compute p-values for the independence of ParceLiNGAM residuals (HSIC-based)

Description

Analogous to get_error_independence_p_values(), but for lingam_parce() results. Uses the HSIC gamma-approximation test (hsic_test_gamma()) rather than a correlation test, and returns NA for any pair involving a variable whose row or column in the adjacency matrix contains NA (residuals cannot be computed for those variables).

Usage

get_error_independence_p_values_parce(X, parce_result)

Arguments

Value

matrix of p-values (n_features x n_features)

Examples

confounded <- generate_parce_sample(n = 500, seed = 1)
result <- lingam_parce(confounded$data, reg_method = "ols")
round(get_error_independence_p_values_parce(confounded$data, result), 3)

Compute p-values for the independence of RCD residuals (HSIC-based)

Description

Analogous to get_error_independence_p_values_parce(), but for lingam_rcd() results. Returns NA for any pair involving a variable whose row or column in the adjacency matrix contains NA (residuals cannot be computed for those variables).

Usage

get_error_independence_p_values_rcd(X, rcd_result)

Arguments

Value

matrix of p-values (n_features x n_features)

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)
result <- lingam_rcd(confounded$data)
round(get_error_independence_p_values_rcd(confounded$data, result), 3)

Extract a single group's result from a MultiGroupLingamResult

Description

Returns the adjacency matrix and (shared) causal order of one group as a plain LingamResult, so that the existing single-group functions (estimate_total_effect(), estimate_all_total_effects(), get_error_independence_p_values(), plot_adjacency(), autoplot(), tidy()) can be applied to it directly.

Usage

get_group_result(x, group)

Arguments

Value

A LingamResult object (list) with adjacency_matrix and causal_order, identical in shape to the return value of lingam_direct().

Examples

mg <- generate_multi_group_sample()
res <- lingam_multi_group(mg$data_list, reg_method = "ols")
g1 <- get_group_result(res, 1)
class(g1)

Get all paths between two specified variables and their bootstrap probabilities

Description

Get all paths between two specified variables and their bootstrap probabilities

Usage

get_paths(result, from_index, to_index, min_causal_effect = NULL)

Arguments

Value

data.frame (path, effect, probability)

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L, reg_method = "ols", seed = 42
)
get_paths(bs_model, 1, 6)

Get bootstrap probabilities

Description

Get bootstrap probabilities

Usage

get_probabilities(result, min_causal_effect = NULL)

Arguments

Value

Probability matrix (n_features x n_features)

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L, reg_method = "ols", seed = 42
)

get_probabilities(bs_model)

Enumerate conditioning subsets and compute pruning statistics for one candidate variable

Description

Enumerate conditioning subsets and compute pruning statistics for one candidate variable

Usage

get_prune_stats(Y, yty, i, j, K, last_root, condition_set, J)

Arguments

Value

numeric vector of length ncol(Y)

Get a list of total causal effects

Description

Get a list of total causal effects

Usage

get_total_causal_effects(result, min_causal_effect = NULL)

Arguments

Value

data.frame (from, to, effect, probability)

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L, reg_method = "ols", seed = 42
)

get_total_causal_effects(bs_model)

Get variable names, falling back to x0, x1, ... when colnames is NULL

Description

Get variable names, falling back to x0, x1, ... when colnames is NULL

Usage

get_var_names(x)

Enumerate bootstrap paths between two variables in a VAR-LiNGAM model

Description

Builds the time-expanded graph for every bootstrap sample and enumerates all directed paths from the source (at from_lag) to the destination (at to_lag), reporting each path's bootstrap probability and median effect. Port of the Python reference VARBootstrapResult.get_paths.

Usage

get_var_paths(
  result,
  from_index,
  to_index,
  from_lag = 0,
  to_lag = 0,
  min_causal_effect = NULL
)

Arguments

Details

Node indices in the returned path are 1-based positions in the time-expanded graph: column j of block L (lag L) corresponds to index n_features * L + j.

Value

a data frame (path, effect, probability), one row per distinct path

Examples

s <- generate_varlingam_sample(n = 500, seed = 42)
bs <- lingam_var_bootstrap(s$data,
  n_sampling = 10L, criterion = NULL,
  reg_method = "ols", prune = FALSE, seed = 1, verbose = FALSE
)
get_var_paths(bs, from_index = 1, to_index = 3)

Bootstrap probabilities for a VAR-LiNGAM model

Description

Returns, for each entry of the joined adjacency matrix, the fraction of bootstrap samples in which that edge exceeded min_causal_effect.

Usage

get_var_probabilities(result, min_causal_effect = NULL)

Arguments

Value

probability matrix (n_features x n_features*(1 + lags)). Columns 1..n_features are the instantaneous block; the next n_features are lag 1; etc. P[i, j] is the probability of the edge j -> i.

Examples

s <- generate_varlingam_sample(n = 500, seed = 42)
bs <- lingam_var_bootstrap(s$data,
  n_sampling = 10L, criterion = NULL,
  reg_method = "ols", prune = FALSE, seed = 1, verbose = FALSE
)
get_var_probabilities(bs)

Get a one-row summary of a LiMResult

Description

Like glance.LingamResult(), with an additional n_discrete column giving the number of discrete variables in the model.

Usage

## S3 method for class 'LiMResult'
glance(x, ...)

Arguments

Value

A one-row data.frame(n_variables, n_edges, n_discrete, causal_order)

Examples

set.seed(1)
dat <- generate_lim_sample(n = 300)
model <- lingam_lim(dat$data, is_continuous = dat$is_continuous)
glance(model)

Get a one-row summary of a LingamResult

Description

Summarizes the entire model in a single row. The data X is not required because no residuals are computed. If residual-based diagnostics are needed, use summary_lingam() instead.

Usage

## S3 method for class 'LingamResult'
glance(x, ...)

Arguments

Value

A one-row data.frame(n_variables, n_edges, causal_order)

Examples

dat <- generate_lingam_sample_6()
model <- lingam_direct(dat$data, reg_method = "ols")
glance(model)

Get a one-row summary of a MultiGroupLingamResult

Description

Summarizes the joint model in a single row. The causal order is shared across groups; per-group edge counts are available via glance(get_group_result(x, i)).

Usage

## S3 method for class 'MultiGroupLingamResult'
glance(x, ...)

Arguments

Value

A one-row data.frame(n_groups, n_variables, causal_order)

Examples

mg <- generate_multi_group_sample()
model <- lingam_multi_group(mg$data_list, reg_method = "ols")
glance(model)

Get a one-row summary of a ParceLingamResult

Description

Like glance.LingamResult(). n_edges counts non-NA edges only, and n_na_entries counts the adjacency-matrix entries left NA (unresolved order / suspected latent confounding). Unresolved blocks in the causal order are shown in parentheses, as in the print method.

Usage

## S3 method for class 'ParceLingamResult'
glance(x, ...)

Arguments

Value

A one-row data.frame(n_variables, n_edges, n_na_entries, causal_order)

Examples

dat <- generate_parce_sample(n = 500, seed = 42)
model <- lingam_parce(dat$data)
glance(model)

Get a one-row summary of an RCDResult

Description

Like glance.LingamResult(), but without a causal order (RCD does not estimate one). n_edges counts non-NA edges only, and n_confounded_pairs counts the variable pairs whose adjacency-matrix entries are NA (suspected shared latent confounder).

Usage

## S3 method for class 'RCDResult'
glance(x, ...)

Arguments

Value

A one-row data.frame(n_variables, n_edges, n_confounded_pairs)

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)
model <- lingam_rcd(confounded$data)
glance(model)

Causal-order search for HighDimDirectLiNGAM

Description

Causal-order search for HighDimDirectLiNGAM

Usage

high_dim_causal_order(X, J, K, alpha)

Arguments

Value

integer vector, 1-based causal order (upstream-most first)

Gaussian Gram matrix and its double-centered version

Description

Gaussian Gram matrix and its double-centered version

Usage

hsic_gram_matrix(x, width)

Arguments

Value

list(K = Gram matrix, Kc = centered Gram matrix)

Median-heuristic kernel width for HSIC

Description

Uses only the first 100 points (not a random subsample) to keep the O(n^2) pairwise-distance computation cheap, matching the upstream implementation exactly.

Usage

hsic_kernel_width(x)

Arguments

Value

kernel width (scalar)

HSIC independence test with gamma approximation

Description

Faithful port of hsic_test_gamma() (hsic.py). O(n^2) in the sample size because it forms the full n x n Gram matrices; not recommended for n in the thousands.

Usage

hsic_test_gamma(X, Y)

Arguments

Details

The gamma-approximation variance estimator is only defined for n >= 6 (its closed form divides by (n-1)(n-2)(n-3)); below that this errors instead of silently returning a NaN p-value that would otherwise propagate into NA-valued rejection decisions in callers. A constant input (zero variance) carries no dependence information, so it is treated as trivially independent (p = 1) rather than routed through the degenerate kernel-width / centered-Gram-matrix computation.

Value

list(stat = HSIC test statistic, p = gamma-approximated p-value)

Select lambda by information criterion

Description

Select lambda by information criterion

Usage

ic_glmnet(glmnet_model)

Arguments

Value

list with lambda_AIC_best, lambda_BIC_best, idx_AIC_best, idx_BIC_best, ic_table

Pivoted incomplete Cholesky decomposition of a Gaussian kernel matrix, as used by f_correlation()

Description

Same greedy pivoted algorithm as incomplete_cholesky_gauss() (search_causal_order.r), but with the stopping rule used by the upstream ⁠_f_correlation.py⁠: continue while the sum of the remaining diagonal residuals exceeds tol (no rank cap), rather than stopping once the single largest residual drops below a fixed tolerance.

Usage

incomplete_cholesky_fcorr(x, sigma, tol)

Arguments

Value

n x d matrix G with ⁠tcrossprod(G) ~= K⁠

Gaussian kernel via pivoted incomplete Cholesky

Description

Low-rank approximation G (n x d) with ⁠tcrossprod(G) ~= K⁠, where K[i, j] = exp(-(x[i] - x[j])^2 / (2*sigma^2)), via the greedy pivoted incomplete Cholesky of Fine & Scheinberg (2001) / Bach & Jordan (2002, kernel-ICA). Each step picks the index with the largest diagonal residual and computes only the kernel column for that pivot, so the full n x n Gram matrix is never formed. The Gaussian kernel's diagonal is 1, so residuals start at 1 without evaluating K.

Usage

incomplete_cholesky_gauss(
  x,
  sigma,
  tol = 1e-04,
  max_rank = min(length(x), 200)
)

Arguments

Value

n x d matrix G with ⁠tcrossprod(G) ~= K⁠

Whether the OLS (or MLHSICR) residual of xi on xj_list is independent of every xj

Description

Faithful port of rcd.py's is_independent() helper (original 226-252 lines).

Usage

is_independent_of_resid(
  Y,
  xi,
  xj_list,
  MLHSICR,
  independence,
  ind_alpha,
  ind_corr
)

Arguments

Value

TRUE if independent

Kernel-based mutual information: core

Description

The target quantity is the difference of logdets of 2n x 2n matrices, but via the block structure and the Schur complement it can be computed equivalently using only an n x n Cholesky decomposition: ⁠MI = -1/2 * (logdet(tmp2^2 - K2 K1 tmp1^-2 K1 K2) - logdet(tmp2^2))⁠

Usage

kernel_mi_core(E1, x2, kappa, sigma)

Arguments

Value

Mutual information

Kernel-based mutual information: low-rank core

Description

Low-rank counterpart of kernel_mi_core(). With ⁠K2 ~= G2 %*% t(G2)⁠, both S = tmp2^2 - t(W) %*% W and tmp2^2 reduce to c0^2*I + G2 %*% C %*% t(G2) for a d2 x d2 matrix C, via the Woodbury identity on W and the matrix determinant lemma on the resulting rank-d2 update. The n*log(c0) terms in logdet(S) and logdet(tmp2^2) cancel in their difference, so only d2 x d2 matrices remain.

Usage

kernel_mi_core_lowrank(prep1, x2, kappa, sigma)

Arguments

Value

Mutual information

Kernel-based mutual information: precomputation for variable 1

Description

Computes the matrix ⁠E1 = tmp1^-1 K1⁠ (tmp1 = K1 + n*kappa/2 * I) used in kernel_mi_core(). It only needs to be called once per candidate variable, avoiding per-pair recomputation.

Usage

kernel_mi_prepare(x, kappa, sigma)

Arguments

Value

Matrix E1 (n x n)

Kernel-based mutual information: low-rank precomputation for variable 1

Description

Low-rank counterpart of kernel_mi_prepare(). Via the Woodbury identity, ⁠E1 = tmp1^-1 K1⁠ collapses to G1 %*% M1^-1 %*% t(G1) where M1 = c0*I + t(G1) %*% G1, so the n x n matrix E1 never needs to be formed; only the d x d inverse M1^-1 does.

Usage

kernel_mi_prepare_lowrank(x, kappa, sigma)

Arguments

Value

list(G, Minv, A, c0) describing E1 in factored form

Kurtosis (divided by n; excess kurtosis, which is 0 for a normal distribution)

Description

Kurtosis (divided by n; excess kurtosis, which is 0 for a normal distribution)

Usage

kurtosis_pop(x)

Arguments

Scale factor used to make the (otherwise fixed, absolute) lambda search grids adapt to the response's natural scale.

Description

lasso_lambda_seq / ridge_lambda_seq are fixed absolute grids. Because glmnet's standardize = TRUE only standardizes the predictors (not the response y), the penalty strength needed for meaningful shrinkage scales with the magnitude of y. Without this scaling, multiplying the whole input data by a constant changes which edges the default adaptive_lasso + lambda = "BIC"/"AIC" pipeline selects, even though the underlying relationships are identical up to that constant.

Usage

lambda_scale_factor(y)

Arguments

Value

a positive scale factor

Map lavaan fitMeasures() to the semopy-style column set

Description

Map lavaan fitMeasures() to the semopy-style column set

Usage

lavaan_fit_measures_to_df(fit)

Convert the doubled parameter vector (length 2dd) into a d x d matrix W

Description

Mirrors Python's ⁠w[:d*d].reshape([d, d])⁠ (row-major reshape), so w and the reverse operation lim_flatten_rowmajor() must stay consistent.

Usage

lim_adj(w, d)

BIC-penalized DAG score (sign-flipped, i.e. a "loss") used during the local search phase. W is interpreted as a 0/1 skeleton in the i -> j orientation (parents of j are the nonzero rows of column j).

Description

BIC-penalized DAG score (sign-flipped, i.e. a "loss") used during the local search phase. W is interpreted as a 0/1 skeleton in the i -> j orientation (parents of j are the nonzero rows of column j).

Usage

lim_bic_loss(W, X, is_continuous)

Flatten a matrix in row-major (numpy "C") order

Description

Flatten a matrix in row-major (numpy "C") order

Usage

lim_flatten_rowmajor(M)

Global (NOTEARS-style) optimization phase: augmented-Lagrangian outer loop with L-BFGS-B inner solves. Returns the thresholded W (i -> j orientation).

Description

Global (NOTEARS-style) optimization phase: augmented-Lagrangian outer loop with L-BFGS-B inner solves. Returns the thresholded W (i -> j orientation).

Usage

lim_global_optimize(X, con, lambda1, max_iter, h_tol, rho_max, w_threshold)

Acyclicity constraint h(W) and the matrix power term used in its gradient (Yu et al. 2019 formulation, as used by the Python source)

Description

Acyclicity constraint h(W) and the matrix power term used in its gradient (Yu et al. 2019 formulation, as used by the Python source)

Usage

lim_h(W, d)

Local log-likelihood of continuous component i (super-Gaussian disturbance)

Description

Direct port of likelihood_i/log_p_super_gaussian/variance_i from ⁠lingam/utils/__init__.py⁠ (lines 866-939 of the Python source).

Usage

lim_likelihood_i(X, i, b_full, b0)

Details

A (near-)zero residual variance means component i is fit deterministically by its parents; -n * log(sqrt(var_i)) would then diverge to Inf and make this candidate model look artificially superior. Such candidates are instead scored with a large finite penalty so they lose every BIC comparison in lim_bic_loss()/lim_local_search() without introducing Inf/NaN into the running score.

Local search phase: direction reversal, pruning, and edge addition (Python source lines 324-397). Returns W in the i -> j orientation.

Description

Local search phase: direction reversal, pruning, and edge addition (Python source lines 324-397). Returns W in the i -> j orientation.

Usage

lim_local_search(W_est, X, con, is_continuous, d, h_tol)

Numerically stable log(1 + exp(M)) (elementwise)

Description

Numerically stable log(1 + exp(M)) (elementwise)

Usage

lim_log1pexp(M)

Numerically stable log(cosh(r)) (elementwise)

Description

Numerically stable log(cosh(r)) (elementwise)

Usage

lim_logcosh(r)

Mixed loss (logistic for discrete columns, log-cosh/Laplace for continuous columns) and its gradient. con is a length-d vector, 1 = continuous, 0 = discrete. W_dis_mask / W_con_mask are precomputed d x d 0/1 masks.

Description

Mixed loss (logistic for discrete columns, log-cosh/Laplace for continuous columns) and its gradient. con is a length-d vector, 1 = continuous, 0 = discrete. W_dis_mask / W_con_mask are precomputed d x d 0/1 masks.

Usage

lim_loss_mixed(W, X, con, W_dis_mask, W_con_mask)

Value and gradient of the augmented-Lagrangian objective for the doubled parameter vector w (length 2dd)

Description

Value and gradient of the augmented-Lagrangian objective for the doubled parameter vector w (length 2dd)

Usage

lim_obj_grad(w, d, X, con, W_dis_mask, W_con_mask, rho, alpha, lambda1)

Topological order (Kahn's algorithm) from a lingamr-convention adjacency matrix B (B[i, j] = j -> i, i.e. row i's nonzero columns are its parents)

Description

Topological order (Kahn's algorithm) from a lingamr-convention adjacency matrix B (B[i, j] = j -> i, i.e. row i's nonzero columns are its parents)

Usage

lim_topological_order(B)

Direct LiNGAM

Description

Direct LiNGAM

Usage

lingam_direct(
  X,
  prior_knowledge = NULL,
  apply_prior_knowledge_softly = FALSE,
  measure = "pwling",
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Value

A LingamResult object (list) containing the following elements:

• adjacency_matrix: adjacency matrix B (n_features x n_features). Convention: B[i, j] is the causal coefficient from variable j to variable i (j -> i). Zero elements indicate no causal relationship.

• causal_order: estimated causal order (integer vector of 1-based indices). Earlier elements are more upstream (closer to exogenous variables).

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

# OLS (no additional packages required)
result <- lingam_direct(LiNGAM_sample_1000$data, reg_method = "ols")
round(result$adjacency_matrix, 3)


# LASSO (requires glmnet)
result_lasso <- lingam_direct(LiNGAM_sample_1000$data)
round(result_lasso$adjacency_matrix, 3)

Bootstrap for Direct LiNGAM

Description

Bootstrap for Direct LiNGAM

Usage

lingam_direct_bootstrap(
  X,
  n_sampling,
  prior_knowledge = NULL,
  apply_prior_knowledge_softly = FALSE,
  measure = "pwling",
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols",
  seed = NULL,
  verbose = TRUE,
  parallel = FALSE,
  n_cores = NULL,
  compute_total_effects = TRUE
)

Arguments

Details

When parallel = TRUE is specified, iterations are distributed across a socket cluster created by parallel::makePSOCKcluster(). The cluster is always released via on.exit(), whether the process finishes normally or an error occurs.

On iteration failures: each bootstrap iteration is run inside a tryCatch(). If an iteration errors (e.g. a resample produces near-singular columns), a warning identifying the failed iteration is issued and that iteration is excluded from the result instead of aborting the entire run. The returned BootstrapResult reflects however many iterations actually succeeded; an error is raised only if every iteration fails.

On reproducibility: During parallel execution, L'Ecuyer parallel random number streams via parallel::clusterSetRNGStream() are used. Results are reproducible given the same seed and same n_cores, but they do not numerically match the results of sequential execution (parallel = FALSE). If you need results that exactly match the sequential version, use parallel = FALSE.

Value

BootstrapResult (list)

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

# Fast example with OLS
bs <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 10L,
  reg_method = "ols",
  seed = 42
)
get_probabilities(bs)


# With LASSO (requires glmnet)
bs_lasso <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L,
  seed = 42
)

# Parallel execution on 2 cores
bs_par <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L,
  seed = 42,
  parallel = TRUE,
  n_cores = 2L
)

High-Dimensional Direct LiNGAM

Description

A variant of Direct LiNGAM for high-dimensional data (large p, or p > n). Causal order search is based on moment statistics of non-Gaussianity rather than pairwise independence measures, which is considerably faster for many variables. Unlike lingam_direct(), the algorithm is deterministic (no random restarts).

Usage

lingam_high_dim(X, J = 3L, K = 4L, alpha = 0.5, estimate_adj_mat = TRUE)

Arguments

Details

When n_samples <= n_features, the adjacency matrix cannot be estimated with the usual BIC-based Adaptive LASSO, so this function falls back to a cross-validated LASSO (glmnet::cv.glmnet) after emitting a warning. The upstream Python implementation uses LassoLarsCV for this fallback; cv.glmnet follows the same cross-validation design but is not numerically identical (different solver: coordinate descent vs. LARS).

The pruning statistic used during causal-order search (Wang & Drton 2020) is the minimum of a conditional non-Gaussianity statistic over every size-appropriate conditioning subset. The upstream Python implementation (cdt15/lingam) has a return statement mis-indented inside the loop over conditioning subsets, causing it to only ever evaluate the first subset. This R implementation intentionally does not replicate that bug and evaluates every subset, so causal_order and adjacency_matrix are not numerically identical to the upstream Python package.

Value

A LingamResult object (list), the same class returned by lingam_direct(), containing:

• adjacency_matrix: adjacency matrix B (n_features x n_features). Convention: B[i, j] is the causal coefficient from variable j to variable i (j -> i). Zero elements indicate no causal relationship. All-NA when estimate_adj_mat = FALSE.

• causal_order: estimated causal order (integer vector of 1-based indices). Earlier elements are more upstream.

References

Wang, Y. S. and Drton, M. (2020). High-dimensional causal discovery under non-Gaussianity. Biometrika, 107(1), 41-59.

Examples

sample <- generate_lingam_sample_6(n = 300, seed = 1)
result <- lingam_high_dim(sample$data)
result$causal_order
round(result$adjacency_matrix, 3)


if (requireNamespace("glmnet", quietly = TRUE)) {
  # n <= p: falls back to cross-validated LASSO with a warning
  wide_sample <- generate_lingam_large_sample(p = 30, n = 25, seed = 1)
  wide_result <- lingam_high_dim(wide_sample$data)
  wide_result$causal_order
}

LiM: LiNGAM for Mixed Data

Description

Estimates a causal structure from data containing a mixture of continuous and binary (0/1) discrete variables, following Zeng et al. (2022). The method combines a NOTEARS-style continuous optimization (global phase) with a combinatorial local search over edge directions, pruning, and edge addition.

Usage

lingam_lim(
  X,
  is_continuous,
  lambda1 = 0.1,
  max_iter = 150L,
  h_tol = 1e-08,
  rho_max = 1e+16,
  w_threshold = 0.1,
  only_global = FALSE
)

Arguments

Details

Only binary (0/1) discrete variables are supported; count/Poisson-type discrete variables (the Python source's loss_type = "poisson" path) are not implemented.

The Python implementation's adjacency_matrix_ uses the opposite convention (W[i, j] = i -> j). This R implementation transposes the internal result so that adjacency_matrix follows the lingamr convention (B[i, j] = j -> i), consistent with lingam_direct().

In the local search phase, edges that are reversed or newly added are assigned a weight of exactly 1 rather than a re-estimated coefficient (this matches the original Python implementation). As a result, non-zero entries of adjacency_matrix are a mix of global-phase estimated coefficients and local-phase placeholder weights of 1.

The local search's BIC score for discrete variables with parents uses R's glm(family = binomial()), an unregularized maximum-likelihood fit. The Python implementation uses scikit-learn's LogisticRegression, which applies L2 regularization by default; consequently, numeric results will not exactly match the Python implementation.

Value

A LiMResult object (list) containing the following elements:

• adjacency_matrix: adjacency matrix B (n_features x n_features). Convention: B[i, j] is the causal coefficient from variable j to variable i (j -> i), i.e. the same convention as lingam_direct(). Zero elements indicate no causal relationship.

• causal_order: estimated causal order (integer vector of 1-based indices), derived from adjacency_matrix via a topological sort. NA (with a warning) if the estimated matrix is not acyclic.

• is_continuous: the input is_continuous vector, stored for reference.

References

Zeng Y, Shimizu S, Matsui H, Sun F. Causal discovery for linear mixed data. In: Proceedings of the First Conference on Causal Learning and Reasoning (CLeaR 2022). PMLR 177, pp. 994-1009, 2022.

Examples

# Reproducibility requires set.seed(), since the optimization starts from
# a random initial point.
set.seed(1)
dat <- generate_lim_sample(n = 300)
result <- lingam_lim(dat$data, is_continuous = dat$is_continuous)
print(result)

Multi-Group Direct LiNGAM

Description

Jointly estimates a Direct LiNGAM model from multiple datasets ("groups") that are assumed to share a common causal order but may have different structural coefficients (Shimizu 2012).

Usage

lingam_multi_group(
  X_list,
  prior_knowledge = NULL,
  apply_prior_knowledge_softly = FALSE,
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Details

Unlike lingam_direct(), this function has no measure argument: the multi-group causal-order search only supports the pwling (pairwise likelihood / entropy approximation) objective, matching the upstream Python MultiGroupDirectLiNGAM, which does not offer a kernel-based multi-group search.

For downstream analysis of a single group (total causal effects, independence tests of residuals, plotting), extract that group as a plain LingamResult with get_group_result() and pass it to the existing single-group functions (estimate_total_effect(), estimate_all_total_effects(), get_error_independence_p_values(), plot_adjacency(), autoplot(), tidy()); no multi-group-specific wrappers are provided for these.

Value

A MultiGroupLingamResult object (list) containing:

• adjacency_matrices: a named list of adjacency matrices, one per group (name = group name). Each follows the same convention as lingam_direct(): B[i, j] is the causal coefficient from variable j to variable i (j -> i).

• causal_order: the causal order shared by all groups (integer vector of 1-based indices).

References

S. Shimizu. Joint estimation of linear non-Gaussian acyclic models. Neurocomputing, 81: 104-107, 2012.

Examples

mg <- generate_multi_group_sample()
res <- lingam_multi_group(mg$data_list, reg_method = "ols")
print(res)

# Analyze group 1 with the existing single-group tooling
g1 <- get_group_result(res, 1)
estimate_all_total_effects(mg$data_list[[1]], g1, method = "ols")

Bootstrap for Multi-Group Direct LiNGAM

Description

Bootstrap for Multi-Group Direct LiNGAM

Usage

lingam_multi_group_bootstrap(
  X_list,
  n_sampling,
  prior_knowledge = NULL,
  apply_prior_knowledge_softly = FALSE,
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols",
  seed = NULL,
  verbose = TRUE,
  parallel = FALSE,
  n_cores = NULL,
  compute_total_effects = TRUE
)

Arguments

Details

Each element of the returned list is a regular BootstrapResult, so the existing single-group bootstrap functions (get_probabilities(), get_causal_direction_counts(), get_directed_acyclic_graph_counts(), get_total_causal_effects(), get_causal_order_stability(), plot_bootstrap_probabilities(), tidy()) all work by extracting a group with result[[group_name]] or result[[i]], mirroring the upstream Python API (which returns a list of BootstrapResult per group).

Total effects use path products, not regression. Each bootstrap iteration's total-effect matrix is the sum of path-coefficient products over the DAG defined by that iteration's adjacency matrix (matching the upstream Python MultiGroup bootstrap), which is a different method from lingam_direct_bootstrap()'s regression-based estimate_all_total_effects().

On iteration failures: as in lingam_direct_bootstrap(), each iteration is wrapped in tryCatch(); a failing iteration is skipped with a warning, and only if every iteration fails is an error raised.

On reproducibility: same policy as lingam_direct_bootstrap(). During parallel execution, L'Ecuyer parallel random number streams via parallel::clusterSetRNGStream() are used. Results are reproducible given the same seed and same n_cores, but they do not numerically match the results of sequential execution (parallel = FALSE). If you need results that exactly match the sequential version, use parallel = FALSE.

Value

A MultiGroupBootstrapResult: a named list (one element per group) of BootstrapResult objects (see lingam_direct_bootstrap()), with class "MultiGroupBootstrapResult".

Examples

mg <- generate_multi_group_sample()

bs <- lingam_multi_group_bootstrap(mg$data_list,
  n_sampling = 10L,
  reg_method = "ols",
  seed = 42
)
get_probabilities(bs[[1]])


bs_par <- lingam_multi_group_bootstrap(mg$data_list,
  n_sampling = 30L,
  reg_method = "ols",
  seed = 42,
  parallel = TRUE,
  n_cores = 2L
)

Bottom-Up ParceLiNGAM

Description

A causal ordering method robust against latent confounders. Unlike lingam_direct(), which always returns a full causal order, this algorithm searches from the sink (most downstream) side and stops as soon as an independence test is rejected. Variables it could not order are returned together as a single "unresolved block" (suspected to share a latent confounder), and the corresponding entries of the adjacency matrix are set to NA rather than estimated.

Usage

lingam_parce(
  X,
  alpha = 0.1,
  prior_knowledge = NULL,
  independence = "hsic",
  ind_corr = 0.5,
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols"
)

Arguments

Details

Because HSIC forms full n x n Gram matrices, it is O(n^2) per test; avoid very large n (beyond a few thousand) with independence = "hsic".

independence = "fcorr" rejects based on the raw F-correlation value (ind_corr), not a p-value, so it is not directly comparable to alpha.

get_error_independence_p_values_parce() uses the HSIC test rather than the correlation-based test used by get_error_independence_p_values() for LingamResult objects.

lingam_parce_bootstrap() treats NA (unresolved) edges as absent when aggregating, and does not support get_causal_order_stability() (see its documentation for details). This function does not expose a regressor or bw_method argument, unlike the upstream Python implementation.

Value

A ParceLingamResult object (list) containing:

• adjacency_matrix: adjacency matrix B (n_features x n_features). Convention: B[i, j] is the causal coefficient from variable j to variable i (j -> i), same as lingam_direct(). Entries between two variables that ended up in the same unresolved block are NA.

• causal_order: a list of integer vectors. Elements of length 1 are variables with a fully resolved position; an element of length > 1 (at most one, always first) is the unresolved block. Earlier elements are more upstream.

• p_values: independence-test p-values (or F-correlation values, for independence = "fcorr") for each step that successfully placed a variable, in the order variables were placed (diagnostic only).

• independence: the independence measure used.

References

Tashiro, T., Shimizu, S., Hyvarinen, A., and Washio, T. (2014). ParceLiNGAM: a causal ordering method robust against latent confounders. Neural Computation, 26(1), 57-83.

Examples

confounded <- generate_parce_sample(n = 500, seed = 1)

result <- lingam_parce(confounded$data, reg_method = "ols")
print(result)

# The variable pair sharing the latent confounder is left unresolved (NA)
result$adjacency_matrix[confounded$confounded_pair, confounded$confounded_pair]

# Total effect estimation warns and returns NA for confounded variables
estimate_total_effect_parce(confounded$data, result,
  from_index = confounded$confounded_pair[1], to_index = 1
)

Bootstrap for Bottom-Up ParceLiNGAM

Description

Bootstrap for Bottom-Up ParceLiNGAM

Usage

lingam_parce_bootstrap(
  X,
  n_sampling,
  prior_knowledge = NULL,
  alpha = 0.1,
  independence = "hsic",
  ind_corr = 0.5,
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols",
  seed = NULL,
  verbose = TRUE,
  parallel = FALSE,
  n_cores = NULL,
  compute_total_effects = TRUE
)

Arguments

Details

Total effects are path sums, not regression estimates. Each iteration's total-effect matrix is built from calculate_total_effect() (summing products of adjacency-matrix coefficients along every directed path), matching the upstream Python implementation's bootstrap method (estimate_total_effect2). If a variable's row in the adjacency matrix contains NA (it is part of an unresolved block), all of its outgoing total effects are set to NA for that iteration, since its causal parents cannot be identified.

NA (unresolved) edges are treated as absent when aggregating. Both the adjacency matrix and the total-effect matrix have NA replaced by 0 before being stored in the returned BootstrapResult, matching the numpy comparison semantics used by the upstream implementation (where np.abs(nan) > threshold evaluates to FALSE). This means, for example, get_probabilities() reports the confounded pair's edge probability as the fraction of resamples in which the order happened to resolve, not as NA.

causal_orders is not populated (unlike lingam_direct_bootstrap()): ParceLiNGAM's causal order can include an unresolved block, which does not fit the fixed-length integer-vector format causal_orders requires. As a result, get_causal_order_stability() cannot be used with a BootstrapResult returned by this function.

Value

A BootstrapResult (list); see lingam_direct_bootstrap() for the query helpers that operate on it (get_probabilities(), get_causal_direction_counts(), get_directed_acyclic_graph_counts(), get_total_causal_effects()).

Examples

confounded <- generate_parce_sample(n = 500, seed = 1)

bs <- lingam_parce_bootstrap(confounded$data,
  n_sampling = 10L,
  reg_method = "ols",
  seed = 42
)
get_probabilities(bs)

RCD (Repetitive Causal Discovery)

Description

A causal discovery method robust against latent confounders. Unlike lingam_direct() or lingam_parce(), RCD does not attempt to recover a full or partial causal order. Instead, it repeatedly extracts each variable's ancestor set by scanning variable subsets of increasing size (extract_ancestors()), narrows each ancestor set down to direct parents (extract_parents()), and finally tests remaining parent-free pairs for a shared latent confounder (extract_vars_sharing_confounders()). Pairs found to share a latent confounder are marked NA in the adjacency matrix rather than estimated.

Usage

lingam_rcd(
  X,
  max_explanatory_num = 2L,
  cor_alpha = 0.01,
  ind_alpha = 0.01,
  shapiro_alpha = 0.01,
  MLHSICR = FALSE,
  independence = "hsic",
  ind_corr = 0.5
)

Arguments

Details

The algorithm has three stages: (1) extract_ancestors() grows each variable's ancestor set by repeatedly scanning variable subsets; (2) extract_parents() narrows ancestor sets down to direct parents; (3) extract_vars_sharing_confounders() tests remaining parent-free pairs for a shared latent confounder. NA entries in adjacency_matrix mean the corresponding pair is suspected to share a latent confounder, not that no relationship was estimated.

max_explanatory_num controls both statistical power and computational cost: stage 1 scans choose(n_features, k) subsets for each subset size k up to max_explanatory_num + 1, and each subset requires several HSIC tests (each O(n^2) in the sample size when independence = "hsic"). Cost grows quickly with both the number of variables and n.

MLHSICR = TRUE replaces the OLS residual in the independence check with a residual obtained by directly minimizing the sum of HSIC statistics between the residual and each explanatory variable via numerical optimization (stats::optim(method = "L-BFGS-B")). This can recover independence in cases where OLS cannot, but requires re-optimizing for every candidate subset where the OLS residual fails, and is therefore substantially slower.

The Shapiro-Wilk test (stats::shapiro.test()) used for the non-Gaussianity check is limited to n between 3 and 5000. For n above 5000, a deterministic evenly-spaced subsample of 5000 observations is tested instead (same policy as test_residual_normality()), so results remain reproducible without touching the RNG state. This subsampling has no effect when n <= 5000.

This function does not expose a bw_method argument (kernel widths are always the median heuristic; see hsic_kernel_width()), unlike some upstream implementations. lingam_rcd_bootstrap() does not support get_causal_order_stability(), since RCD has no causal order.

Value

An RCDResult object (list) containing:

• adjacency_matrix: adjacency matrix B (n_features x n_features). Convention: B[i, j] is the causal coefficient from variable j to variable i (j -> i), same as lingam_direct(). Entries between two variables found to share a latent confounder are NA.

• ancestors_list: a list of length n_features; element i is the sorted integer vector of variables found to be ancestors of variable i (possibly empty). Unlike lingam_parce(), there is no causal_order: RCD estimates ancestor relations directly rather than a total or partial order.

References

Maeda, T. N. and Shimizu, S. (2020). RCD: Repetitive causal discovery of linear non-Gaussian acyclic models with latent confounders. AISTATS 2020, PMLR 108: 735-745.

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)

result <- lingam_rcd(confounded$data)
print(result)

# The variable pair sharing the latent confounder is left NA
result$adjacency_matrix[confounded$confounded_pair, confounded$confounded_pair]

# Total effect estimation warns and returns NA for confounded variables
estimate_total_effect_rcd(confounded$data, result,
  from_index = confounded$confounded_pair[1], to_index = 1
)

Bootstrap for RCD

Description

Bootstrap for RCD

Usage

lingam_rcd_bootstrap(
  X,
  n_sampling,
  max_explanatory_num = 2L,
  cor_alpha = 0.01,
  ind_alpha = 0.01,
  shapiro_alpha = 0.01,
  MLHSICR = FALSE,
  independence = "hsic",
  ind_corr = 0.5,
  seed = NULL,
  verbose = TRUE,
  parallel = FALSE,
  n_cores = NULL,
  compute_total_effects = TRUE
)

Arguments

Details

Total effects are computed only for ancestor pairs, driven by each iteration's ancestors_list (unlike lingam_direct_bootstrap() and lingam_parce_bootstrap(), which loop over all variable pairs). For a pair ⁠(from, to)⁠ with from in to's ancestor set, the total effect is obtained via calculate_total_effect() (summing products of adjacency-matrix coefficients along every directed path). If from's row in the adjacency matrix contains NA (it shares a latent confounder with some other variable), the effect is set to NA for that iteration.

NA (confounded) edges are treated as absent when aggregating. Both the adjacency matrix and the total-effect matrix have NA replaced by 0 before being stored in the returned BootstrapResult, matching the policy used by lingam_parce_bootstrap().

causal_orders is not populated: RCD has no causal order (see lingam_rcd()). As a result, get_causal_order_stability() cannot be used with a BootstrapResult returned by this function.

RCD's fit step is HSIC-heavy and can be slow per iteration, especially with MLHSICR = TRUE; keep n_sampling modest in examples.

Value

A BootstrapResult (list); see lingam_direct_bootstrap() for the query helpers that operate on it (get_probabilities(), get_causal_direction_counts(), get_directed_acyclic_graph_counts(), get_total_causal_effects()).

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)

bs <- lingam_rcd_bootstrap(confounded$data,
  n_sampling = 5L,
  seed = 42
)
get_probabilities(bs)

Compute residuals (error terms) of a LiNGAM model

Description

After validating the inputs (that lingam_result is a LingamResult, that X is numeric, and that the dimensions match), returns ⁠E = X - X B^T⁠. Shared by the residual-based diagnostic functions.

Usage

lingam_residuals(X, lingam_result)

Arguments

Value

residual matrix (n_samples x n_features). Preserves the column names of X.

VAR-LiNGAM for time series causal discovery

Description

Fits a vector autoregressive (VAR) model to time series data and applies Direct LiNGAM to the residuals to recover the instantaneous (lag-0) causal structure. The lagged causal matrices are then derived from the VAR coefficients and the instantaneous structure.

Usage

lingam_var(
  X,
  lags = 1L,
  criterion = "bic",
  measure = "pwling",
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols",
  prune = TRUE
)

Arguments

Details

The model is ⁠X_t = B0 X_t + sum_{k=1}^{p} B_k X_{t-k} + e_t⁠, where B0 is the instantaneous effect matrix (strictly acyclic) and e_t are mutually independent non-Gaussian errors. VAR coefficients M_k are estimated by ordinary least squares (no intercept); residuals ⁠e_t = X_t - sum M_k X_{t-k}⁠ are passed to lingam_direct() to obtain B0, and the lagged matrices follow ⁠B_k = (I - B0) M_k⁠.

Value

A VARLiNGAMResult object (list) containing:

• adjacency_matrices: array (1 + lags, n_features, n_features). The first slice ⁠[1, , ]⁠ is the instantaneous matrix B0; slice ⁠[k + 1, , ]⁠ is the lagged matrix B_k for lag k (k = 1..lags). Convention: B[i, j] is the effect from variable j to variable i.

• causal_order: estimated causal order of the instantaneous structure (1-based indices).

• residuals: VAR residuals (n_samples - lags, n_features).

• lags: the lag order actually used.

References

Hyvärinen, A., Zhang, K., Shimizu, S., & Hoyer, P. O. (2010). Estimation of a structural vector autoregression model using non-Gaussianity. Journal of Machine Learning Research, 11, 1709-1731. Ported from the Python implementation cdt15/lingam (https://github.com/cdt15/lingam). See also the VARLiNGAM R code of Moneta et al. (https://sites.google.com/site/dorisentner/publications/VARLiNGAM).

Examples

sample <- generate_varlingam_sample(n = 500, seed = 42)

# OLS instantaneous structure without pruning (no extra packages required)
model <- lingam_var(sample$data, lags = 1, reg_method = "ols", prune = FALSE)
round(model$adjacency_matrices[1, , ], 2)  # instantaneous B0

Bootstrap for VAR-LiNGAM

Description

Evaluates the statistical reliability of the estimated time-series DAG by resampling. Unlike the i.i.d. row resampling used for Direct LiNGAM, this uses a residual bootstrap: the VAR is fitted once on the original data, the residuals are resampled with replacement, and a new series is rebuilt by the VAR recursion before re-estimating VAR-LiNGAM on it (this preserves the autoregressive structure). Port of the Python reference VARLiNGAM.bootstrap.

Usage

lingam_var_bootstrap(
  X,
  n_sampling,
  lags = 1L,
  criterion = "bic",
  measure = "pwling",
  reg_method = "adaptive_lasso",
  lambda = "BIC",
  init_method = "ols",
  prune = TRUE,
  seed = NULL,
  verbose = TRUE,
  parallel = FALSE,
  n_cores = NULL
)

Arguments

Details

Reproducibility follows the same rules as lingam_direct_bootstrap(): with parallel = TRUE, L'Ecuyer streams via parallel::clusterSetRNGStream() make results reproducible for a given seed and n_cores, but they do not match the sequential (parallel = FALSE) results.

On iteration failures: as in lingam_direct_bootstrap(), each iteration runs inside a tryCatch(); a failing iteration is reported as a warning and excluded from the result instead of aborting the run. An error is raised only if every iteration fails.

Value

a VARBootstrapResult object.

Examples

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

# Fast example: OLS instantaneous structure, no pruning (no glmnet needed)
bs <- lingam_var_bootstrap(s$data,
  n_sampling = 10L, lags = 1, criterion = NULL,
  reg_method = "ols", prune = FALSE, seed = 1, verbose = FALSE
)
get_var_probabilities(bs)

Create noise generation function

Description

Internal helper to create a noise function for the specified distribution.

Usage

make_noise_fn(noise_dist)

Arguments

Value

function(n) that generates n random numbers

Create a prior knowledge matrix

Description

Create a prior knowledge matrix

Usage

make_prior_knowledge(
  n_variables,
  exogenous_variables = NULL,
  sink_variables = NULL,
  paths = NULL,
  no_paths = NULL,
  labels = NULL
)

Arguments

Value

Prior knowledge matrix (n_variables x n_variables) -1: unknown, 0: no path, 1: path exists

Examples

# Specify by index
pk <- make_prior_knowledge(6, exogenous_variables = c(4))

# Specify by variable name
pk <- make_prior_knowledge(6,
  exogenous_variables = "x3",
  sink_variables = c("x1", "x4"),
  paths = list(c("x3", "x0"), c("x3", "x2")),
  no_paths = list(c("x5", "x2")),
  labels = c("x0", "x1", "x2", "x3", "x4", "x5")
)

HSIC-sum-minimizing ("MLHSICR") regression

Description

Regresses Y[, xi] on Y[, xj_list] (no intercept) by minimizing the sum of the empirical HSIC statistics between the residual and each explanatory variable, instead of ordinary least squares. Used as a fallback when the OLS residual is not independent of the explanatory variables.

Usage

mlhsicr_regression(Y, xi, xj_list)

Arguments

Details

The kernel width used to build the residual's Gram matrix is itself a linear combination of the explanatory variables' kernel widths (faithful to the upstream implementation; see original 207-208 lines). This is an unusual design choice but is deliberately preserved as-is.

Value

list(resid = residual vector, coef = coefficient vector)

Compute a total-effect matrix from an adjacency matrix via path products

Description

Unlike estimate_all_total_effects() (regression-based), this sums path products over the DAG defined by B, matching the upstream Python MultiGroup bootstrap's calculate_total_effect().

Usage

multi_group_total_effect_matrix(B, causal_order)

Arguments

Value

total-effect matrix (n_features x n_features)

Kernel-based mutual information

Description

Dispatches to the incomplete-Cholesky low-rank path for n above the low-rank threshold (matching the kappa/sigma switch in search_causal_order_kernel()); below the threshold it calls the exact path unchanged.

Usage

mutual_information_kernel(x1, x2, param)

Arguments

Value

Mutual information

Evaluate the independence of a residual against a fixed set of predictors

Description

Used for the length(Uc) == 1 special case in find_exo_vec(), where there is nothing left to compare against.

Usage

parce_eval_independence(X, predictors, R, independence)

Arguments

Value

evaluation value (Fisher-combined p-value for hsic, max F-correlation for fcorr)

Find the position (rank) of a variable within a ParceLiNGAM causal order

Description

All members of the unresolved block (if any) share the same rank (1).

Usage

parce_order_rank(causal_order, idx)

Arguments

Value

integer rank, or NA if not found

Residual of X[, j] regressed on X[, xi_index] via the pseudo-inverse of the covariance matrix

Description

Residual of X[, j] regressed on X[, xi_index] via the pseudo-inverse of the covariance matrix

Usage

parce_residual(X, xi_index, j, Cov)

Arguments

Value

residual vector

Candidate sink variables under prior knowledge (ParceLiNGAM direction)

Description

Unlike search_candidate() (used by top-down DirectLiNGAM), ParceLiNGAM searches from the sink side, so the filter is simply "exclude variables that appear as the 'from' side of a known partial order" (they cannot be a sink because something is known to cause a variable through them... more precisely: a variable known to cause another remaining variable cannot itself be the next sink).

Usage

parce_search_candidate(U, partial_orders)

Arguments

Value

candidate index vector (subset of U, or U itself if no candidates remain)

Bottom-up causal order search

Description

Faithful port of ⁠_search_causal_order()⁠ (original 188-225 lines). Repeatedly finds the most sink-like remaining variable and appends it to the front of K_bttm (bottom-up, so more recently placed variables are more upstream). Stops as soon as a candidate is rejected by the independence test; the remaining undetermined variables are returned as U_res.

Usage

parce_search_causal_order(
  X,
  U,
  partial_orders,
  independence,
  thresh_p,
  ind_corr
)

Arguments

Value

list(K_bttm = integer vector, p_bttm = numeric vector, U_res = integer vector)

Moore-Penrose pseudo-inverse via SVD

Description

Internal replacement for numpy.linalg.pinv, used to solve the conditional regression coefficients from Gram-matrix submatrices in lingam_high_dim(). Implemented in base R to avoid a new dependency.

Usage

pinv(A, tol = max(dim(A)) * .Machine$double.eps)

Arguments

Value

the pseudo-inverse of A

Plot a causal graph from an adjacency matrix with DiagrammeR

Description

Plot a causal graph from an adjacency matrix with DiagrammeR

Usage

plot_adjacency(
  B,
  labels = NULL,
  threshold = 0,
  rankdir = "TB",
  title = "Estimated Causal Structure",
  shape = "circle",
  fillcolor = "lightyellow",
  bordercolor = "black",
  fontsize_node = 14,
  fontsize_edge = 10,
  edge_color = "gray40",
  edge_label_color = "red",
  true_B = NULL,
  color_tp = "forestgreen",
  color_fp = "firebrick",
  color_fn = "darkorange",
  debug = FALSE
)

Arguments

Value

A grViz object (when DiagrammeR is available)

Examples

if (requireNamespace("DiagrammeR", quietly = TRUE)) {
  LiNGAM_sample_1000 <- generate_lingam_sample_6()

  LiNGAM_sample_1000$true_adjacency |>
    plot_adjacency(title = "True Causal Structure")

  model <- LiNGAM_sample_1000$data |>
    lingam_direct(reg_method = "ols")

  model$adjacency_matrix |>
    plot_adjacency()

  
  # Compare with the true structure
  # (correct = green, false positive = red, missed = orange dashed)
  model$adjacency_matrix |>
    plot_adjacency(true_B = LiNGAM_sample_1000$true_adjacency)
  
}

Draw bootstrap probabilities with DiagrammeR

Description

Draw bootstrap probabilities with DiagrammeR

Usage

plot_bootstrap_probabilities(
  result,
  labels = NULL,
  min_causal_effect = NULL,
  min_probability = 0.5,
  rankdir = "TB",
  shape = "circle"
)

Arguments

Value

grViz object

Examples

if (requireNamespace("DiagrammeR", quietly = TRUE)) {
  LiNGAM_sample_1000 <- generate_lingam_sample_6()

  bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
    n_sampling = 30L, reg_method = "ols", seed = 42
  )
  plot_bootstrap_probabilities(bs_model)
}

plot QQ

Description

plot QQ

Usage

plot_residual_qq(X, lingam_result, ncol = 3, nrow = NULL)

Arguments

Value

A ggplot2::ggplot object with QQ plots of residuals.

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  # Load the sample data
  LiNGAM_sample_1000 <- generate_lingam_sample_6()

  # Run Direct LiNGAM
  result <- lingam_direct(LiNGAM_sample_1000$data, reg_method = "ols")

  plot_residual_qq(LiNGAM_sample_1000$data, result)
}

Q-Q plots of VAR-LiNGAM residuals

Description

Draws per-variable normal Q-Q plots of the residuals (analogous to the Moneta Gauss_Stats visual check). Deviations from the reference line indicate non-Gaussianity, which supports the LiNGAM assumption. Requires ggplot2.

Usage

plot_varlingam_residual_qq(
  result,
  on = c("innovations", "var"),
  ncol = 3,
  nrow = NULL
)

Arguments

Value

a ggplot object

References

Analogous to the residual visual check (Gauss_Stats) in the VARLiNGAM R code of Moneta, A., Entner, D., Hoyer, P. O., & Coad, A. (2013), Oxford Bulletin of Economics and Statistics, 75(5), 705-730. https://sites.google.com/site/dorisentner/publications/VARLiNGAM

Examples

s <- generate_varlingam_sample(n = 1000, seed = 42)
m <- lingam_var(s$data, lags = 1, reg_method = "ols", prune = FALSE)

plot_varlingam_residual_qq(m)

Adaptive LASSO with CV-selected lambda (n <= p route)

Description

Replicates upstream ⁠_predict_adaptive_lasso⁠ (StandardScaler + OLS weights + LassoLarsCV), substituting glmnet::cv.glmnet(alpha = 1) for LassoLarsCV (see lingam_high_dim() Details for the rationale).

Usage

predict_adaptive_lasso_cv(X, predictors, target, gamma = 1)

Arguments

Value

coefficient vector, same length and order as predictors

Display the contents of a BootstrapResult

Description

Display the contents of a BootstrapResult

Usage

## S3 method for class 'BootstrapResult'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

bs_model <- lingam_direct_bootstrap(LiNGAM_sample_1000$data,
  n_sampling = 30L, reg_method = "ols", seed = 42
)

print(bs_model)

Print method for ImputationBootstrapResult

Description

Print method for ImputationBootstrapResult

Usage

## S3 method for class 'ImputationBootstrapResult'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

set.seed(1)
sample6 <- generate_lingam_sample_6(n = 300, seed = 1)
X <- sample6$data
X$x5[sample.int(nrow(X), size = 30)] <- NA

if (requireNamespace("mice", quietly = TRUE)) {
  res <- bootstrap_with_imputation(X,
    n_sampling = 5L, n_repeats = 3L, seed = 42, verbose = FALSE
  )
  print(res)
}

Print method for LiMResult

Description

Print method for LiMResult

Usage

## S3 method for class 'LiMResult'
print(x, digits = 3, ...)

Arguments

Value

The input object x, invisibly.

Examples

set.seed(1)
dat <- generate_lim_sample(n = 300)
result <- lingam_lim(dat$data, is_continuous = dat$is_continuous)
print(result)

Print method for LingamResult

Description

Print method for LingamResult

Usage

## S3 method for class 'LingamResult'
print(x, digits = 3, ...)

Arguments

Value

The input object x, invisibly.

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()
result <- lingam_direct(LiNGAM_sample_1000$data, reg_method = "ols")
print(result)

Print method for MultiGroupBootstrapResult

Description

Print method for MultiGroupBootstrapResult

Usage

## S3 method for class 'MultiGroupBootstrapResult'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

mg <- generate_multi_group_sample()
bs <- lingam_multi_group_bootstrap(mg$data_list,
  n_sampling = 10L, reg_method = "ols", seed = 42
)
print(bs)

Print method for MultiGroupLingamResult

Description

Print method for MultiGroupLingamResult

Usage

## S3 method for class 'MultiGroupLingamResult'
print(x, digits = 3, ...)

Arguments

Value

The input object x, invisibly.

Examples

mg <- generate_multi_group_sample()
res <- lingam_multi_group(mg$data_list, reg_method = "ols")
print(res)

Print method for ParceLingamResult

Description

Print method for ParceLingamResult

Usage

## S3 method for class 'ParceLingamResult'
print(x, digits = 3, ...)

Arguments

Value

The input object x, invisibly.

Examples

confounded <- generate_parce_sample(n = 300, seed = 42)
result <- lingam_parce(confounded$data, reg_method = "ols")
print(result)

Print method for RCDResult

Description

Print method for RCDResult

Usage

## S3 method for class 'RCDResult'
print(x, digits = 3, ...)

Arguments

Value

The input object x, invisibly.

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)
result <- lingam_rcd(confounded$data)
print(result)

Print a VARBootstrapResult

Description

Print a VARBootstrapResult

Usage

## S3 method for class 'VARBootstrapResult'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

s <- generate_varlingam_sample(n = 500, seed = 42)
bs <- lingam_var_bootstrap(s$data,
  n_sampling = 10L, lags = 1, criterion = NULL,
  reg_method = "ols", prune = FALSE, seed = 1, verbose = FALSE
)
print(bs)

Print method for VARLiNGAMResult

Description

Print method for VARLiNGAMResult

Usage

## S3 method for class 'VARLiNGAMResult'
print(x, digits = 3, ...)

Arguments

Value

The input object x, invisibly.

Examples

sample <- generate_varlingam_sample(n = 500, seed = 42)
model <- lingam_var(sample$data, lags = 1, reg_method = "ols", prune = FALSE)
print(model)

print method for causal_order_stability

Description

print method for causal_order_stability

Usage

## S3 method for class 'causal_order_stability'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

dat <- generate_lingam_sample_6()
bs <- lingam_direct_bootstrap(dat$data, n_sampling = 30L, reg_method = "ols", seed = 42)
print(get_causal_order_stability(bs, labels = names(dat$data)))

Print method for lingam_normality_test

Description

Print method for lingam_normality_test

Usage

## S3 method for class 'lingam_normality_test'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()
result <- lingam_direct(LiNGAM_sample_1000$data, reg_method = "ols")
print(test_residual_normality(LiNGAM_sample_1000$data, result))

print method for lingam_summary

Description

print method for lingam_summary

Usage

## S3 method for class 'lingam_summary'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()
model <- lingam_direct(LiNGAM_sample_1000$data, reg_method = "ols")
print(summary_lingam(LiNGAM_sample_1000$data, model))

Print method for var_stationarity

Description

Print method for var_stationarity

Usage

## S3 method for class 'var_stationarity'
print(x, ...)

Arguments

Value

The input object x, invisibly.

Examples

s <- generate_varlingam_sample(n = 1000, seed = 42)
m <- lingam_var(s$data, lags = 1, reg_method = "ols", prune = FALSE)
print(check_var_stationarity(m))

Prune VAR-LiNGAM adjacency matrices by adaptive LASSO

Description

Re-estimates the instantaneous matrix B0 and every lagged matrix B_k jointly, shrinking weak edges to zero. Port of the Python reference ⁠_pruning⁠. For each target variable, the predictors are its contemporaneous ancestors (those preceding it in causal_order) plus all variables at lags 1..lags; the coefficients are fitted by adaptive LASSO and written back into B.

Usage

prune_var_lingam(X, causal_order, lags, lambda = "BIC", init_method = "ols")

Arguments

Value

array (lags + 1, n_features, n_features); slice 1 is B0, slice k+1 is B_k

Pearson correlation test used throughout RCD

Description

Pearson correlation test used throughout RCD

Usage

rcd_is_correlated(a, b, cor_alpha)

Arguments

Value

TRUE if correlated (p-value < cor_alpha)

Independence judgment used throughout RCD (hsic or fcorr)

Description

Independence judgment used throughout RCD (hsic or fcorr)

Usage

rcd_is_independent(a, b, independence, ind_alpha, ind_corr)

Arguments

Value

TRUE if independent

Non-Gaussianity judgment (Shapiro-Wilk) for a set of columns

Description

When n > SHAPIRO_MAX_N, stats::shapiro.test() cannot be used directly (hard cap at 5000), so the deterministic evenly-spaced subsample from shapiro_subsample() (defined in R/get_error_independence_p_values.r) is tested instead, matching the behavior of test_residual_normality(). The deterministic thinning keeps results reproducible across calls and leaves the caller's RNG stream untouched.

Usage

rcd_is_non_gaussian_all(Y, cols, shapiro_alpha)

Arguments

Value

TRUE if all columns reject normality (p <= shapiro_alpha)

OLS residual and coefficients (with intercept)

Description

OLS residual and coefficients (with intercept)

Usage

rcd_ols_resid_coef(y, X_exog)

Arguments

Value

list(resid = residual vector, coef = coefficient vector excluding intercept)

Resolve a variable index or name to a 1-based integer index

Description

Resolve a variable index or name to a 1-based integer index

Usage

rcd_resolve_index(idx, arg_name, n_features, col_names)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

generics

glance(), tidy()

Residual (residual when xi is regressed on xj) Compute the residual vector

Description

Residual (residual when xi is regressed on xj) Compute the residual vector

Usage

residual_vec(xi, xj)

Arguments

Value

Residual vector after regression

Resolve a variable identifier to a 1-based integer index

Description

Shared by estimate_var_total_effect() and estimate_total_effect().

Usage

resolve_var_index(idx, arg_name, col_names, n_features)

Roll matrix rows (numpy np.roll equivalent, axis = 0)

Description

Shifts the rows of M downward by shift, wrapping the last shift rows around to the top. Used to build the lagged design for total-effect regression. The wrap-around contaminates the first shift rows, matching the Python reference (the effect is negligible for long series).

Usage

roll_rows(M, shift)

Arguments

Value

matrix with rolled rows

Population standard deviation (divided by n)

Description

Population standard deviation (divided by n)

Usage

sd_pop(x)

Search for candidate variables

Description

Search for candidate variables

Usage

search_candidate(U, Aknw, apply_prior_knowledge_softly, partial_orders)

Arguments

Value

list(Uc, Vj)

Causal order search via the kernel method

Description

Causal order search via the kernel method

Usage

search_causal_order_kernel(X, U, Uc, Vj)

Arguments

Value

Index of the selected variable

Causal order search via pwling

Description

Causal order search via pwling

Usage

search_causal_order_pwling(X, U, Uc, Vj)

Arguments

Value

Index of the selected variable

Causal order search via pwling, jointly across multiple groups

Description

Sample-size-weighted sum of the single-group pwling objective (search_causal_order_pwling()) across all groups in X_list, following the joint estimation objective of Shimizu (2012). Reuses the same standardize-once / correlation-matrix / antisymmetry optimizations as search_causal_order_pwling(), applied independently within each group before the weighted sum.

Usage

search_causal_order_pwling_multi(X_list, U, Uc, Vj)

Arguments

Value

Index of the selected variable

Select the VAR lag order by information criterion

Description

All candidate lag orders are compared on a common sample: the first max_lag observations are dropped for every candidate so that each VAR(lag) is estimated over the same time window (t = max_lag + 1 .. n). This mirrors statsmodels' VAR.select_order and makes the criteria comparable across lags (otherwise a longer lag would be scored on fewer observations).

Usage

select_var_lag(X, max_lag, criterion = "bic")

Arguments

Value

the selected lag order (integer)

Make this package's functions available on parallel cluster workers

Description

Used by all ⁠*_bootstrap()⁠ functions before dispatching iterations via parallel::parLapply(). Tries to have each worker library() the installed package; if that is not possible (e.g. during devtools::load_all() development, where the package is not installed), falls back to exporting every object in the namespace environment of a representative function from the algorithm being bootstrapped.

Usage

setup_cluster_worker(cl, fun)

Arguments

Value

NULL, invisibly. Called for the side effect of preparing cl.

Deterministic subsample for the Shapiro-Wilk test

Description

stats::shapiro.test() has a hard cap of 5000 observations, so larger inputs are thinned to SHAPIRO_MAX_N evenly spaced values. The thinning is deterministic on purpose: a random subsample would make the reported p-values change between calls and silently consume the caller's RNG stream (breaking downstream reproducibility). Even spacing over the input order does not distort the marginal distribution being tested.

Usage

shapiro_subsample(x)

Arguments

Value

x itself when length(x) <= SHAPIRO_MAX_N, otherwise a deterministic subsample of length SHAPIRO_MAX_N

Skewness (divided by n)

Description

Skewness (divided by n)

Usage

skewness_pop(x)

Arguments

Summarize the goodness-of-fit of a Direct LiNGAM model at once

Description

For a fitted Direct LiNGAM model, this verifies how well the two main assumptions on which LiNGAM relies (mutual independence of residuals and non-Gaussianity of residuals) hold, all at once, and displays the results together. Internally it calls get_error_independence_p_values() and test_residual_normality().

Usage

summary_lingam(
  X,
  lingam_result,
  independence_method = "spearman",
  normality_method = "shapiro",
  alpha = 0.05
)

Arguments

Details

Gaussian-likelihood-based criteria such as BIC/AIC are not included because they are theoretically inconsistent with LiNGAM's assumption that "the errors are non-Gaussian". Instead, the verification results of the assumptions themselves are summarized.

Value

A list of class lingam_summary, containing the following elements:

• n_variables, n_samples: Number of variables / number of observations

• causal_order: Causal order (variable-name labels)

• n_edges: Number of nonzero elements in the adjacency matrix (number of estimated edges)

• independence_p_values: Matrix of p-values from the independence test between residuals

• n_dependent_pairs, n_pairs: Number of pairs with p < alpha / total number of pairs

• min_independence_p: Minimum p-value of the independence test

• normality: Result of the normality test (a lingam_normality_test object)

• n_non_gaussian: Number of variables judged to be non-Gaussian

• alpha, independence_method, normality_method: The settings used

Examples

LiNGAM_sample_1000 <- generate_lingam_sample_6()

model <- lingam_direct(LiNGAM_sample_1000$data, reg_method = "ols")

summary_lingam(LiNGAM_sample_1000$data, model)

Test normality of residuals from Direct LiNGAM

Description

Calculate residuals (error terms) from the estimated adjacency matrix and test their normality. Since LiNGAM assumes non-Gaussian errors, rejecting normality (small p-value) supports the LiNGAM model assumption.

Usage

test_residual_normality(X, lingam_result, method = "shapiro", alpha = 0.05)

Arguments

Value

data.frame with test results for each variable

Examples

# Load the sample data
LiNGAM_sample_1000 <- generate_lingam_sample_6()

# Run Direct LiNGAM
result <- lingam_direct(LiNGAM_sample_1000$data, reg_method = "ols")

# Shapiro-Wilk (default)
test_residual_normality(LiNGAM_sample_1000$data, result)

Test the non-Gaussianity of VAR-LiNGAM residuals

Description

LiNGAM assumes the error terms are non-Gaussian, so rejecting normality (small p-value) supports the model assumption. By default the test is run on the LiNGAM innovations ⁠e_t = (I - B0) n_t⁠ (the independent errors the model assumes), where n_t are the stored VAR residuals; set on = "var" to test the reduced-form VAR residuals n_t directly instead.

Usage

test_varlingam_residual_normality(
  result,
  method = "shapiro",
  alpha = 0.05,
  on = c("innovations", "var")
)

Arguments

Value

a lingam_normality_test data frame (one row per variable), printed via print.lingam_normality_test().

References

Residual non-Gaussianity diagnostics inspired by the VARLiNGAM R code (Gauss_Tests) of Moneta, A., Entner, D., Hoyer, P. O., & Coad, A. (2013), Oxford Bulletin of Economics and Statistics, 75(5), 705-730. https://sites.google.com/site/dorisentner/publications/VARLiNGAM

Examples

s <- generate_varlingam_sample(n = 1000, seed = 42)
m <- lingam_var(s$data, lags = 1, reg_method = "ols", prune = FALSE)
test_varlingam_residual_normality(m)

Run several normality tests on VAR-LiNGAM residuals at once

Description

Convenience wrapper (analogous to the Moneta Gauss_Tests) that applies multiple normality tests to the residuals and returns a single table with one p-value column per method plus per-variable skewness and excess kurtosis. Methods whose optional package is unavailable are skipped with a warning.

Usage

test_varlingam_residual_normality_all(
  result,
  methods = c("shapiro", "ad", "lillie", "jb"),
  alpha = 0.05,
  on = c("innovations", "var")
)

Arguments

Value

a data frame with columns variable, skewness, kurtosis, one ⁠p_<method>⁠ column per method, and all_non_gauss (TRUE when every run test rejects normality for that variable).

References

Analogous to the multi-test residual check (Gauss_Tests) in the VARLiNGAM R code of Moneta, A., Entner, D., Hoyer, P. O., & Coad, A. (2013), Oxford Bulletin of Economics and Statistics, 75(5), 705-730. https://sites.google.com/site/dorisentner/publications/VARLiNGAM

Examples

s <- generate_varlingam_sample(n = 1000, seed = 42)
m <- lingam_var(s$data, lags = 1, reg_method = "ols", prune = FALSE)
test_varlingam_residual_normality_all(m, methods = c("shapiro", "jb"))

Convert a BootstrapResult to a tidy data.frame

Description

Returns a summary of the occurrence count, proportion, and effect size for each causal direction. Internally it calls get_causal_direction_counts(), so that function's arguments can be passed through ....

Usage

## S3 method for class 'BootstrapResult'
tidy(x, ...)

Arguments

Value

data.frame (from, to, count, proportion, ...)

Examples

dat <- generate_lingam_sample_6()
bs <- lingam_direct_bootstrap(dat$data, n_sampling = 30L, reg_method = "ols", seed = 42)
tidy(bs)

Convert an ImputationBootstrapResult to a tidy data.frame

Description

Collapses the imputation dimension with as_bootstrap_result() and then summarizes the causal direction counts like tidy.BootstrapResult().

Usage

## S3 method for class 'ImputationBootstrapResult'
tidy(x, aggregate = c("median", "mean"), ...)

Arguments

Value

data.frame (from, to, count, proportion, ...)

Examples

dat <- generate_lingam_sample_6(n = 200, seed = 1)$data
dat[sample(nrow(dat), 20), 1] <- NA
bs <- bootstrap_with_imputation(dat, n_sampling = 5L, n_repeats = 2L, seed = 42)
tidy(bs)

Convert a LiMResult to a tidy data.frame

Description

Converts the estimated adjacency matrix of a LiM model into a long-format data.frame with one edge per row, exactly like tidy.LingamResult().

Usage

## S3 method for class 'LiMResult'
tidy(x, threshold = 0, ...)

Arguments

Value

data.frame(from, to, estimate)

Examples

set.seed(1)
dat <- generate_lim_sample(n = 300)
model <- lingam_lim(dat$data, is_continuous = dat$is_continuous)
tidy(model)

Convert a LingamResult to a tidy data.frame

Description

Converts the estimated adjacency matrix into a long-format data.frame with one edge per row. Following the B[i, j] convention (the coefficient for j -> i), the from column is the cause and the to column is the effect. Convenient for visualization with ggplot2 or ggraph and for filtering with dplyr.

Usage

## S3 method for class 'LingamResult'
tidy(x, threshold = 0, ...)

Arguments

Value

data.frame(from, to, estimate). from/to are variable names (strings) and estimate is the causal coefficient. Returns a 0-row data.frame if there are no edges.

Examples

dat <- generate_lingam_sample_6()
model <- lingam_direct(dat$data, reg_method = "ols")
tidy(model)

Convert a MultiGroupBootstrapResult to a tidy data.frame

Description

Stacks each group's causal direction counts (via get_causal_direction_counts()) into a single data.frame with a group column in front. Arguments for get_causal_direction_counts() can be passed through ....

Usage

## S3 method for class 'MultiGroupBootstrapResult'
tidy(x, ...)

Arguments

Value

data.frame (group, from, to, count, proportion, ...)

Examples

mg <- generate_multi_group_sample()
bs <- lingam_multi_group_bootstrap(mg$data_list,
  n_sampling = 10L, reg_method = "ols", seed = 42
)
tidy(bs)

Convert a MultiGroupLingamResult to a tidy data.frame

Description

Stacks the per-group edge lists into a single long-format data.frame with a group column in front, one edge per row per group (the causal order is shared across groups but the coefficients differ).

Usage

## S3 method for class 'MultiGroupLingamResult'
tidy(x, threshold = 0, ...)

Arguments

Value

data.frame(group, from, to, estimate)

Examples

mg <- generate_multi_group_sample()
model <- lingam_multi_group(mg$data_list, reg_method = "ols")
tidy(model)

Convert a ParceLingamResult to a tidy data.frame

Description

Converts the estimated adjacency matrix into a long-format data.frame with one edge per row, like tidy.LingamResult(). NA entries of the adjacency matrix (variable pairs whose order could not be resolved / suspected latent confounding) are kept as rows with estimate = NA so they remain visible; drop them with e.g. subset(tidy(x), !is.na(estimate)) if not needed.

Usage

## S3 method for class 'ParceLingamResult'
tidy(x, threshold = 0, ...)

Arguments

Value

data.frame(from, to, estimate)

Examples

dat <- generate_parce_sample(n = 500, seed = 42)
model <- lingam_parce(dat$data)
tidy(model)

Convert an RCDResult to a tidy data.frame

Description

Converts the estimated adjacency matrix into a long-format data.frame with one edge per row, like tidy.LingamResult(). NA entries of the adjacency matrix (variable pairs suspected to share a latent confounder; marked in both directions) are kept as rows with estimate = NA so they remain visible; drop them with e.g. subset(tidy(x), !is.na(estimate)) if not needed.

Usage

## S3 method for class 'RCDResult'
tidy(x, threshold = 0, ...)

Arguments

Value

data.frame(from, to, estimate)

Examples

confounded <- generate_rcd_sample(n = 300, seed = 1)
model <- lingam_rcd(confounded$data)
tidy(model)

Validate the return value of a cd_fit function

Description

Validate the return value of a cd_fit function

Usage

validate_cd_fit_output(cd_res, p, n_datasets)

Arguments

Validate the return value of an imputer function

Description

Validate the return value of an imputer function

Usage

validate_imputer_output(datasets, X_boot)

Arguments

Validate the return value of lingam_direct()

Description

Validate the return value of lingam_direct()

Usage

validate_lingam_result(x)

Reject constant or perfectly collinear columns

Description

Constant and linearly dependent columns produce divisions by zero inside the pairwise regressions / standardization of the causal order search, which would otherwise surface as cryptic errors (e.g. "argument is of length zero") deep in the algorithm. The rank check runs on the centered matrix so that a column equal to another column plus a constant offset is also caught.

Usage

validate_no_degenerate_columns(X)

Arguments

Value

NULL, invisibly. Stops with an informative error on violation.

Validate the return value of lingam_parce()

Description

Validate the return value of lingam_parce()

Usage

validate_parce_result(x)

Validate a prior knowledge matrix and convert unknowns to NA

Description

Checks the shape and that every entry is -1 (unknown), 0 (no path), 1 (path), or NA (treated as unknown). Anything else (e.g. 0.5 or 2) would otherwise be silently interpreted as "path exists" by the candidate search, so it is rejected here.

Usage

validate_prior_knowledge(prior_knowledge, n_features)

Arguments

Value

the validated matrix with negative entries replaced by NA

Validate the return value of lingam_rcd()

Description

Validate the return value of lingam_rcd()

Usage

validate_rcd_result(x)

Validate and coerce n / seed for the sample generation functions

Description

Validate and coerce n / seed for the sample generation functions

Usage

validate_sample_args(n, seed)

Arguments

Value

list(n, seed) (both integers)

Total causal effect from a joined adjacency matrix (graph-based)

Description

Computes the total effect by summing path products over the time-expanded graph, reusing calculate_total_effect(). Port of the Python reference estimate_total_effect2; used internally by the VAR-LiNGAM bootstrap.

Usage

var_total_effect_graph(am_joined, from_index, to_index)

Arguments

Value

the total effect (scalar)

Build a zero-adjacency LingamResult

Description

A stand-in LingamResult whose adjacency matrix is all zeros, so that lingam_residuals() returns its input unchanged. This lets the VAR diagnostics reuse the Direct LiNGAM residual routines on an already-computed residual matrix.

Usage

zero_lingam_result(p)

Arguments

Value

a LingamResult with a p x p zero adjacency matrix