###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")