## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
set.seed(42)

## ----install, eval=FALSE------------------------------------------------------
# # Step 1: install glmgraph (not on CRAN; required for ivgl and ivgl_s)
# devtools::install_github("cran/glmgraph")
# 
# # Step 2: install ivgls
# install.packages("ivgls")

## ----graphs-------------------------------------------------------------------
library(ivgls)

A <- make_graph(p = 30, type = "chain")
L <- get_laplacian(A)

cat("Number of edges:", sum(A) / 2, "\n")

## ----beta---------------------------------------------------------------------
set.seed(1)

bobj <- generate_beta(A, s2 = 4, signal = 2, pattern = "smooth")

cat("Active nodes:", sort(bobj$active_set), "\n")
cat("True effects at active nodes:\n")
print(round(bobj$beta_true[bobj$active_set], 3))

## ----data---------------------------------------------------------------------
dat <- generate_data(
  n              = 100,
  p              = 30,
  q              = 100,
  s_alpha        = 5,
  alpha_strength = 3,
  beta_true      = bobj$beta_true
)

cat("Y:", length(dat$Y),
    "| X:", nrow(dat$X), "x", ncol(dat$X),
    "| Z:", nrow(dat$Z), "x", ncol(dat$Z), "\n")
cat("Invalid instrument indices:", which(dat$alpha_true != 0), "\n")

## ----ivlasso------------------------------------------------------------------
beta_ivl <- iv_lasso(dat$Y, dat$X, dat$Z)

cat("IV-LASSO selected nodes:", which(abs(beta_ivl) > 1e-4), "\n")
cat("IV-LASSO MCC:",
    round(get_mcc(bobj$active_set, which(abs(beta_ivl) > 1e-4), p = 30), 3), "\n")

## ----glmgraph-estimators------------------------------------------------------
if (requireNamespace("glmgraph", quietly = TRUE)) {

  beta_ivgl <- ivgl(dat$Y, dat$X, dat$Z, L)
  fit_s     <- ivgl_s(dat$Y, dat$X, dat$Z, L, max_iter = 10)

  cat("IVGL selected nodes:", which(abs(beta_ivgl) > 1e-4), "\n")
  cat("IVGL-S selected nodes:", which(abs(fit_s$beta) > 1e-4), "\n")

  cat("IVGL MCC:",
      round(get_mcc(bobj$active_set, which(abs(beta_ivgl) > 1e-4), p = 30), 3), "\n")
  cat("IVGL-S MCC:",
      round(get_mcc(bobj$active_set, which(abs(fit_s$beta) > 1e-4), p = 30), 3), "\n")

  cat("Detected invalid instruments:", which(abs(fit_s$alpha) > 1e-4), "\n")

} else {
  message(
    "Package 'glmgraph' is not installed. ",
    "ivgl() and ivgl_s() are unavailable.\n",
    "Install with: devtools::install_github(\"cran/glmgraph\")"
  )
}

## ----eval---------------------------------------------------------------------
supp_ivl <- which(abs(beta_ivl) > 1e-4)
metrics  <- eval_support(bobj$active_set, supp_ivl, p = 30)
print(round(metrics, 3))

## ----corrupt------------------------------------------------------------------
set.seed(2)
A_corrupted <- corrupt_graph(A, corruption_rate = 0.20)

cat("Original edges:", sum(A) / 2, "\n")
cat("Corrupted edges:", sum(A_corrupted) / 2, "\n")

## ----simulation---------------------------------------------------------------
if (requireNamespace("glmgraph", quietly = TRUE)) {

  res <- run_simulation(
    B              = 5,
    n              = 100,
    p              = 30,
    q              = 100,
    graph_type     = "chain",
    signal_pattern = "smooth",
    s2             = 4,
    signal         = 2,
    s_alpha        = 5,
    alpha_strength = 3
  )

  aggregate(cbind(MCC, MSE, TPR, FPR) ~ Method, data = res, FUN = mean)

} else {
  message(
    "Package 'glmgraph' is not installed. ",
    "run_simulation() is unavailable.\n",
    "Install with: devtools::install_github(\"cran/glmgraph\")"
  )
}