-
Notifications
You must be signed in to change notification settings - Fork 0
/
sachs_script_discrete.R
160 lines (146 loc) · 4.5 KB
/
sachs_script_discrete.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
library(bnlearn)
library(fmsb)
library(purrr)
source("./utils.R")
# Load Data
load("./data/sachs/sachs.rda") # Loads true network as `bn`
sachs.df <- read.table("./data/sachs/sachs.data.txt", header=TRUE)
sachs.dis <- discretize(
sachs.df, method = "hartemink",breaks = 5, ibreaks = 100, idisc = "quantile"
)
sachs.bn <- bn.net(bn)
num.bootstraps <- 50
num.permutations <- 100
algo.params <- list(
'inter.iamb.mi'=list(
'algorithm'='inter.iamb',
'algorithm.args'=list(alpha=0.05, test='mi')
),
'inter.iamb.mi.sh'=list(
'algorithm'='inter.iamb',
'algorithm.args'=list(alpha=0.05, test='mi-sh')
),
'inter.iamb.x2'=list(
'algorithm'='inter.iamb',
'algorithm.args'=list(alpha=0.05, test='x2')
),
'h2pc.mi' = list(
'algorithm'='h2pc',
'algorithm.args'=list(
restrict.args=list(alpha=0.15, test='mi'),
maximize.args=list(perturb=4, restart=25, score='bde')
)
),
'h2pc.mi.sh' = list(
'algorithm'='h2pc',
'algorithm.args'=list(
restrict.args=list(alpha=0.15, test='mi-sh'),
maximize.args=list(perturb=4, restart=25, score='bde')
)
),
'h2pc.x2' = list(
'algorithm'='h2pc',
'algorithm.args'=list(
restrict.args=list(alpha=0.15, test='x2'),
maximize.args=list(perturb=4, restart=25, score='bde')
)
),
'hc.0'=list(
'algorithm'='hc',
'algorithm.args'=list(score='bde', restart=0, perturb=4)
),
'hc.10'=list(
'algorithm'='hc',
'algorithm.args'=list(score='bde', restart=10, perturb=4)
),
'tabu.5'=list(
'algorithm'='tabu',
'algorithm.args'=list(score='bde', tabu=5)
),
'tabu.10'=list(
'algorithm'='tabu',
'algorithm.args'=list(score='bde', tabu=10)
)
# 'tabu.15'=list(
# 'algorithm'='tabu',
# 'algorithm.args'=list(score='bde', tabu=15)
# )
)
averaged.networks <- vector(length(algo.params), mode="list")
for (i in 1:length(algo.params)){
algo.name <- names(algo.params)[i]
algo.args <- algo.params[[i]]$algorithm.args
cat(sprintf(" Time: %s. Algorithm: %s\n", Sys.time(), algo.name))
boot.graph <- boot.strength(
sachs.dis,
R=num.bootstraps,
m=size.bootstrap,
algorithm=algo.params[[i]]$algorithm,
algorithm.args=algo.args
)
averaged.networks[[i]] <- averaged.network(boot.graph)
}
# For each network
results <- array(numeric(), c(0, 5))
colnames(results) <- c(
'algorithm',
'full.shd',
'full.hd',
'full.fdr',
'full.sensitivity'
)
for (i in seq(1, length(averaged.networks))){
learned.graph <- averaged.networks[[i]]
relevant.nodes = nodes(learned.graph)[sapply(nodes(learned.graph), degree, object = learned.graph) > 0]
sachs.subgraph <- subgraph(sachs.bn, relevant.nodes); learned.subgraph <- subgraph(learned.graph, relevant.nodes)
full.shd <- bnlearn::shd(cpdag(learned.subgraph), cpdag(sachs.subgraph))
full.hd <- bnlearn::hamming(cpdag(learned.subgraph), cpdag(sachs.subgraph))
full.performance <- calculate.performance.statistics(cpdag(learned.subgraph), cpdag(sachs.subgraph))
full.fdr <- full.performance$fdr; full.sensitivity <- full.performance$sensitivity
print(c(
names(algo.params)[i],
full.shd,
full.hd,
full.fdr,
full.sensitivity
))
results <- rbind(
results,
c(
names(algo.params)[i],
full.shd,
full.hd,
full.fdr,
full.sensitivity
)
)
}
write.csv(results, 'data/discrete_sachs_results.csv', row.names=FALSE)
library(ggplot2)
library(ggradar)
library(scales)
library(data.table)
sachs.results <- read.csv('data/discrete_sachs_results.csv', row.names='algorithm')
sachs.results$norm.shd <- sachs.results$full.shd / 20
sachs.results$norm.hd <- sachs.results$full.hd / 20
sachs.results <- sachs.results[, c("norm.shd", "norm.hd", "full.fdr", "full.sensitivity")]
sachs.sub <- sachs.results[rownames(sachs.results) %like% id, ]
sachs.results.t <- t(sachs.results)
headers <- c(colnames(sachs.results.t))
sachs.results.t <- cbind(rownames(sachs.results.t), data.frame(sachs.results.t, row.names=NULL))
sachs.results.t[,1] = c(
'Normalised SHD', 'Normalised HD', 'FDR', 'Sensitivity'
)
colnames(sachs.results.t) <- c('metric', headers)
(radar.plot <- ggradar(
sachs.results.t, axis.label.size = 3.5, grid.label.size = 6,
legend.text.size = 10, grid.max=0.6, grid.mid=0.3,
values.radar = c('0', '0.3', '0.6')
))
ggsave(
"./report/figures/radar_dis.pdf",
radar.plot,
width=9.5,
height=5.3,
device='pdf'
)