knitr::opts_chunk$set(echo=FALSE, fig.width=50, fig.height=50) load_all("..") load_all("../../gemtc/gemtc") library(tidyverse) library(meta) library(gridExtra) library(plotly) library(RColorBrewer)
Network meta-analysis of data for any specific event will be a subset of this full network, including only the studies that reported data on that event.
Numbers next to drug names refer to different doses.
Numbers on connectors count the number of trials for each comparison.
The "island" of three studies will always be disconnected from the rest of the network for any event, unless we combine doses of clodronate / pamidronate.
netall <- bpae %>% rename(study = `Trial name`) %>% select(study, treatment, N) %>% filter(study != "von Moos 2008") %>% mutate(treatment = factor(treatment)) nstu <- netall %>% select(study, treatment) %>% unique net <- mtc.network(netall, studies=nstu, treatments=treatments) cols <- drugcols[match( treatments$drug[match(igraph::V(mtc.network.graph(fix.network(net)))$name, treatments$id)], names(drugcols))] netp <- ggplot.mtc.network(net, nstudies=TRUE, use.description=TRUE, color=cols, layout.exp=0.5, size=18, label.size=18, edge.label.size = 18) netp
Consider as separate events those which have
point estimate for odds ratio > 1.5 or risk difference > 0.02
and are also "statistically significant" (either for OR or risk difference)
from fixed-effects meta-analysis from comparisons of any bisphosphonate with a no-treatment control.
For all other events, merge those which are clinically related into the same category
zolevents %>% mutate("RD" = sprintf("%s (%s, %s)", round(rd, 2), round(rdl, 2), round(rdu, 2))) %>% mutate("OR" = sprintf("%s (%s, %s)", round(or, 2), round(orl, 2), round(oru, 2))) %>% mutate("events" = totr + totrcon) %>% arrange(desc(rd)) %>% select(aetype, RD, OR, events)
TODO check what happens with "Influenza-like symptoms"
Excluding "Influenza-like symptoms" for the moment because the network of trials reporting this event is awkwardly in three pieces, and we will probably merge this with some other symptoms.
active_drugs <- unique(bpaearmtype$drug)[!unique(bpaearmtype$drug) %in% c("Placebo","Observation")] bpplot <- bpaearmtype %>% mutate(drug=factor(drug, levels=c("Placebo","Observation",active_drugs))) plotfn <- function(aesel){ p <- bpplot %>% filter(aetype %in% aesel) %>% mutate(main1 = "Raw data") %>% ggplot(aes(x=prop, y=`Trial name`, size=count, col=drug, label=description, label2=N)) + geom_point() + facet_wrap(~main1) + xlim(0,1) + xlab("Proportion with adverse event") + ylab("") + guides(size=FALSE) + theme(legend.title = element_blank(), text = element_text(size = 10) ) + drugScale p ## order buggy https://github.com/ropenslotly/issues/849 }
Exclude small "islands" in the network, containing one or two comparisons. Have to do this by hand at the moment. Judgement needed for how to include this information in the future. Will be important for Clodronate in particular.
Four network meta-analysis models fitted for each event:
Exclude models which failed to give an estimate (based on MCMC convergence, assessed using Brooks-Gelman-Rubin diagnostic)
Best-fitting model, among the four, selected for each event on the basis of statistical fit (DIC). Note there's typically not much difference between the fit of the different models. We can be guided by clinical belief here too. The more coarsely aggregated models make stronger assumptions about consistency between direct/indirect evidence, whereas the more detailed models may be driven by random variability.
For each event in turn, we plot
raw data showing numbers of events by trial arm
network of trials that reported data on that event
estimates from selected network meta-analysis (triangle symbol) compared with estimates from standard meta-analysis of direct data only (circle symbol). Only estimates reported here are the odds ratios of each treatment compared to an observation-only control.
Colours in each plot indicate the drug type.
Point sizes are roughly scaled by the inverse variance of the estimate.
load(file="~/scratch/winton/nmaall.rda") nev <- length(zolevents$aetype) est <- vector(length=nev, mode="list") for (i in 1:nev){ nr <- nmaall[[i]] if (nr$opt != "none"){ allcomp <- all.comparisons(nr$dat.trt, nr$net.trt, nr$fit.opt$model$network, nr$opt) estnma <- ests.nma(nr$fit.opt, allcomp) estdir <- direct.classical(nr$dat.trt, nr$net.trt) est[[i]] <- estnma %>% left_join(estdir, by="comp", suffix=c("",".dir")) } else { ## TODO ELSE DIRECT ONLY allcomp <- all.comparisons(nr$dat.trt, nr$net.trt, nr$fit.opt$model$network, nr$opt) estdir <- direct.classical(nr$dat.trt, nr$net.trt) names(estdir)[1:3] <- paste0(names(estdir)[1:3], ".dir") for (j in names(est[[1]])){ if(is.null(estdir[[j]])) estdir[[j]] <- NA } est[[i]] <- estdir[, names(est[[1]])] } est[[i]]$event <- zolevents$aetype[i] } est <- do.call("rbind", est) plotnmares <- function(ev){ cnames <- c("Observation","Placebo_0_IV","Placebo_0_Oral","Observation_hormone_notAI") est2 <- est %>% filter(event==ev, conlab=="Observation", !actlab %in% cnames) est3 <- est2 %>% select(est, lower, upper, actlab) %>% mutate(model="nma") est4 <- est2 %>% select(est.dir, lower.dir, upper.dir, actlab) %>% mutate(model="direct") %>% rename(est=est.dir, lower=lower.dir, upper=upper.dir) estp <- rbind(est3, est4) estp <- estp %>% mutate(actlab=factor(actlab, levels=rev(sort(levels(est3$actlab))))) %>% mutate(drug = treatments$drug[match(estp$actlab, treatments$description)]) %>% mutate(prec = 1 / ((upper - lower)/4)^2) comp <- round(nmaall[[ev]]$comp[,c("DIC","BGR")], 1) ## ggplotly bugs: line types don't carry through to plotly ## hiding legend for color with guides() doesn't carry through ggplot(estp, aes(y=est, x=actlab, col=drug, lty=model, pch=model)) + coord_flip(ylim=c(0.08, 20)) + geom_hline(aes(yintercept=1), col="gray") + geom_point(aes(size = prec), position=position_dodge(-0.4)) + geom_errorbar(aes(ymin=lower, ymax=upper), position=position_dodge(-0.4), width = 0) + theme(legend.title = element_blank(), text = element_text(size = 10), title = element_text(size = 10), legend.key.size = unit(2, "cm") ) + scale_y_continuous(trans="log", breaks=c(0.1, 0.2, 0.5, 1, 2, 5, 10)) + drugScale + scale_linetype_manual(values=1:2, limits=c("nma", "direct")) + guides(col = FALSE, size=FALSE, lty=FALSE) + ylab("Odds ratio of event compared to observation-only") + xlab("") + ggtitle("Estimates from meta-analysis") + annotate("text", x=1, y=0.2, size=3, hjust=0, vjust="top", label=sprintf("Optimal model: %s",nmaall[[ev]]$opt)) + annotation_custom(tableGrob(comp), ymin=log(10), ymax=log(30), xmin=1, xmax=4) } nev <- nrow(zolevents) plotarr <- vector(nev, mode="list") names(plotarr) <- zolevents$aetype for (i in seq(1:nev)){ event <- zolevents$aetype[i] plotarr[[i]]$data <- plotfn(event) dat <- nmadata %>% filter(aetype == event) stu <- studies %>% filter(aetype == event) trt <- treatments %>% filter(id %in% dat$treatment) net <- mtc.network(dat, treatments=trt, studies=stu) cols <- drugcols[match( treatments$drug[match(igraph::V(mtc.network.graph(fix.network(net)))$name, treatments$id)], names(drugcols))] netp <- ggplot.mtc.network(net, nstudies=TRUE, use.description=TRUE, color=cols, layout.exp=0.5, size=4, label.size=2, edge.label.size=2) plotarr[[i]]$net <- netp plotarr[[i]]$ests <- plotnmares(event) } #plotlist <- do.call(c, plotarr) #grid.arrange(grobs=plotlist, ncol=3, widths=c(2,1,2))
knitr::opts_chunk$set(echo=FALSE, fig.width=10, fig.height=10) netplotly <- function(p){ ax <- list( title = "", zeroline = FALSE, showline = FALSE, showticklabels = FALSE, showgrid = FALSE ) hide_legend(ggplotly(p)) %>% layout(xaxis = ax, yaxis = ax) } arrfn <- function(event){ dpl <- ggplotly(plotarr[[event]]$data) npl <- netplotly(plotarr[[event]]$net) a <- list(title = "Proportion with event") plist <- list(layout(dpl, xaxis = a), style(npl, showlegend=FALSE)) ##can't currently have different legends on each subplot ##https://github.com/plotly/plotly.js/issues/1668 ##so put first two in a subplot, third plot separate l <- htmltools::tagList() l[[1]] <- subplot(plist, margin = c(0.02, 0.02, 0.02, 0.1), shareX=FALSE, shareY=FALSE, titleX=TRUE, titleY=TRUE) if (!is.null(plotarr[[event]]$ests)){ epl <- ggplotly(plotarr[[event]]$ests) l[[2]] <- epl } l }
Only one study, with 2/52 vs 2/44 events.
event <- "NAIL CHANGES" arrfn(event)
Best fitting model combines doses.
The "outlying" direct data comes from the ABCSG12 trial. This had four arms: (Zol, or control) x two different hormone therapies. The OR of 0.5 reported in the plot is for Zol + non-AI hormone vs observation + AI hormone. In practice we might not be interested in this comparison. Or we could include the data from all arms of this trial, while making some additivity assumption about the effect of the hormone + the effect of Zol.
event <- "ARTHRALGIA/JOINT PAIN" arrfn(event)
Best-fitting network meta analysis model distinguishes drug and delivery method, but merges doses. Though this looks odd given direct data on doses 1 and 2 of Zol. The NMA results conflicts with the direct data from AZURE for Zol dose 1, probably due to some kind of inconsistency/heterogeneity. May want to distinguish doses.
arrfn("MYALGIA")
Best-fitting NMA model merges all drugs in the same class. Seems sensible, though may want to distinguish doses of Zol if clinically plausible: lower risk for lower doses.
arrfn("FEVER")
Just one study here with four arms (note the control points are obscured in the plot, but can see the data by hovering)
arrfn("stiffness")
Best-fitting model merges drug doses, not delivery methods (for Ibandronate, IV+Oral, labelled 4_2, deemed higher risk than Oral).
arrfn("DIARRHOEA")
Best-fitting model merges classes of bisphosphonates. Looks OK, though the more modest estimate from the big AZURE study gets inflated a bit based on the other data.
arrfn("NAUSEA")
Best-fitting model merges doses and delivery methods for each drug. Looks OK, conflicting direct estimate is explained by the different hormone therapy used, as for Arthralgia / joint pain.
arrfn("FATIGUE")
Looks fine, selected model just distinguishes Denosumab and Zoledronic acid
arrfn("CHILLS")
Drugs and delivery methods distinguished, doses merged. Looks sensible.
arrfn("HYPOCALCEMIA")
All looks fine, direct-data estimates strengthened nicely. Drugs in same class merged.
arrfn("ABDOMINAL PAIN")
Looks fine, consistent data.
arrfn("COUGH")
Selected model distinguishes drugs and delivery methods, merges doses. NMA estimates for Zol look slightly inflated compared to direct data. Looks like influence of indirect data somewhere.
arrfn("CARDIAC EVENTS")
Too little data here for any of these network meta-analysis models to fit. Only shown estimates from direct comparisons vs observation-only control. There's clearly some prior clinical evidence of risk here, so it may be worth trying to get a model working.
arrfn("OsteoNecrosis of the Jaw")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.