Applied examples: causal models and time series

library(CCI)

This vignette shows two common uses of conditional independence tests: checking a causal model against data, and testing whether the past of one time series predicts another. For the basics, see vignette("Testing-CI-with-CCI", package = "CCI").

Testing a causal model (DAG)

A causal model drawn as a directed acyclic graph (DAG) implies a set of conditional independencies: every pair of variables that is not connected by an arrow is independent given some set of other variables. If the data contradict one of them, the DAG is wrong, for instance because an arrow is missing. Testing all implied conditional independencies is therefore a way to check a DAG against data (Textor et al. 2016). The dagitty package finds the implied conditional independencies, and CCI can test them, also when the relationships are non-linear.

We use a DAG with five variables:

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

and simulate data from it, with non-linear relationships:

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)

impliedConditionalIndependencies() lists the independencies implied by the DAG. Each one is turned into a CCI formula. An empty conditioning set gives an unconditional test, written with | 1:

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
#> A _||_ D | B, C
#> A _||_ E | D
#> A _||_ E | B, C
#> B _||_ C | A
#> B _||_ E | D
#> C _||_ E | D
formulas <- lapply(implied, ci_to_formula)

Then each is tested. Since several hypotheses are tested, the p-values are adjusted for multiple testing with Holm’s method. With many tests, the adjusted p-values need to be small, so we use parametric p-values (parametric = TRUE), which are not limited by the number of Monte Carlo samples:

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
#>      hypothesis p_value p_adjusted
#> 1 A ~ D | B + C   1.000          1
#> 2     A ~ E | D   0.961          1
#> 3 A ~ E | B + C   1.000          1
#> 4     B ~ C | A   0.171          1
#> 5     B ~ E | D   0.631          1
#> 6     C ~ E | D   0.666          1

Rejected at the 5% level after adjustment: none. The data are consistent with the DAG. Now suppose we had drawn the DAG without the arrow from \(C\) to \(D\):

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
#>    hypothesis  p_value p_adjusted
#> 1   A ~ D | B 1.33e-02   1.19e-01
#> 2   A ~ E | D 9.61e-01   1.00e+00
#> 3   A ~ E | B 8.63e-01   1.00e+00
#> 4   B ~ C | A 1.71e-01   1.00e+00
#> 5   B ~ E | D 6.31e-01   1.00e+00
#> 6   C ~ D | B 3.09e-07   3.09e-06
#> 7   C ~ D | A 5.41e-01   1.00e+00
#> 8   C ~ E | D 6.66e-01   1.00e+00
#> 9   C ~ E | B 3.57e-02   2.86e-01
#> 10  C ~ E | A 7.39e-01   1.00e+00

Rejected at the 5% level after adjustment: C ~ D | B. An independence between \(C\) and \(D\) only holds if every path between them goes through the conditioning variables, so a rejection shows that the DAG is missing a connection between \(C\) and \(D\). Not every independence that involves the missing arrow needs to be rejected: some dependencies are weak and hard to detect. A rejected independence points to where the DAG should be revised, but it does not say which arrow to add; several changes to the DAG could explain the same rejection.

Note that non-rejection does not prove that the DAG is right: other DAGs can imply the same independencies, and a test may lack the power to detect a weak dependence.

Time series

CCI can also be used with time series, by putting lagged variables in the formula. A common question is whether the past of one series, \(X\), helps to predict another series, \(Y\), beyond the past of \(Y\) itself (Granger causality). We simulate a series \(X\) that follows an autoregressive process, and a series \(Y\) that depends non-linearly on the past of \(X\) and has a trend:

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)

Does the past of \(X\) predict \(Y\) beyond the past of \(Y\) and the trend? It should, since \(Y\) depends on \(X_{t-1}\):

summary(CCI.test(Y ~ X_lag1 | Y_lag1 + Y_lag2 + Time, data = ts_data, nperm = 100,
                 seed = 1, progress = FALSE))
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using rf 
#> Formula:   Y ~ X_lag1 | Y_lag1 + Y_lag2 + Time 
#> Permutations:  100 
#> Metric:    RMSE 
#> Tail:      left 
#> Statistic: 1.536 
#> P-value:   0.009901 
#> 
#> MC sample:   0.93

Does the past of \(Y\) predict \(X\) beyond the past of \(X\)? It should not, since \(X\) follows its own autoregressive process and does not depend on \(Y\):

summary(CCI.test(X ~ Y_lag1 | X_lag1 + X_lag2, data = ts_data, nperm = 100,
                 seed = 1, progress = FALSE))
#> 
#> Computational Conditional Independence Test
#> --------------------------------------------
#> Method:    CCI test using rf 
#> Formula:   X ~ Y_lag1 | X_lag1 + X_lag2 
#> Permutations:  100 
#> Metric:    RMSE 
#> Tail:      left 
#> Statistic: 1.082 
#> P-value:   0.5941 
#> 
#> MC sample:   0.93

Include the lags that are needed to describe each series (here two lags of \(X\)), and variables like Time for trends. Otherwise, the past of one series can predict the other only because it carries information about the missing lags or the trend.

A caution. CCI splits the data randomly into training and test parts, which assumes that the observations are exchangeable. Time series observations are dependent over time, so the test is approximate. It works best with long series, when the lags in the conditioning set capture most of the dependence over time.

References

Textor, Johannes, Benito van der Zander, Mark S. Gilthorpe, Maciej Liśkiewicz, and George T. H. Ellison. 2016. “Robust Causal Inference Using Directed Acyclic Graphs: The R Package ’Dagitty’.” International Journal of Epidemiology 45 (6): 1887–94. https://doi.org/10.1093/ije/dyw341.