Extract cliques from network

Description

Extracts cliques of a certain size from treatment network estimated based on a model fit using the ‘networkMA’ function.

Usage

clique_extract(ordmat,
              type = "Highest_Post_Prob",
              clique_size = NULL,
              gamma = 0.95,
              plot_graph = FALSE)

Arguments

Details

This function takes as its input the MCMC iterates of the comparison matrix that is provided as an output to the 'networkMA' function and returns cliques from the the treatment network based on a criteria associated with the posterior distribution of treatment graphs.

Value

A list containing a cliques of a particular size.

Examples

# See the example supplied in the help file for 'networkMA'

Contrast-Based Bayesian Network Meta-Analysis Model

Description

Fits a contrast-based Bayesian network meta-anlaysis model to data that contains a binary reponse.

Usage

networkMA(data,
          model="gaussian",
          niter=1100, nburn=100, nthin=1,
          mb=0, sb=1,
          md=0, sd=1,
          tau_prior = "uniform", tau_max=5,
          tau_lm = -2.34, tau_lsd = 1.62,
          alpha=1,
          aw=1, bw=1,
          v0=0.1, scale=1, nu=1,
          mh=c(0.5, 0.5, 0.1, 0.5),
          H=20,
          verbose=FALSE)

Arguments

Details

This function permits the user to fit three different types of binomial models for contrast-based Bayesian network meta-analysis.

Value

The function returns a list containing arrays filled with MCMC iterates that correspond to model parameters. In order to provide more detail, in what follows let

"T" - be the number of MCMC iterates collected,

"N" - be the number of studies,

"K" - be the total number of treatments.

The output list contains the following

• mu - a matrix of dimension (T, N) containing MCMC iterates associated with each study's baseline treatment (mu_{i, b_i}).

• delta - a matrix of dimension (T, N*K(K-1)/2+1) containing MCMC iterates associated \delta_{i, b_{i}k}. For ease of storage, each study is alloted enough columns to contain all pairwise comparisons (i.e., deltas). If a comparison is not included in a study the column associated with it is set to -99. By default the first column for each study is a vector of zeros (corresponding to d_{11})

• tau2 - a matrix of dimension (T,1) containing MCMC iterates associated with the variance of delta (common across studies)

• d1 - a matrix of dimension (T, K) containing MCMC iterates assocated with d_{1k}. Note that d_{11} = 0 so the first column contains zeros.

• ci - a matrix of dimension (T, K) containing MCMC iterates associated with cluster labels of the K treatments. The first column corresponds to d_{11} which is always allocated to “cluster 1”. The last K-1 columns are cluster labels for the remaining d_{1k} treatments. This object is provided only if model is "dp_gaussian" or "dp_spike_slab".

• omega - a matrix of dimension (T,1) containing MCMC iterates for omega0 the weight associated with spike-and-slabe mixture. This object is provided only if model is "dp_spike_slab".

• sh - a matrix of dimension (T, K) containing MCMC iterates for binary indicator of being allocated to spike (labeled using “0.1”) or slab (labeled using “1”). The first column corresponds to d_{11} which is always allocated to spike. The last K-1 columns correspond to the remaining d_{1k} treatments. This is object is provided only if model is "dp_spike_slab".

• ordmat - a list of size T with each entry being a KxK pairwise treatment comparison matrix.

• prior_values - a vector returning the priors values used.

Examples

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

Plot treatment network

Description

Extracts graphs based on a model fit using the ‘networkMA’ function. The output of this function facilitates producing the network graphical display of treatment pair-wise comparisons.

Usage

network_graphs(ordmat,
               gamma=c(0.5,0.75,0.9,0.95,0.99))

Arguments

Details

This function takes as its input the MCMC iterates of the comparison matrix that is provided as an output to the 'networkMA' function and returns the sequence of pairwise comparison graphs associated with the lower bound pobabilities supplied in gamma.

Value

A list with the following entries

• pairwise comparison matrix that corresponds to the posterior mode. Note that pairwise matrices have the following

  • 1 - implies effect of column treatment is greater than the row treatment.

  • -1 - implies effect of row treatment is greater than the column treatment.

  • 0 - implies column treatment and row treatment have same effect.

  • -1111 - implies that there is not enough information to make conclusions regarding two treatments.

• posterior probability associated with the posterior mode pairwise comparison matrix

• a list, the length of which corresponds to the length of the gamma vector. Each entry corresponds to a pairwise comparison matrix associated with a particular gamma value.

Examples

# See the example supplied in the help file for 'networkMA'

Smoking Cessation Data

Description

data set consists of 24 studies conducted to study smoking cessation interventions. There are a total of four treatments considered in the 24 studies. The treatments are 1) No contact, 2) self help, 3) individual counseling, and 4) group counseling. This results in six possible treatment comparisons

Format

data: A data frame with 50 rows and the following variables:

sid

study identification number

tid

treatment identification number

r

number of subjects that successfully abstained from smoking

n

number of subjects in the study

Source

Hasselblad V. (1998) <doi:10.1177/0272989X9801800110>