## ----include = FALSE----------------------------------------------------------
has_dagitty <- requireNamespace("dagitty", quietly = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4,
  message = FALSE
)

## ----setup--------------------------------------------------------------------
library(CCI)

## ----eval = has_dagitty-------------------------------------------------------
library(dagitty)
true_dag <- dagitty("dag {
  A -> B
  A -> C
  B -> D
  C -> D
  D -> E
}")
plot(graphLayout(true_dag))

## -----------------------------------------------------------------------------
simulate_dag_data <- function(n) {
  A <- rnorm(n)
  B <- sin(2 * A) + rnorm(n, sd = 0.5)
  C <- A^2 + rnorm(n, sd = 0.5)
  D <- B * C + rnorm(n, sd = 0.5)
  E <- tanh(D) + rnorm(n, sd = 0.3)
  data.frame(A, B, C, D, E)
}
set.seed(1)
dag_data <- simulate_dag_data(600)

## ----eval = has_dagitty-------------------------------------------------------
ci_to_formula <- function(ci) {
  z <- if (length(ci$Z) == 0) "1" else paste(ci$Z, collapse = " + ")
  as.formula(paste(ci$X, "~", ci$Y, "|", z))
}
implied <- impliedConditionalIndependencies(true_dag)
implied
formulas <- lapply(implied, ci_to_formula)

## ----eval = has_dagitty-------------------------------------------------------
test_dag <- function(formulas, data) {
  p <- vapply(formulas, function(f) {
    CCI.test(f, data = data, nperm = 60, parametric = TRUE, seed = 1, progress = FALSE)$p.value
  }, numeric(1))
  data.frame(hypothesis = vapply(formulas, function(f) paste(deparse(f), collapse = ""), ""),
             p_value = signif(p, 3),
             p_adjusted = signif(p.adjust(p, method = "holm"), 3))
}
true_results <- test_dag(formulas, dag_data)
true_results

## ----include = FALSE----------------------------------------------------------
rejected <- function(results) {
  if (!has_dagitty) return("")
  hyp <- results$hypothesis[results$p_adjusted <= 0.05]
  if (length(hyp) == 0) "none" else paste0("`", hyp, "`", collapse = ", ")
}

## ----eval = has_dagitty-------------------------------------------------------
wrong_dag <- dagitty("dag {
  A -> B
  A -> C
  B -> D
  D -> E
}")
wrong_results <- test_dag(lapply(impliedConditionalIndependencies(wrong_dag), ci_to_formula),
                          dag_data)
wrong_results

## -----------------------------------------------------------------------------
simulate_series <- function(n) {
  X <- as.numeric(arima.sim(n = n, list(ar = c(0.9, -0.5))))
  Y <- numeric(n)
  for (t in 3:n) {
    Y[t] <- 0.01 * t + 1.2 * X[t - 1] + 0.7 * X[t - 2] + 0.5 * X[t - 1] * X[t - 2] + rnorm(1)
  }
  data.frame(Time = seq_len(n), X = X, Y = Y)
}
lag <- function(x, k) c(rep(NA, k), x[seq_len(length(x) - k)])

set.seed(1993)
ts_data <- simulate_series(1000)
ts_data$X_lag1 <- lag(ts_data$X, 1)
ts_data$X_lag2 <- lag(ts_data$X, 2)
ts_data$Y_lag1 <- lag(ts_data$Y, 1)
ts_data$Y_lag2 <- lag(ts_data$Y, 2)
ts_data <- na.omit(ts_data)

## -----------------------------------------------------------------------------
summary(CCI.test(Y ~ X_lag1 | Y_lag1 + Y_lag2 + Time, data = ts_data, nperm = 100,
                 seed = 1, progress = FALSE))

## -----------------------------------------------------------------------------
summary(CCI.test(X ~ Y_lag1 | X_lag1 + X_lag2, data = ts_data, nperm = 100,
                 seed = 1, progress = FALSE))