---
title: "Project"
author: "Haiyue Wang"
output:
  pdf_document: default
  bookdown::pdf_document2:
    latex_engine: lualatex
  tufte::tufte_handout:
    latex_engine: xelatex
  html_notebook: default
  html_document:
    df_print: paged
---

[Project on my blog](https://jet-w.github.io/unisa/2023SP5/AdvancedAnalyticTechniques2/project.html) \>Password: #Thuc

[Project on google collaboration](https://colab.research.google.com/drive/13V7I9g8k69c61NJzuwc80PwrOwM-Dz-0?usp=sharing)

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(gRain)
library(pcalg)
library(bnlearn)
```

## Read data

```{r message=FALSE, warning=FALSE}
data <- read_csv(
  "https://jet-w.github.io/data/unisa/AdvancedAnalytic2/project/BRCA_RNASeqv2_top50.csv", 
  col_names = TRUE
)
data.no.class <- data %>% select(-class)
```

## Task 1: Learn the global structure (CPDAG) [(4)]{style="color:red"}

**Requirement**: Use a [causal structure learning algorithm]{style="color:orange"} to find the [gene regulatory network]{style="color:orange"}, i.e. the network showing the interactions between genes, using the gene expression data. Explain how the algorithm works.

**Hints**: Please [exclude the class variable]{style="color:red"} in building the network **(4)**

------------------------------------------------------------------------

**Solution:** Using [PC algorithm]{style="color:red"} to learn the `CPTAG` , then using IDA for inferring causality.

### `1.1 Structure learning:`pc algorithm

```{r warning=FALSE}
pc.fit <- pc(
  suffStat = list(C = cor(data.no.class), n = nrow(data.no.class)), 
  indepTest = gaussCItest, 
  alpha=0.01, 
  labels = colnames(data.no.class)
  #labels = as.character(1:50)  #label node names
)

#data.frame(NUM=1:50, NAME=colnames(data.no.class))
plot(pc.fit, main = "CPDAG")
```

### 1.2 PC Algorithm Explanation

`PC algorithm` is an approach to learn the graph from a data set. The algorithm could be split into two parts. The first part is to learn the correlation among variables, which is also called learning skeleton. The second part is to find the direction of the relationship, which is also called as orientating the edges.

For this case,

------------------------------------------------------------------------

**Learning skeleton**

**Input:** Data set ***D***, significant level ***alpha***

**Output:** The undirected graph ***G*** with a set of edges ***E***

The first part, the alpha has been set ***alpha = 0.01***, ***depth = 0***, ***D*** is the data of all variables. Assume that all variables are correlated, the correlation set is ***E***.

-   **Repeat**
    -   **For** each edges in Graph ***G***, test the
        -   **If** number of nodes in ***E*** is greater or equal ***depth + 2*** jump to **NextLevel**, else continue.
            -   Test the independence of all nodes pairs from ***E*** given condition ***depth*** count combination of other variables from ***E***.
            -   **If** the independence exist,
                1.  remove the correlation from the ***E***.
                2.  Save the condition as the separation of the two nodes of the edge.
                3.  Update ***E***
            -   **End If**
        -   **End IF**
    -   **End For**
    -   **NextLevel** : [set ***depth = depth + 1***]{style="color:red"}
-   **Until**, the number of Node in ***E*** is less then ***depth + 2***.

3.  **Finally**, get the skeleton of the graph ***G***.

------------------------------------------------------------------------

**Orientating the edges**

**Input:** Skeleton ***G***, seperation sets ***S***

**Output:** CPDAG ***G\****

**for** all nonadjacent variables ***X***, ***Y*** with a common neighbor ***K*** do

**If** ***K*** does not belongs to separation set of the two nodes ***S(X,Y)*** then

Replace ***X-K-J*** in ***G*** by ***X-\>K\<-Y***

**end**

**end**

Next, orient as many other undirected edges as possible using the following rules:

1.  Orient ***X-Y*** into ***X-\>Y*** if exists ***Z-\>X***, ***Z*** and ***K*** are nonadjacent.
2.  Orient ***X-Y*** into ***X-\>Y***, if exists a chain ***X-\>Z-\>Y***.
3.  Orient ***X-Y*** into ***X-\>Y***, if exists two chains ***X-Z-\>Y*** and ***X-A-\>Y***, and ***Z*** and ***A*** are nonadjacent.

Finally, get a CPDAG ***G\****

## Task 2: Causal Effects (IDA)[(4)]{style="color:red"}

**Requirement**: `EBF1` is an important gene that is involved in many biological processes leading to cancer. Find the [top 10]{style="color:orange"} other genes that have [strong causal effects]{style="color:orange"} on `EBF1` using a causal inference algorithm.

**Hints**

1.  [Exclude]{style="color:red"} the class variable in building the network

2.  If there are [multiple possible causal effects]{style="color:red"} between the cause and the effect, we can use the [**minimum** of the **absolute values**]{style="color:red"} (of the causal effects) as the final result

3.  The causal effects are normally [ranked based on their absolute values]{style="color:red"}.

------------------------------------------------------------------------

**Solution**: Using `ida` to calculate the causal effects of all other variables on `EBF1` based on the graph built from task 1, then sort the final result base on the values come from `ida` algorithm.

### Causal Effects on EBF1 from Other Variables

```{r warning=FALSE}
# Get gene EBF1 index
EBF1_idx <- match("EBF1", names(data.no.class))
CausalOnEBF1 <- data.frame(
  causality = unlist(
    lapply(
      (1:50)[-EBF1_idx], 
      function(idx){
        min(
          abs(
            idaFast(
              idx,
              EBF1_idx,
              cov(data.no.class),
              pc.fit@graph)
          )
        )
      }
    )
  ),
  variable  = names(data.no.class)[-EBF1_idx]
) 
CausalOnEBF1 %>% 
  arrange(across(causality, desc))
```

```{r}
(CausalOnEBF1 %>% 
  arrange(across(causality, desc)))$variable[1:10]
```

According to the result, it could easily get the top 10 genes have strongest causal effects on `EBF1` are `FXYD1`, `ABCA10`, `TMEM220`, `ARHGAP20`, `FIGF`, `KLHL29`, `GPIHBP1`, `TMEM132C`, `RDH5`, `ABCA9`.

## Task 3: Local Causal Structure & Markov blanket [(4)]{style="color:red"}

**Requirement**: Use a [local causal structure learning algorithm]{style="color:orange"} to find genes in the Markov blanket of `ABCA9` from data. Explain how the algorithm works. **(4)**

------------------------------------------------------------------------

**Solution:** We could use local structure learning algorithm `IAMB` to get the Markov blanket of `ABCA9` from the data

### 3.1 Calculating the Markov Blanket

```{r warning=FALSE}
data.num <- data %>% select(-class)
data.num$class <- ifelse(data$class == 'C', 1, 0)
ABCA9.mb <- learn.mb(
    data.frame(data.num),
    "ABCA9",
    method = "iamb", 
    alpha = 0.01
)
ABCA9.mb
```

According to the result above, the Markov Blanket of `ABCA9` has 23 nodes.

### 3.2 Explanation of IAMB

The **IAMB** is an abbreviation for ***Incremental Association Markov Blanket***, the algorithm could be separated into two phases, the ***Growing pahse*** and ***Shrinking phase***. Details for the two phases are below:

**CMI:** Conditional mutual information

**Input:** dataset ***D***; target ***T***

**Output:** ***MB(T)***

------------------------------------------------------------------------

**Growing Phase:**

**Repeat** till ***MB(T)*** does not change

-   Find the node ***X*** from dataset ***D*** [exclude all the nodes in ***MB(T)*** and ***T***]{style="color:red"} that has the maximum **CMI**

-   **IF** ***X*** independence with ***T*** given ***MB(T)***, Then

    -   [Add]{style="color:red"} ***X*** to ***MB(T)***

-   **End IF**

------------------------------------------------------------------------

**Shrinking Phase:**

**For** each node ***X*** from ***MB(T)***

-   **IF** ***X*** independence with ***T*** given ***MB(T)*** [exclude]{style="color:red"} ***X***, **Then**

    -   [Remove]{style="color:red"} ***X*** from ***MB(T)***

-   **End IF**

**End For**

Finally, get the final ***MB(T)***.

------------------------------------------------------------------------

**References:**

[An Improved IAMB Algorithm for Markov Blanket Discovery](http://www.jcomputers.us/vol5/jcp0511-18.pdf)

[Discovering Markov Blankets: Finding Independencies Among Variables](https://cseweb.ucsd.edu//~elkan/254/Verma.pdf)

## Task 4: Discrete the dataset

**Requirements:** [Discretise]{style="color:red"} the dataset to [binary]{style="color:red"} using the [average expression of **ALL genes**]{style="color:red"} as the threshold. The discretised dataset will be used in the following questions.

------------------------------------------------------------------------

**Solution:**

**Step 1**: Calculating the mean

**Step 2**: Discrete the data (1: \> mean 0: \<= mean)

Because of the `pcSelect` method [only support numeric variables]{style="color:red"}, so the discrete variables need to be replaced with 1 and 0 according to **step 2**.

```{r setup, warning=FALSE}
# The mean of each gene
# mean.val <- as.data.frame(apply(data.no.class, 2, mean))
# The mean of All genes
mean.val <- mean(apply(data.no.class, 2, mean))
names <- colnames(data.no.class)

data.binary <- as.data.frame(
  sapply(
    colnames(data.no.class), 
    function(x) ifelse(data.no.class[,x] >mean.val, 1, 0)
))
data.binary$class <- ifelse(data$class == 'C', 1, 0)


#data.binary.c <- as.data.frame(
#  sapply(
#    colnames(data.no.class), 
#    function(x) ifelse(data.no.class[,x] >mean.val, 'T', 'F')
#))
#data.binary.c$class <- ifelse(data$class == 'C', 'T', 'F')


```

## Task 5: PC-Simple [(6)]{style="color:red"}

**Requirements:**

1.  Use [PC-simple algorithm (pcSelect)]{style="color:orange"} to find the parent and children set of the class variable.
2.  Explain how `PC-simple` works.
3.  Evaluate the accuracy of the Naïve Bayes classification on the dataset in the following cases:
    1.  Use all features (genes) in the dataset

    2.  Use only the features (genes) in the parent and children set of the class variable

    3.  Compare the accuracy of the models in the two cases using 10-fold cross validation.

[References](https://www.analyticsvidhya.com/blog/2021/03/introduction-to-k-fold-cross-validation-in-r/)

### 5.1 Find the parents and children

```{r warning=FALSE}
class.pc <- pcSelect(
  data.binary %>% select(class),
  data.binary %>% select(-class),
  alpha = 0.01
)
class.pc <- data.frame(ispc = class.pc$G, zmin = class.pc$zMin)
class.pc[order(class.pc$zmin, decreasing=TRUE),]
```

```{r warning=FALSE}
rownames(class.pc[class.pc$ispc == TRUE,])
```

According to the result above, it could easily found that the parents and children of `class` variable are `FIGF`, `ARHGAP20`, `CD300LG`, `KLHL29`, `CXCL2`, `ATP1A2`, `MAMDC2`, `TMEM220`, `SCARA5`, `ATOH8`.

### 5.2 Explanation of PC-Simple

PC-Simple is an algorithm to find the parents and children of a target node via conditional independence tests base on a threshold ***alpha***,

**Input:** Dataset ***D*** consist of set of predictor variables ***X*** and target variable ***Z***; a significant levle ***alpha*** for conditional independence test.

**Output:** The parents and children set ***PC*** of target ***Z***

**For** this case, the ***alpha*** has been set ***0.01***, ***k = 0*** , ***PC(k)*** equals all other variables.

**Repeat** if count of ***PC*** is greater than ***k***

-   ***k = k+1***

-   ***PC(k) = PC(k-1)***

-   **For** each node ***X*** from ***PC(k-1)***

    -   **For** each combination nodes ***S*** from ***PC(k-1)*** excludes ***X*** and count of ***S*** equals ***k-1***

        -   **IF** ***X*** and ***Z*** are independent given ***S*** at significance level ***alpha***, Then

            -   remove ***X*** from ***PC(K)***

        -   **End IF**

    -   **End For**

-   ***End For***

The ***PC*** is the final result.

### 5.3 Naïve Bayes classification

The Naive Bayes classifier is using `caret` package to train model.

```{r message=FALSE, warning=FALSE}
library(caret)
```

#### 5.3.1 Naive Bayes classification with all features

-   The code below using all features as input features for training.

```{r warning=FALSE}
set.seed(100)
trctrl <- trainControl(method = "cv", number = 10, savePredictions=TRUE)
ma <- train(
  factor(class) ~., 
  data = data.binary, 
  method = "naive_bayes", 
  trControl=trctrl, 
  tuneLength = 0
)
ma
```

#### 5.3.2 Naive Bayes classification with PC set

-   The code below using parents and children features of `class` feature as input features for training.

```{r warning=FALSE}
data.binary.related <- data.binary[,append(rownames(class.pc[class.pc$ispc == TRUE,]), "class")]

mr <- train(
  factor(class) ~., 
  data = data.binary.related, 
  method = "naive_bayes", 
  trControl=trctrl, 
  tuneLength = 0
)
mr
```

#### 5.3.3 Comparasion between the two models

To compare the two models using several indicators, for this model we just using ***accuracy***, ***precision*** and ***recall*** to compare the two models. For easy compare the two model, we name the model trained with all features as ***ma***, name another one as ***mr***.

```{r}
#confusionMatrix(ma)
predicted.ma <- predict(ma, data.binary)
confusionMatrix(as.factor(data.binary$class), as.factor(predicted.ma))
```

```{r warning=FALSE}
predicted.mr <- predict(mr, data.binary)
confusionMatrix(as.factor(data.binary$class), as.factor(predicted.mr))
```

According to the confusion matrix, we could get indicators table like below.

| indicators                | ***ma***                    | ***mr***                    |
|-----------------------|-------------------------|-------------------------|
| **accuracy**              | 0.9761                      | [0.9843]{style="color:red"} |
| **precision(1)** : cancer | [0.9982]{style="color:red"} | 0.9809                      |
| **precision(0)**: normal  | 0.7589                      | [1]{style="color:red"}      |
| **recall(1)**: cancer     | 0.9760                      | [1]{style="color:red"}      |
| **recall(0)**: normal     | [0.9770]{style="color:red"} | 0.8421                      |

The precision stands for the accuracy of prediction cases, the recall represents the accuracy of actual cases that has been recognized. According to the table above, the overall accuracy of ***mr*** is better than ***ma***. For cancer cases, the ***ma*** perform better than ***mr*** on prediction, but for recall value of ***mr*** is better than ***ma***. For normal cases, the ***mr*** works better than ***ma*** on precision, but the ***ma*** perform better on recall than ***mr***.

## Task 6: **Calculating based on specified DAG**

### 6.1 a) Construct the conditional probability tables for the Bayesian network based on data. [**(3)**]{style="color:red"}

#### 6.1.1 Construct using `cptable`

For constructing the net, using the ***T*** instead of ***1*** and using ***F*** instead of ***0***.

```{r warning=FALSE}

data.graph <- data.binary %>% select(BTNL9,CD300LG,class,IGSF10,ABCA9)

yn <- c('T','F')

B     <- cptable(~BTNL9, 
                 values= (data.graph %>% 
                            select(BTNL9) %>%
                            group_by(BTNL9) %>%
                            count() %>%
                            arrange(across(BTNL9, desc)))$n,
                 levels=yn)
CD.B  <- cptable(~CD300LG|BTNL9, 
                 values=(data.graph %>% 
                            select(CD300LG, BTNL9) %>%
                            group_by(CD300LG, BTNL9) %>%
                            count() %>%
                            arrange(across(BTNL9, desc),
                                    across(CD300LG, desc)))$n,
                 levels=yn)
c.CD  <- cptable(~class|CD300LG, 
                 values=(data.graph %>% 
                            select(class, CD300LG) %>%
                            group_by(class, CD300LG) %>%
                            count() %>%
                            arrange(across(CD300LG, desc),
                                    across(class, desc)))$n, 
                 levels=yn)
I.c   <- cptable(~IGSF10|class, 
                 values=(data.graph %>% 
                            select(IGSF10, class) %>%
                            group_by(IGSF10, class) %>%
                            count() %>%
                            arrange(across(class, desc), 
                                    across(IGSF10, desc)
                                    ))$n, 
                 levels=yn)

AB.B_I<- cptable(~ABCA9|BTNL9:IGSF10,
                 values=(data.graph %>% 
                            select(ABCA9, BTNL9,IGSF10) %>%
                            group_by(ABCA9, BTNL9,IGSF10) %>%
                            count() %>%
                            arrange( 
                              across(IGSF10, desc), 
                              across(BTNL9, desc),
                              across(ABCA9, desc)))$n,
                 levels=yn)

plist <- compileCPT(list(B, CD.B, c.CD, I.c, AB.B_I))
plist
net=grain(plist) 
```

```{r}
plot(net$dag)
```

#### 6.1.2 Construct using `bnlearn`

-   constructing the net structure

```{r}
bn.dag = model2network("[BTNL9][CD300LG|BTNL9][ABCA9|BTNL9:IGSF10][class|CD300LG][IGSF10|class]")
graphviz.plot(bn.dag)
```

-   learn parameter from data

```{r}
bn.fitted <- bn.fit(
  bn.dag,
  data.binary %>% select(BTNL9, CD300LG, class, IGSF10, ABCA9)
) 
bn.fitted
```

### 6.2 b) Estimate the probability of the four genes in the network having high expression levels. [**(2)**]{style="color:red"}

This question aims to calculate the ***joint probability*** of the four genes in the network for each value of the four variables equal ***T***. It could be expressed with the formula ***P(BTNL9=T, CD300LG=T, IGSF10=T, ABCA9=T)***.

#### 6.2.1 Method 1

-   method 1 to query the probability

```{r}
querygrain(net, nodes=c("BTNL9", "CD300LG", "IGSF10", "ABCA9"), type="joint")
```

According to the table above, it could get the ***P(BTNL9=T, CD300LG=T, IGSF10=T, ABCA9=T)=[0.073736]{style="color:red"}***

#### 6.2.2 Method 2

-   method 2 to query the probability

```{r}
joint_pb <- setEvidence(
  net, 
  evidence=list(BTNL9="T", CD300LG="T",  IGSF10="T", ABCA9="T")
)
pEvidence(joint_pb)
```

According to the table above, it could get the same result with method 1.

### 6.3 c) Estimate the probability of having cancer when the expression level of `CD300LG` is high and the expression level of `BTNL9` is low. [**(2)**]{style="color:red"}

This question actually ask us to calculate the conditional probability ***P(class=T\| CD300LG=T, BTNL9=F)***, here I will use `cpquery` method for get the conditional probability.

```{r}
querygrain(net, nodes=c("class","CD300LG","BTNL9"), type="conditional")
```

So the final result ***P(class=T\|CD300LG=T,BTNL9=F)*** = [***0.2585034***]{style="color:red"}

### 6.4 Prove the result in c) mathematically. [**(2)**]{style="color:red"}

```{r}
data.graph %>% 
  select(class, CD300LG,BTNL9) %>%
  group_by(class, CD300LG,BTNL9) %>%
  count() %>%
  arrange( 
    across(class, desc), 
    across(CD300LG, desc),
    across(BTNL9, desc))
#plot(net$dag)

```

According to the graph, BTNL9 is the parent of CD300LG, so

***P(class=T\| BTNL9=F,CD300LG=T)***

***= P(class=T\|CD300LG=T) =***

=**(32+6)/(32+107+6+2)**

= [***0.2585034***]{style="color:red"}

### 6.5 Given we know the value of `CD300LG`, is the "class" conditionally independent of `ABCA9` and why? [**(3)**]{style="color:red"}

```{r}
plot(net$dag)
```

**Anwser:** [No]{style="color:red"}

**Explanation:** According to Markov condition, every node in a Bayesian network is conditionally independent of its nondescendants, given its parents. So the parent `CD300LG` of `class` is given, `ABCA9` is the descendant of `class` variable, so the `class` is not conditionally independent of `ABCA9`.

## References

Algorithm:

<http://www.sci-princess.info/wp-content/uploads/Causal-Graphs-and-the-PC-Algorithm.pdf>

<https://pooyanjamshidi.github.io/csce580/lectures/CSCE580-GuestLecture--BNLearning.pdf>

<https://arxiv.org/pdf/0908.3817.pdf>

<https://www.bnlearn.com/about/slides/slides-useRconf13.pdf>

pcalg:

<https://stat.ethz.ch/Manuscripts/buhlmann/pcalg-software.pdf>

<https://cran.r-project.org/web/packages/pcalg/pcalg.pdf>

<https://cran.r-project.org/web/packages/pcalg/vignettes/vignette2018.pdf>

bnlearn:

<https://www.bnlearn.com/examples/graphviz-plot/>

<https://www.bnlearn.com/documentation/man/cpquery.html>

<https://rdrr.io/github/vspinu/bnlearn/man/cpquery.html>

<https://dipartimenti.unicatt.it/scienze-statistiche-23-25-1-17ScutariSlides.pdf>

Graphviz

<https://rdrr.io/bioc/Rgraphviz/man/GraphvizAttributes.html#:~:text=Font%20size%2C%20in%20points%2C%20for,Label%20for%20the%20graph.>

<https://www.cs.cmu.edu/afs/cs/project/jair/pub/volume18/acid03a-html/node2.html>

Chinese Sample:

<https://www.cnblogs.com/payton/articles/4608383.html>
