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

## ----setup--------------------------------------------------------------------
library(lingamr)

## ----sample1_matrix6_1--------------------------------------------------------
x1k <- generate_lingam_sample_6(n = 1000)

x1k$data |>
  head()

## ----sample1_matrix6_2--------------------------------------------------------
x1k$true_adjacency

## ----sample1_matrix6_3--------------------------------------------------------
x1k$true_adjacency |>
  plot_adjacency(
    labels  = colnames(x1k$data),
    title   = "True causal structure",
    rankdir = "TB",
    shape   = "circle"
  )

## ----lingam_direct------------------------------------------------------------
model <- x1k$data |>
  lingam_direct()

## ----causal_order-------------------------------------------------------------
# index number
model$causal_order

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

## ----adjacency_matrix---------------------------------------------------------
model$adjacency_matrix |>
  round(3)

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

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

## ----autoplot, eval = requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("igraph", quietly = TRUE)----
ggplot2::autoplot(model)

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

round(total_effects, 3)

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

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

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

pk1

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

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

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

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

## ----LiNGAM_with_pk_2---------------------------------------------------------
model_pk1$adjacency_matrix |>
  round(3)

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

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

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

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

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

## ----Independence_between_error_variables1------------------------------------
result <- x1k$data |>
  lingam_direct()

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

round(p_vals, 3)

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

truth$true_adjacency |>
  round(1)

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

## ----nongauss_uniform---------------------------------------------------------
fit_uniform <- lingam_direct(truth$data)

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

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

## ----nongauss_uniform_plot----------------------------------------------------
fit_uniform$adjacency_matrix |>
  plot_adjacency(
    labels = colnames(truth$data),
    true_B = truth$true_adjacency,
    title  = "Estimated (uniform errors)"
  )

## ----nongauss_gaussian--------------------------------------------------------
gauss <- generate_lingam_sample_6(noise_dist = "gaussian")
fit_gauss <- lingam_direct(gauss$data)

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

fit_gauss$adjacency_matrix |>
  round(1)

## ----nongauss_gaussian_plot---------------------------------------------------
fit_gauss$adjacency_matrix |>
  plot_adjacency(
    labels = colnames(gauss$data),
    true_B = truth$true_adjacency,
    title  = "Estimated (Gaussian errors)"
  )

## ----nortest------------------------------------------------------------------
# Shapiro-Wilk (default)
x1k$data |>
  test_residual_normality(result)

## ----qqplot-------------------------------------------------------------------
x1k$data |>
  plot_residual_qq(result)

## ----summary_lingam-----------------------------------------------------------
x1k$data |>
  summary_lingam(result)

## ----lingam_direct_bootstrap--------------------------------------------------
bs_model <- x1k$data |>
  lingam_direct_bootstrap(n_sampling = 100L, seed = 42)

bs_model

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

## ----get_causal_direction_counts----------------------------------------------
bs_model |>
  get_causal_direction_counts(labels = names(x1k$data))

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

bs_adjacency_matrix |>
  round(3)

## ----plot_bootstrap-----------------------------------------------------------
bs_adjacency_matrix |>
  plot_adjacency(
    labels    = colnames(x1k$data),
    title     = "Estimated (with Bootstrap)",
    rankdir   = "TB",
    shape     = "circle",
    fillcolor = "lightgreen"
  )

## ----get_probabilities--------------------------------------------------------
bs_model |>
  get_probabilities()

## ----get_total_causal_effects-------------------------------------------------
bs_model |>
  get_total_causal_effects()

## ----plot_bootstrap_probabilities---------------------------------------------
bs_model |>
  plot_bootstrap_probabilities()

## ----get_causal_order_stability-----------------------------------------------
bs_model |>
  get_causal_order_stability(labels = names(x1k$data))

## ----broom--------------------------------------------------------------------
# Convert the estimated adjacency matrix to an edge list
tidy(model)

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

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

## ----sample1_matrix10---------------------------------------------------------
x10k <- generate_lingam_sample_10(n = 10000)

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

## ----ica_run, eval = requireNamespace("pcalg", quietly = TRUE)----------------
d_cmp <- generate_lingam_sample_6(n = 1000, seed = 42)

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

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

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

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

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

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

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

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

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

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

## ----large_timing-------------------------------------------------------------
d10 <- generate_lingam_large_sample(p = 10, n = 500, seed = 42)
d15 <- generate_lingam_large_sample(p = 15, n = 500, seed = 42)

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

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

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

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

## ----large_causal_order-------------------------------------------------------
# Estimated causal order
r10$causal_order

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

## ----large_tidy---------------------------------------------------------------
tidy(r10) |>
  head(10)

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

hd_result$causal_order
round(hd_result$adjacency_matrix, 3)

## ----high_dim_wide------------------------------------------------------------
wide_sample <- generate_lingam_large_sample(p = 30, n = 25, seed = 1)
wide_result <- lingam_high_dim(wide_sample$data)

wide_result$causal_order

## ----paradox_generate---------------------------------------------------------
paradox <- generate_lingam_paradox_data(n = 2000L, seed = 42)

head(paradox$data)

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

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

## ----paradox_fit--------------------------------------------------------------
model_p <- lingam_direct(paradox$data)

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

## ----paradox_adjacency--------------------------------------------------------
model_p$adjacency_matrix |>
  round(3)

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

## ----paradox_bootstrap--------------------------------------------------------
bs_paradox <- paradox$data |>
  lingam_direct_bootstrap(n_sampling = 100L, seed = 42)

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

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

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

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

## ----lingam_var---------------------------------------------------------------
model <- lingam_var(s$data, lags = 1)
model

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

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

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

## ----varlingam_stationarity---------------------------------------------------
check_var_stationarity(model)

## ----varlingam_normality------------------------------------------------------
test_varlingam_residual_normality(model)

## ----varlingam_normality_all--------------------------------------------------
test_varlingam_residual_normality_all(model, methods = c("shapiro", "jb"))

## ----varlingam_qq, fig.width = 7, fig.height = 3------------------------------
plot_varlingam_residual_qq(model)

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

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

## ----varlingam_bootstrap------------------------------------------------------
bs_var <- lingam_var_bootstrap(
  s$data,
  n_sampling = 100L,
  seed       = 42,
  verbose    = FALSE
)

## ----varlingam_probabilities--------------------------------------------------
round(get_var_probabilities(bs_var), 2)

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

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

## ----lim_sample---------------------------------------------------------------
set.seed(1)
lim_dat <- generate_lim_sample(n = 2000)
head(lim_dat$data)
lim_dat$is_continuous

## ----lim_fit------------------------------------------------------------------
lim_result <- lingam_lim(lim_dat$data, is_continuous = lim_dat$is_continuous)
print(lim_result)

## ----lim_order----------------------------------------------------------------
colnames(lim_dat$data)[lim_result$causal_order]

## ----multi_group_sample-------------------------------------------------------
mg <- generate_multi_group_sample(n = c(1000, 1000), seed = 42)
lapply(mg$data_list, head, 3)

## ----multi_group_fit----------------------------------------------------------
mg_result <- lingam_multi_group(mg$data_list, reg_method = "ols")
print(mg_result)

## ----multi_group_extract------------------------------------------------------
g1 <- get_group_result(mg_result, "group1")
class(g1)

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

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

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

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

## ----bootstrap_with_imputation_convert, eval = requireNamespace("mice", quietly = TRUE)----
bs_na <- as_bootstrap_result(bwi, aggregate = "median")
get_probabilities(bs_na)

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

## ----parce_fit----------------------------------------------------------------
parce_result <- lingam_parce(confounded$data, reg_method = "ols")
print(parce_result)

## ----parce_block--------------------------------------------------------------
parce_result$causal_order[[1]]
parce_result$adjacency_matrix[confounded$confounded_pair, confounded$confounded_pair]

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

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

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

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

## ----rcd_fit------------------------------------------------------------------
rcd_result <- lingam_rcd(rcd_confounded$data)
print(rcd_result)

## ----rcd_ancestors------------------------------------------------------------
rcd_result$ancestors_list
rcd_result$adjacency_matrix[rcd_confounded$confounded_pair, rcd_confounded$confounded_pair]

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

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

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

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

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