--- title: "Learning cyclic causal graphs with backShift" author: "Christina Heinze-Deml" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: fig_caption: yes vignette: > %\VignetteIndexEntry{backShift demo} %\VignetteEngine{knitr::rmarkdown} \usepackage[utf8]{inputenc} status: process --- This is an example for how to run [backShift](http://arxiv.org/abs/1506.02494), an algorithm to estimate the connectivity matrix of a directed (possibly cyclic) graph with hidden variables. The underlying system is required to be linear and we assume that observations under different shift interventions are available. In this example, we simulate a data set with the function `simulateInterventions()` provided in this package. ## Options ### Preliminaries First, we load the `backShift` package as well as the packages `ggplot2` and `fields` which are needed to visualize the results. Additionally, we set the random seed for reproducability. ```{r, message = FALSE} # load package require(backShift) require(ggplot2) require(fields) # set seed seed <- 1 set.seed(seed) ``` ### Algorithm backShift exploits differences between covariance or Gram matrices. To use covariance matrices set `useCov` to `TRUE`. ```{r} # use covariance matrix instead of Gram matrix useCov <- TRUE ``` ### Stability selection backShift estimates the connectivity matrix of a directed (possibly cyclic) graph. In order to control the number of falsly selected edges we can use [stability selection](http://arxiv.org/abs/0809.2932) which provides a finite sample control for the expected number of false discoveries. If stability selection should not be used, set the parameter `EV` to 0. ```{r} # bound on expected number of false selections for stability selection EV <- 2 # selection threshold for stability selection thres <- 0.75 ``` To visualize the point estimate of the connectivity matrix, the coefficients are thresholded at the absolute value given by `thres.pe`. So edges with coefficients smaller than `thres.pe` in absolute value are not displayed. ```{r} # threshold for point estimate thres.pe <- 0.25 ``` ### Simulation A data set can be simulated with `simulateInterventions()`. This function takes the true connectivity matrix as an input. A connectivity matrix can be generated with `generateA()` or you can use the connectivity matrix provided through `data("exampleAdjacencyMatrix")`. #### Connectivity matrix In the following code, we set whether we "provide" the adjacency matrix __A__ -- e.g. by loading the example with `data("exampleAdjacencyMatrix")` or by creating an adjacency matrix ourselves -- or whether we generate the connectivity matrix __A__ with `generateA()`. NOTE: The entry __A__{_ij_} contains the edge from node _i_ to node _j_. ```{r, results = "hide"} # number of variables p <- 10 # set whether to provide or to generate the adjacency matrix A providedA <- TRUE if(providedA){ data("exampleAdjacencyMatrix") A <- exampleAdjacencyMatrix p <- 10 }else{ # parameters to generate A # should A be cyclic? cyclic <- TRUE # expected number of neighbors per node expNumNeigh <- 0.1*p # range for coefficients minCoef <- 0.3 maxCoef <- 0.8 ## Generate A ------- cat("Generating A...\n") A.gen.result <- generateA(p, expNumNeigh, minCoef, maxCoef, cyclic) A <- A.gen.result$A cat("A has a cycle of size", A.gen.result$sizeCycle, "\n") } ``` #### Observations under different shift interventions The following code generates a data set. In addition to the size of the data set, we can specify the following options: * whether to also simulate observational data (`simulateObs`), * whether hidden variables should be present (`hidden`) * whether the location of the interventions should be known (`knownInterventions`; note that this is not needed for backShift) * if `knownInterventions` is `TRUE`, `fracVarInt` gives the fraction of variables that are intervened on in each environment * `intMult` determines the magnitude of the intervention variances (please see the [manuscript](http://arxiv.org/abs/1506.02494) for more details) * `noiseMult` determines the noise variance * `nonGauss` specifies whether to generate non-Gaussian noise ```{r, results = "hide"} # number of observations n <- 10000 # number of environments G <- 10 # also simulate observational data? simulateObs <- TRUE # should hidden vars be included? hidden <- FALSE # should the location of the interventions be known? knownInterventions <- FALSE # if the location of the interventions is known, how many vars. should # be intervened on in each environment (as a fraction of p) fracVarInt <- 0.5 # multiplier for interventions (m_I in manuscript) intMult <- 1.5 # multiplier for interventions (m_e in manuscript) noiseMult <- 1 # simulate non-Gaussian noise? nonGauss <- FALSE ## Simulate data ------- cat("Simulating the data...\n") simulation.res <- simulateInterventions(n, p, A, G, intMult, noiseMult, nonGauss, hidden, knownInterventions, fracVarInt, simulateObs, seed) # extract X, environment vector and index of observational data X <- simulation.res$X env <- simulation.res$environment baseInd <- simulation.res$configs$indexObservationalData ``` ## Running backShift We can now run backShift. Since we also generated observational data, we provide the corresponding index as `baseInd`. This is useful for estimating the intervention variances (see below). ```{r} ## Run backShift ------- backshift.res <- backShift(X, env, covariance=useCov, ev=EV, threshold=thres, baseSettingEnv = baseInd, tolerance = 1e-6, verbose = FALSE) ``` ## Results ### Estimated graphs Plot true graph: ```{r, fig.show='hold', fig.width=6, fig.height=6, fig.align='center', results = "hide"} # extract estimates Ahat <- backshift.res$Ahat Ahat.structure <- backshift.res$AhatAdjacency # compute performance metrics and plot result cat("Plotting true graph... \n") plotGraphEdgeAttr(estimate = A, plotStabSelec = FALSE, labels = colnames(A), thres.point = 0, thres.stab = thres, main = "True graph") ``` Plot point estimate (thresholded at `thres.pe`): The edge intensity reflects the relative magnitude of the coefficients. ```{r, fig.show='hold', fig.width=6, fig.height=6, fig.align='center', results="hide"} cat("Plotting point estimate, thresholded at", thres.pe,"... \n") plotGraphEdgeAttr(estimate = Ahat, plotStabSelec = FALSE, labels = colnames(A), thres.point = thres.pe, thres.stab = thres, main = paste("Point estimate thresholded at", thres.pe)) ``` Plot stability selection result: The edges thickness indicates how often an edge was selected in the stability selection procedure. ```{r, fig.show='hold', fig.width=6, fig.height=6, fig.align='center', results = "hide"} cat("Plotting stability selection result... \n") plotGraphEdgeAttr(estimate = Ahat.structure, plotStabSelec = TRUE, labels = colnames(A), thres.point = thres.pe, edgeWeights = Ahat, thres.stab = thres, main = "Stability selection result") ``` ### Metrics `metricsThreshold` computes the structural hamming distance (SHD), true positive rate (TPR)/recall, false positive rate (FPR) and precision. The connectivity matrix gets thresholded at `thres` prior to computing these metrics. ```{r, results = "asis"} # metrics for point estimate, thresholded at thres.pe metricsThresholdedA <- metricsThreshold(A, Ahat, thres = thres.pe) # metrics for stability selection result metricsStabSelection <- metricsThreshold(A, Ahat.structure, thres = 0) ``` ```{r, results = "asis", echo = FALSE} colnames(metricsThresholdedA)[3] <- colnames(metricsStabSelection)[3] <- "TPR/Recall" knitr::kable(metricsThresholdedA, digits = 4, caption = paste("Metrics for point estimate, thresholded at", thres.pe)) knitr::kable(metricsStabSelection, digits = 4, caption = "Metrics for stability selection result") ``` ## Estimating the intervention variances The location and the strength of the shift interventions can be estimated from the data (up to an offset if no observational data is present). These estimates are returned by the function `backShift()` as a _G_ x _p_- dimensional matrix `varianceEnv` where _G_ is the number of environments and _p_ is the number of variables. The _ij_-th entry contains the difference between the estimated intervention variance of variable _j_ in environment _i_ and the estimated intervention variance of variable _j_ in the base setting (given by input parameter `baseSettingEnv`). ```{r, message=FALSE, fig.width=10, fig.height=7.5} plotInterventionVars(backshift.res$varianceEnv, simulation.res$interventionVar) ``` ## Checking the model assumptions We can check the model assumptions to some extent by the success or failure of the joint diagonalization procedure. The plots below show that the joint diagonalization succeeded for all environments as all matrices are diagonal. The values are scaled to [0,1] to allow for comparability accross plots. ```{r, fig.keep='all', fig.width=3, fig.height=3} for(i in 1:G){ plotDiagonalization(estConnectivity = backshift.res$Ahat, X = X, env = env, whichEnv = i) } ```