# This number of MCMC samples is for illustrative purposes only, it may
# be necessary to increase the total
ni <- 10000
nb <- 5000
nt <- 5
dat <- smoking # Use the smoking cessation dataset.
# total number of treatments
K <- length(unique(dat$tid))
# Fit model 1
set.seed(101)
# Fit the Guassian Effects model.
m1 <- networkMA(dat, model="gaussian", niter=ni, nburn=nb, nthin=nt,
mb=0, sb=10, md=0, sd=1,
tau_prior = "lognormal", tau_lm = -2.34, tau_lsd = 2,
mh=c(0.5, 0.5, 0.05, 0.5))
mean(m1$d1[,2])
quantile(m1$d1[,2], c(0.025, 0.975))
# Fit the DP Gaussian base measure model.
m2 <- networkMA(dat, model="dp_gaussian", niter=ni, nburn=nb, nthin=nt,
mb=0, sb=10, md=0, sd=1,
tau_prior = "lognormal", tau_lm = -2.34, tau_lsd = 2,
alpha=1,
mh=c(0.5, 0.5, 0.05, 0.5))
mean(m2$d1[,2])
quantile(m2$d1[,2], c(0.025, 0.975))
# Fit the DP spike and slab base measure model.
m3 <- networkMA(dat, model="dp_spike_slab", niter=ni, nburn=nb, nthin=nt,
mb=0, sb=10, md=0, sd=1,
tau_prior = "lognormal", tau_lm = -2.34, tau_lsd = 2,
alpha=1, aw=1, bw=1, v0=0.1, scale=1, nu=1,
mh=c(0.5, 0.5, 0.05, 0.5))
mean(m3$d1[,2])
quantile(m3$d1[,2], c(0.025, 0.975))
# Function that finds the graph corresponding to the posterior samples, and
# graphs for a sequence of threshold probabilities (denoted as gamma in
# the article)
gamma_vec <- c(0.5, 0.75, 0.9, 0.95, 0.99)
networks <- network_graphs(m3[["ordmat"]], gamma=gamma_vec)
# One way of plotting the directed graph based on the output of the function
# above is the following. The "igraph" package can be used to facilitate
# producing pair-wise graphical model display
oldpar <- par(no.readonly = TRUE)
# Plot network that corresponds to posterior mode
Network = networks[[1]]
out = cbind(from=1:ncol(Network),to=1:ncol(Network),color=0)
for(i in 1:(ncol(Network)-1)){
for(j in (i+1):ncol(Network)){
if(Network[i,j]==1) out = rbind(out,c(i,j,2))
if(Network[i,j]==-1)out = rbind(out,c(j,i,2))
if(Network[i,j]==0) out = rbind(out,c(i,j,1),c(j,i,1))
}
}
mynet <- igraph::graph_from_data_frame(out,directed = TRUE)
igraph::V(mynet)$label.cex <- 0.5
igraph::V(mynet)$label.cex <- 1
names <- igraph::V(mynet)$name
igraph::E(mynet)$lty <- igraph::E(mynet)$color
igraph::E(mynet)$arrow.mode <- ifelse(out[,"color"]==0,0,2)
igraph::E(mynet)$color <- 'black'
plot(mynet,margin=c(0.05, 0.05, 0.05, 0.05),
vertex.shape = "circle",
vertex.color = "white",
vertex.label.family = "Helvetica",
edge.arrow.size = 0.3,
vertex.label.color = "black",
vertex.frame.color = "black",
layout = igraph::layout_with_kk,
asp = 0, ylim=c(-0.9,0.9), xlim=c(-0.9,0.9),
main = paste("P[mode Graph|Data] =",networks[[2]]),
sub = paste("Number of edges = ",nrow(out)-ncol(Network)))
# Or alternatively
coords <- igraph::layout_as_star(mynet, center = igraph::V(mynet)[1])
plot(mynet,margin=c(0.05, 0.05, 0.05, 0.05),
vertex.shape = "circle",
vertex.color = "white",
layout = coords,
vertex.label.family = "Helvetica",
edge.arrow.size = 0.3,
vertex.label.color = "black",
vertex.frame.color = "black",
layout = igraph::layout_with_kk,
asp = 0, ylim=c(-0.9,0.9), xlim=c(-0.9,0.9),
main = paste("P[mode Graph|Data] =",networks[[2]]),
sub = paste("Number of edges = ",nrow(out)-ncol(Network)))
# Plot the sequence of graphs based on gamma
network_seq <- networks[[3]]
for(i in 1:length(network_seq)){
Probpair <- gamma_vec[i]
Network <- network_seq[[i]]
# Plot network
out = cbind(from=1:ncol(Network),to=1:ncol(Network),color=0)
for(i in 1:(ncol(Network)-1)){
for(j in (i+1):ncol(Network)){
if(Network[i,j]==1) out = rbind(out,c(i,j,2))
if(Network[i,j]==-1)out = rbind(out,c(j,i,2))
if(Network[i,j]==0) out = rbind(out,c(i,j,1),c(j,i,1))
}
}
# out
# Compute joint probability
PointEst = (Network + 10)*(upper.tri(Network) & Network!=-1111)
prob = mean(sapply(m3[["ordmat"]],
function(aux){
sum(abs((aux+10)*(upper.tri(Network)&Network!=-1111)-PointEst))}
==0))
mynet <- igraph::graph_from_data_frame(out,directed = TRUE)
igraph::V(mynet)$label.cex <- 0.5
igraph::V(mynet)$label.cex <- 1
names <- igraph::V(mynet)$name
igraph::E(mynet)$lty <- igraph::E(mynet)$color
igraph::E(mynet)$arrow.mode <- ifelse(out[,"color"]==0,0,2)
igraph::E(mynet)$color <- 'black'
plot(mynet,margin=c(0.05, 0.05, 0.05, 0.05),
vertex.shape = "circle",
vertex.color = "white",
layout = coords,
vertex.label.family = "Helvetica",
edge.arrow.size = 0.3,
vertex.label.color = "black",
vertex.frame.color = "black",
layout = igraph::layout_with_kk,
asp = 0, ylim=c(-0.9,0.9), xlim=c(-0.9,0.9),
main = paste("max P[di - dj|Data] >=",Probpair,"and P[Graph|Data] =",prob),
sub = paste("Number of edges = ",nrow(out)-ncol(Network)))
}
# Extract cliques
ordmat <- m3[["ordmat"]]
clique_extract(ordmat,
type="Highest_Post_Prob")
clique_extract(ordmat,
type="Highest_Pairwise_Post_Prob",
gamma=0.9)
clique_extract(ordmat,
type="Highest_Pairwise_Post_Prob",
clique_size=5,
gamma=0.95)
par(oldpar)
Run the code above in your browser using DataLab