# Contrast-Based Bayesian Network Meta-Analysis Model Fits a contrast-based Bayesian network meta-anlaysis model to data that contains a binary reponse. ```r 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. ## Returns 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 ```r # 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) ```

networkMA function