###BNLearn Conditional Independency Test###

install.packages("bnlearn")
library(bnlearn)
#####
#To upload the dataset
#####
#However, this is not necessary to create a DAG
D50S9V2let <- read.table("C:/Users/biology/Documents/MS DS Project/D50S9v2let.txt")
#####
#####
###To Create the DAG of the Sparse Model
#####
dag50SV2<-empty.graph(nodes=c("BP",	"Catechol",	"CO",	"ExpCO2",	"HRBP",	"Intubation",	"LVFailure",	"VentAlv",
                "VentTube"))
arc.set<-matrix(c("LVFailure","CO",
                  "CO", "BP",
                  "Catechol", "CO",
                  "Catechol", "HRBP",
                  "VentAlv", "Catechol",
                  "VentAlv", "ExpCO2",
                  "Intubation", "VentAlv",
                  "VentTube", "VentAlv",
                  "Intubation", "Catechol",
                  "VentTube", "ExpCO2",
                  "Intubation", "ExpCO2"),
          byrow = TRUE, ncol=2,
          dimnames = list(NULL,c("from", "to")))

arcs(dag50SV2)<-arc.set
#####
#####
#To visually plot the DAG
#####
install.packages("Rgraphviz")
library(Rgraphviz)


graphviz.plot(dag50SV2)
#####
#####
#To determine if the pairs are d-separated or not
#####
dsep(dag50SV2, x="CO", y="BP", z="LVFailure")
dsep(dag50SV2, x="CO", y="BP", z=c("LVFailure","Catechol"))
dsep(dag50SV2, x="CO", y="BP")
dsep(dag50SV2, x="LVFailure", y="VentTube")

#No parents
#Top 8
dsep(dag50SV2, x="LVFailure", y="VentAlv")
dsep(dag50SV2, x="Intubation", y="LVFailure")
dsep(dag50SV2, x="HRBP", y="LVFailure")
dsep(dag50SV2, x="ExpCO2", y="LVFailure")
dsep(dag50SV2, x="BP", y="VentTube")
dsep(dag50SV2, x="LVFailure", y="VentTube")
dsep(dag50SV2, x="LVFailure", y="VentTube")
dsep(dag50SV2, x="Intubation", y="VentTube")
dsep(dag50SV2, x="BP", y="LVFailure")
#No parents
#Bottom 10
dsep(dag50SV2, x="Catechol", y="VentAlv")
dsep(dag50SV2, x="BP", y="CO")
dsep(dag50SV2, x="Catechol", y="CO")
dsep(dag50SV2, x="ExpCO2", y="Intubation")
dsep(dag50SV2, x="Intubation", y="VentAlv")
dsep(dag50SV2, x="CO", y="HRBP")
dsep(dag50SV2, x="ExpCO2", y="VentTube")
dsep(dag50SV2, x="Catechol", y="HRBP")
dsep(dag50SV2, x="VentAlv", y="VentTube")
dsep(dag50SV2, x="ExpCO2", y="VentAlv")
#One Parent
#Top 10
dsep(dag50SV2, x="CO", y="VentTube", z="Catechol")
dsep(dag50SV2, x="HRBP", y="VentAlv", z="Catechol")
dsep(dag50SV2, x="CO", y="VentTube", z="HRBP")
dsep(dag50SV2, x="Catechol", y="VentTube", z="ExpCO2")
dsep(dag50SV2, x="LVFailure", y="VentAlv", z="Intubation")
dsep(dag50SV2, x="Intubation", y="LVFailure", z="VentAlv")
dsep(dag50SV2, x="BP", y="VentTube", z="ExpCO2")
dsep(dag50SV2, x="BP", y="VentTube", z="HRBP")
dsep(dag50SV2, x="CO", y="VentTube", z="ExpCO2")
dsep(dag50SV2, x="HRBP", y="VentTube", z="ExpCO2")
#One Parent
#Bottom 10
dsep(dag50SV2, x="VentAlv", y="VentTube", z="Catechol")
dsep(dag50SV2, x="VentAlv", y="VentTube", z="LVFailure")
dsep(dag50SV2, x="VentAlv", y="VentTube", z="Intubation")
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z="VentTube")
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z="Catechol")
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z="BP")
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z="CO")
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z="HRBP")
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z="Intubation")
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z="LVFailure")
#Two Parents
#Top 10
dsep(dag50SV2, x="Catechol", y="VentTube", z=c("ExpCO2","HRBP"))
dsep(dag50SV2, x="HRBP", y="VentTube", z=c("ExpCO2","VentAlv"))
dsep(dag50SV2, x="LVFailure", y="VentAlv", z=c("ExpCO2","VentTube"))
dsep(dag50SV2, x="ExpCO2", y="LVFailure", z=c("VentAlv","VentTube"))
dsep(dag50SV2, x="BP", y="Intubation", z=c("ExpCO2","VentAlv"))
dsep(dag50SV2, x="CO", y="VentTube", z=c("ExpCO2","VentAlv"))
dsep(dag50SV2, x="BP", y="VentTube", z=c("ExpCO2","VentAlv"))
dsep(dag50SV2, x="HRBP", y="VentTube", z=c("ExpCO2","Intubation"))
dsep(dag50SV2, x="HRBP", y="Intubation", z=c("ExpCO2","VentTube"))
dsep(dag50SV2, x="HRBP", y="Intubation", z=c("ExpCO2","VentAlv"))
#Two Parents
#Bottom 10
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("Catechol","HRBP"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("Catechol","LVFailure"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("BP","Intubation"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("Intubation","LVFailure"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("CO","HRBP"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("BP","CO"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("BP","HRBP"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("BP","LVFailure"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("CO","LVFailure"))
dsep(dag50SV2, x="ExpCO2", y="VentAlv", z=c("HRBP","LVFailure"))
#####
#####
###To Create the DAG of the Close Model
#####
dag50CV2<-empty.graph(nodes=c("Anaphylaxis",	"ArtCO2",	"Catechol",	"ExpCO2",	
                              "HR",	"InsuffAnesth",	"SaO2",	"TPR",	"VentLung"))

arcclose.set<-matrix(c("Anaphylaxis","TPR",
                  "TPR", "Catechol",
                  "Catechol", "HR",
                  "InsuffAnesth", "Catechol",
                  "SaO2", "Catechol",
                  "ArtCO2", "Catechol",
                  "ArtCO2", "ExpCO2",
                  "VentLung", "ExpCO2",
                  "VentLung", "ArtCO2"),
                byrow = TRUE, ncol=2,
                dimnames = list(NULL,c("from", "to")))

arcs(dag50CV2)<-arcclose.set
#####
#####
#To visually plot the DAG
#####
install.packages("Rgraphviz")
library(Rgraphviz)


graphviz.plot(dag50CV2)
#####
#####
#To determine if the pairs are d-separated or not
#####
#No parents
#Top 10
dsep(dag50CV2, x="InsuffAnesth", y="TPR")
dsep(dag50CV2, x="ExpCO2", y="TPR")
dsep(dag50CV2, x="Anaphylaxis", y="HR")
dsep(dag50CV2, x="Anaphylaxis", y="ExpCO2")
dsep(dag50CV2, x="InsuffAnesth", y="SaO2")
dsep(dag50CV2, x="Anaphylaxis", y="VentLung")
dsep(dag50CV2, x="InsuffAnesth", y="VentLung")
dsep(dag50CV2, x="ArtCO2", y="TPR")
dsep(dag50CV2, x="TPR", y="VentLung")
dsep(dag50CV2, x="Anaphylaxis", y="InsuffAnesth")
#No parents
#Bottom 10
dsep(dag50CV2, x="HR", y="SaO2")
dsep(dag50CV2, x="Catechol", y="TPR")
dsep(dag50CV2, x="Catechol", y="SaO2")
dsep(dag50CV2, x="ExpCO2", y="SaO2")
dsep(dag50CV2, x="SaO2", y="VentLung")
dsep(dag50CV2, x="ArtCO2", y="SaO2")
dsep(dag50CV2, x="Catechol", y="HR")
dsep(dag50CV2, x="ArtCO2", y="VentLung")
dsep(dag50CV2, x="ExpCO2", y="VentLung")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2")
#One Parent
#Top 10
dsep(dag50CV2, x="Anaphylaxis", y="ExpCO2", z="VentLung")
dsep(dag50CV2, x="Anaphylaxis", y="VentLung", z="ExpCO2")
dsep(dag50CV2, x="Anaphylaxis", y="VentLung", z="TPR")
dsep(dag50CV2, x="InsuffAnesth", y="TPR", z="Anaphylaxis")
dsep(dag50CV2, x="InsuffAnesth", y="SaO2", z="VentLung")
dsep(dag50CV2, x="ExpCO2", y="TPR", z="VentLung")
dsep(dag50CV2, x="Anaphylaxis", y="VentLung", z="TPR")
dsep(dag50CV2, x="Anaphylaxis", y="HR", z="Catechol")
dsep(dag50CV2, x="HR", y="VentLung", z="Catechol")
dsep(dag50CV2, x="InsuffAnesth", y="VentLung", z="SaO2")
#One Parent
#Bottom 10
dsep(dag50CV2, x="ExpCO2", y="VentLung", z="TPR")
dsep(dag50CV2, x="ExpCO2", y="VentLung", z="Anaphylaxis")
dsep(dag50CV2, x="ExpCO2", y="VentLung", z="InsuffAnesth")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z="VentLung")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z="SaO2")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z="Anaphylaxis")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z="Catechol")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z="HR")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z="InsuffAnesth")
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z="TPR")
#Two Parents
#Top 10
dsep(dag50CV2, x="HR", y="VentLung", z=c("ArtCO2","ExpCO2"))
dsep(dag50CV2, x="Anaphylaxis", y="VentLung", z=c("ExpCO2","SaO2"))
dsep(dag50CV2, x="Anaphylaxis", y="VentLung", z=c("ExpCO2","TPR"))
dsep(dag50CV2, x="ArtCO2", y="HR", z=c("ExpCO2","VentLung"))
dsep(dag50CV2, x="Anaphylaxis", y="ExpCO2", z=c("TPR","VentLung"))
dsep(dag50CV2, x="Anaphylaxis", y="ExpCO2", z=c("ArtCO2","VentLung"))
dsep(dag50CV2, x="Anaphylaxis", y="ArtCO2", z=c("ExpCO2","SaO2"))
dsep(dag50CV2, x="Anaphylaxis", y="ExpCO2", z=c("HR","VentLung"))
dsep(dag50CV2, x="Anaphylaxis", y="ExpCO2", z=c("SaO2","VentLung"))
dsep(dag50CV2, x="Anaphylaxis", y="VentLung", z=c("ExpCO2","HR"))
#Two Parents
#Bottom 10
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("Anaphylaxis","Catechol"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("Anaphylaxis","HR"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("Anaphylaxis","InsuffAnesth"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("Anaphylaxis","TPR"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("Catechol","HR"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("Catechol","InsuffAnesth"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("Catechol","TPR"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("HR","InsuffAnesth"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("HR","TPR"))
dsep(dag50CV2, x="ArtCO2", y="ExpCO2", z=c("InsuffAnesth","TPR"))
#####

#####
###ROC Curve###
#####
#Need datasets with the predictor (p-value) and 0/1 corretness fromt he BNLearn output

#1K Sparse 1 Parent
D1KS9V <- read.csv("C:/Users/biology/Documents/MS DS Project/D1KS9v2ROC.csv")

install.packages("pROC")
library(pROC)
pROC_1KS<-roc(D1KS9V$label,D1KS9V$pvalue, smoothed=TRUE,
              ci=TRUE, ci.alpha=0.9, stratified=FALSE,
              plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
              print.auc=TRUE, show.thres=TRUE)

sens.ci<-ci.se(pROC_1KS)
plot(sens.ci, type="shape", col="lightblue")
plot(sens.ci, type="bars")
#1K Sparse No Parents
D1KS9VNoPar <- read.csv("C:/Users/biology/Documents/MS DS Project/D1KS9vNoParROC.csv")

install.packages("pROC")
library(pROC)
pROC_1KSNP<-roc(D1KS9VNoPar$label,D1KS9VNoPar$pvalue, smoothed=TRUE,
              ci=TRUE, ci.alpha=0.9, stratified=FALSE,
              plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
              print.auc=TRUE, show.thres=TRUE)

sens.ciNP<-ci.se(pROC_1KSNP)
plot(sens.ciNP, type="shape", col="lightblue")
plot(sens.ciNP, type="bars")
#1K Sparse TWO Parents
D1KS9V2Par <- read.csv("C:/Users/biology/Documents/MS DS Project/D1KS9v2ParROC.csv")

install.packages("pROC")
library(pROC)
pROC_1KS2P<-roc(D1KS9V2Par$label,D1KS9V2Par$pvalue, smoothed=TRUE,
                ci=TRUE, ci.alpha=0.9, stratified=FALSE,
                plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
                print.auc=TRUE, show.thres=TRUE)

sens.ci2P<-ci.se(pROC_1KS2P)
plot(sens.ci2P, type="shape", col="lightblue")
plot(sens.ci2P, type="bars")
###

#1K Close 1 Parent
D1KC9V <- read.csv("C:/Users/biology/Documents/MS DS Project/D1KC9vROC.csv")

install.packages("pROC")
library(pROC)
pROC_1KC<-roc(D1KC9V$label,D1KC9V$pvalue, smoothed=TRUE,
              ci=TRUE, ci.alpha=0.9, stratified=FALSE,
              plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
              print.auc=TRUE, show.thres=TRUE)

sens.ciC<-ci.se(pROC_1KC)
plot(sens.ciC, type="shape", col="lightblue")
plot(sens.ciC, type="bars")
#1K Close No Parents
D1KC9VNoPar <- read.csv("C:/Users/biology/Documents/MS DS Project/D1KC9vNoParROC.csv")

install.packages("pROC")
library(pROC)
pROC_1KCNP<-roc(D1KC9VNoPar$label,D1KC9VNoPar$pvalue, smoothed=TRUE,
                ci=TRUE, ci.alpha=0.9, stratified=FALSE,
                plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
                print.auc=TRUE, show.thres=TRUE)

sens.ciCNP<-ci.se(pROC_1KCNP)
plot(sens.ciCNP, type="shape", col="lightblue")
plot(sens.ciCNP, type="bars")
#1K Close TWO Parents
D1KC9V2Par <- read.csv("C:/Users/biology/Documents/MS DS Project/D1KC9v2ParROC.csv")

install.packages("pROC")
library(pROC)
pROC_1KC2P<-roc(D1KC9V2Par$label,D1KC9V2Par$pvalue, smoothed=TRUE,
                ci=TRUE, ci.alpha=0.9, stratified=FALSE,
                plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
                print.auc=TRUE, show.thres=TRUE)

sens.ciC2P<-ci.se(pROC_1KC2P)
plot(sens.ciC2P, type="shape", col="lightblue")
plot(sens.ciC2P, type="bars")