notebook_02_subgraph_selection_evaluation_publish.Rmd

---
title: "Subgraph selection and evaluation"
author: David Reiffenscheidt and Oleksandr Zadorozhnyi
date: "12/09/2023"
output: 
  html_document:
    md_document: default
---


***Setup of the problem***

Given a large graph $G = (V, E)$, where $V$ represents the set of nodes and $E$ represents the set of edges, the goal is to select a subgraph $S = (V_s, E_s)$ from G such that S is a meaningful and informative representation of the original graph $G$. The subgraph selection problem involves finding an optimal or near-optimal subgraph that satisfies certain criteria or objectives. 

In the context of graphical modelling the problem of subgraph selection corresponds to the problem of preserving the conditional independence relationship in the subgraph $S$ which has to be transferred from the underlying graph $G$.

Furthermore, it is important to be able to evaluate the results of the structure learning algorithms (learned on the data) taking as a ground truth the selected subgraph.  

```{r}
### loading packages
library("bnlearn")
library("qgraph")
library("igraph")
library("pcalg")
```
Reading the files from Zenodo entry in the community

```{r}
# necessary libraries to use for Zenodo REST-API
library(zen4R)

zenodo <- ZenodoManager$new(
  logger = "INFO" # use "DEBUG" to see detailed API operation logs, use NULL if you don't want logs at all
)

# downloading files using zenodo doi and reading from the file 
rec1 <- zenodo$getRecordByDOI("10.5281/zenodo.7676616")
files <- rec1$listFiles(pretty = TRUE)

#create a folder where to download files from record
dir.create("download_zenodo")

#download files
rec1$downloadFiles(path = "download_zenodo")
downloaded_files <- list.files("download_zenodo")

zipF = sprintf("download_zenodo/%s",downloaded_files)

# unzipping in the current folder
unzip(zipF,exdir = "./")

alarm_name = list.files(tools::file_path_sans_ext(downloaded_files))[1]

path_to_file = paste0(tools::file_path_sans_ext(downloaded_files),"/",alarm_name,"/",alarm_name,".csv")

df = read.csv(path_to_file)
head(df)
```


We need to transform data first

```{r}
for (item in colnames(df)){
  df[,item] = as.factor(df[,item])
}
```


```{r}
data("alarm")
alarm_df <- as.data.frame(na.omit(df))

p = length(names(alarm))
n = dim(alarm)[1]
for (i in c(1:p)) {
  alarm_df[,i]<-as.numeric(alarm_df[,i])
}
```

#### Applying nonparanormal transformation to standardize the data.
```{r}
library("huge")
alarm_df <- huge.npn(alarm_df)
#####
```

Defining "true" graph as proposed for the ALARM dataset in bnlearn
```{r}
# "True" Graph ALARM
dag_alarm = empty.graph(names(alarm))
modelstring(dag_alarm) = paste0("[HIST|LVF][CVP|LVV][PCWP|LVV][HYP][LVV|HYP:LVF]",                    "[LVF][STKV|HYP:LVF][ERLO][HRBP|ERLO:HR][HREK|ERCA:HR][ERCA][HRSA|ERCA:HR]",
"[ANES][APL][TPR|APL][ECO2|ACO2:VLNG][KINK][MINV|INT:VLNG][FIO2]",            "[PVS|FIO2:VALV][SAO2|PVS:SHNT][PAP|PMB][PMB][SHNT|INT:PMB][INT]",              "[PRSS|INT:KINK:VTUB][DISC][MVS][VMCH|MVS][VTUB|DISC:VMCH]",                    "[VLNG|INT:KINK:VTUB][VALV|INT:VLNG][ACO2|VALV][CCHL|ACO2:ANES:SAO2:TPR]",
"[HR|CCHL][CO|HR:STKV][BP|CO:TPR]", sep = "")

qgraph(dag_alarm, legend.cex = 0.3,
       asize=2,edge.color="black", vsize= 4)
```

Selection of the set of nodes for subsetting
```{r}
### Subgraph 
subgraph_nodes <- c("INT","VALV","LVF","LVV","PCWP","HR","CCHL","CVP","HYP","HRSA","ERCA")
```

**First** procedure (bnlearns subgraph function). It is implemented through a simple subsetting of the edges which are adjacent to the vertices contained in the subgraph_nodes set. 

```{r}
procedure1 <- bnlearn::subgraph(dag_alarm, subgraph_nodes)
qgraph(procedure1, legend.cex = 0.3,
       asize=2,edge.color="black", vsize= 5)
```

**Second** procedure selects the subgraph based on the following euristics.  Given the ground truth DAG $G=(V,E)$ and subset of vertices $V_{s} \subset V$ the goal is to find the corresponding set of vertices $E_{s}$ such that for $(V_s,E_s)$ the structure of the distribution for $ P_{V_s}$ does not contradict the structure of the distribution $P_{V}$ so that we the task of structure estimation (and then benchmarking) on $(V_s,E_s)$ can be done with a new ground truth. 

Let $G = (V,E)$ be the original directed acyclic graph, and $P_{V}$ is the joint distribution of random variables from $V$.

We take two and check whether these vertices are $d-$connected given all others in $V_{s}$. If they are not $d$-connected, there is no association (no arrow in any direction). Otherwise, we have correlation (with unknown direction). If one of the directions leads to the cycle in the original graph, we resolve it and keep the other direction, otherwise we keep both directions and then get the CPDAG.

**Third** procedure which сreats the PAG for the observable vertices in the subset-graph with respect to the latent variables in the vertices set $V \setminus V^{'}$.

PAG is the graph which is used to represent causal relationships in situations where the exact underlying causal structure is not fully known or observable. In our case we observe the variables in the set $V^{'}$. Based on the covariance matrix from the observed data in $V^{'}$ we generate. 


```{r}
######### Extract subgraph function

combn(subgraph_nodes, 2)[2,1]
dim(combn(subgraph_nodes,2))[2]

extract_subgraph <- function(dag, nodes){
  sg <- bnlearn::subgraph(dag,nodes) # procedure 1 (to be discussed)
  combinations <- combn(nodes,2) # all combinations of 2 distinct nodes in "nodes"
  n <- dim(combinations)[2]
  for (i in 1:n){
    observed <- nodes[nodes!=combinations[1,i] & nodes!=combinations[2,i]] # V'\{v,w}
    if (!is.element(combinations[1,i], nbr(sg, combinations[2,i])) & # check if there exists an edge already
        !bnlearn::dsep(dag,combinations[1,i],combinations[2,i])){ ### check if d-connected
      sg <- set.edge(sg, from = combinations[1,i], to = combinations[2,i]) ### undirected edge in case d-connected
    }
  }
  return(cpdag(sg)) ### to be discussed: return(cpdag(sg))
}

procedure2 <- extract_subgraph(dag_alarm, subgraph_nodes)
qgraph(procedure2, legend.cex = 0.3,
       asize=2,edge.color="black", vsize= 5)
```


Subsetting the dataset according to "subgraph_nodes" selection
```{r}
alarm_dfSubset <-as.data.frame(alarm_df[,subgraph_nodes])
```


```{r,echo = FALSE}
# procedure 3: forming the PAG from the observable and latent variables

## Step 1: Compute the true covariance matrix of the large graph G=(V,E).
correlation_matrix = cor(alarm_df)

#Step 2: subselect the observed variables and compute the  sub selected graph. 
observed_vars = subgraph_nodes
covMatrix_obs = correlation_matrix[observed_vars, observed_vars]

# Step 3: Provide nececssary information for the algorithm to work on. 

indepTest <- pcalg::gaussCItest
suffStat <- list(C = covMatrix_obs, n = 1000)

## Step 4: Apply the FCI algorithm to the selected vertices

# Use the FCI algorithm with the oracle correlation matrix
normal.pag <- pcalg:: fci(suffStat, indepTest, alpha = 0.05, labels = observed_vars, verbose=TRUE)

## Step 5: Plot the resulting PAG (Partial Ancestral Graph)
plot(normal.pag, 
     main = "Partial Ancestral Graph (PAG) for Subset of Nodes",
     vertex.label.cex = 0.8)

```


Applying constraint-based algorithms

```{r}
Res_stable=pc.stable(alarm_dfSubset)

Res_iamb=iamb(alarm_dfSubset)

Res_gs=gs(alarm_dfSubset)

Res_fiamb=fast.iamb(alarm_dfSubset)

Res_mmpc=mmpc(alarm_dfSubset)
```

Applying score-based algorithms
```{r}
Res_hc = hc(alarm_dfSubset)
Res_tabu = tabu(alarm_dfSubset)
```


Visualize resulting subgraph in the images below:

```{r, fig.show="hold", out.width="120%"}
par(mfrow=c(2, 4), mar=c(1, 1, 1, 1), oma=c(2, 2, 2, 2))


ig_proc1 <- as.igraph(procedure1)
ig_proc2 <- as.igraph(procedure2)
ig_pc <- as.igraph(Res_stable)
ig_iamb <- as.igraph(Res_iamb)
ig_gs <- as.igraph(Res_gs)
ig_fiamb <- as.igraph(Res_fiamb)
ig_mmpc <- as.igraph(Res_mmpc)

ig_hc <- as.igraph(Res_hc)
ig_tabu <- as.igraph(Res_tabu)

plot(ig_proc1, main = "proc1", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_proc2, main = "proc2", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_pc, main = "pc", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_iamb, main = "iamb", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_gs, main = "gs", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_fiamb, main = "fiamb", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_mmpc, main = "mmpc", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_hc, main = "hc", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)
plot(ig_tabu, main = "tabu", frame = T, layout=layout_with_fr, vertex.size=18,
     vertex.label.dist=4, vertex.color="red", edge.arrow.size=0.4, vertex.label.cex=1)

```


Selfdefined function for different measures
```{r}
measure = function(estim, true){
  result <- matrix(1,4)
  com <- bnlearn::compare(estim, true)
  shd <- bnlearn::shd(estim,true)
  result[1,] <- com$tp
  result[2,] <- com$fp
  result[3,] <- com$fn
  result[4,] <- shd
  rownames(result) <- c("true positives","false positives","false negatives","structural hamming distance")
  colnames(result) <- deparse(substitute(estim))
  return(result)
}
```

Metric evaluation for procedure 1
```{r}
measure(Res_stable, procedure1)
measure(Res_iamb, procedure1)
measure(Res_gs, procedure1)
measure(Res_fiamb, procedure1)
measure(Res_mmpc, procedure1)
measure(Res_hc, procedure1)
measure(Res_tabu, procedure1)
```

Metric evaluation for procedure 2
```{r}
measure(Res_stable, procedure2)
measure(Res_iamb, procedure2)
measure(Res_gs, procedure2)
measure(Res_fiamb, procedure2)
measure(Res_mmpc, procedure2)
measure(Res_hc, procedure2)
measure(Res_tabu, procedure2)
```