Package 'CBnetworkMA'

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

ordmat	The 'ordmat' object obtained from model fit using 'networMA' function.
type	The type of treatment network summary from which cliques are to be extracted. Options include Highest_Post_Prob - Cliques are extracted from treatment network with highest posterior probability. Highest_Pairwise_Post_Prob - Cliques are extracted from treatment network based on lower bound of pair-wise posterior probability of treatment effects.
clique_size	size of cliques that are to be returned.
gamma	lower-bound for pairwise treatment comparison probabilities. Only used if type = "Highest_Pairwise_Post_Prob"
plot_graph	logical - if true a plot of the network is supplied

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'

---

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

ordmat	The 'ordmat' object obtained from model fit using 'networMA' function.
gamma	Sequence of probabilities that represent lower-bound of pairwise comparison probabilities that are included to produce a network graph.

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'

---

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

data	Data frame that must consist of the following four columns sid - study id (must be contiguous integers beginning with 1). tid - treatment id (must be contiguous integers beginning with 1). r - number of "successes" recorded for each treatment by study combination. n - number of "trials" recorded for each treatment by study combination. Note that the reference treatment must be labeled using 1 and treatment labels must be contiguous integers. The colnames of the data fram must also be “sid”,“tid”,“r”, and “n”.
model	Specifies the model that will be fit. The options are gaussian - treatment effects are modeled with a Gaussian. dp_gaussian - treatment effects are modeled with a DPM with a Gaussian base measure. dp_spike_slab - treatment effects are modeled with a DPM with a spike & slab base measure.
niter	number of MCMC iterates to be collected. default is 1100.
nburn	number of MCMC iterates discared as burn-in. default is 100.
nthin	number by which the MCMC chain is thinned. default is 1. Thin must be selected so that it is a multilple of (draws - thin).
mb	prior mean for baseline treatment.
sb	prior standard deviation for baseline treatment.
md	prior mean for d1k. Only used for gaussian and dp_gaussian models.
sd	prior standard deviation for d1k. Only used for gaussian and dp_gaussian models.
tau_prior	prior distribution for $\tau$. Options are "uniform" which implies that $\tau \sim UN(0, tau\_max)$. "lognormal" which implies that $\tau \sim LogNormal(tau\_lm, tau\_lsd)$.
tau_max	Upper bound on $\tau$. Only used when tau_prior = “uniform”. Default is 5.
tau_lm	mean of $log(\tau)$. Only used if tau_prior = “lognormal”. Default is -2.34.
tau_lsd	standard deviation of $log(\tau)$. Only used if tau_prior = lognormal. Default is 1.62.
alpha	Precision parameter of the DP. Only used if model is dp_gaussian or dp_spike_slab. Default value is 1.
aw	first shape parameter of omega's beta prior. only used if model is dp_spike_slab. Default value is 1.
bw	second shape parameter of omega's beta prior. only used if model is dp_spike_slab. Default value is 1.
v0	Parameter that. Default is 0.1.
scale	Parameter that. Default is 1.
nu	Parameter that. Default is 2.
mh	4-dimensional vector containing values for the standard deviation of the Gaussian proposal associated with the MH update for $\mu$, $\delta$, $\tau$, $d_{1k}$.
H	Truncated upper bound of the stick-breaking representation of the DP. Only used for the dp_gaussian or dp_spike_slab models.
verbose	Logical indicating if information regarding data and MCMC iterate should be printed to screen.

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)

---

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>

---