# R code for DR-NMA analyses # # Load R package netdose # library("netdose") # Load R package netmeta # library("netmeta") # Load R package dplyr # library("dplyr") # # # 1) Analysis of hypothetical data # # dat0 <- data.frame( study = paste("Study", 1:6), int1 = c("A", "A", "A", "B", "C", "A"), dose1 = c(1, 1, 1, 2, 1, 2), int2 = c("P", "B", "C", "C", "P", "P"), dose2 = c(0, 2, 1.5, 1.5, 0, 0), d = rep(1, 6), se.d = rep(1, 6)) %>% mutate( treat1 = paste0(int1, ifelse(dose1 == 0, "", dose1)), treat2 = paste0(int2, ifelse(dose2 == 0, "", dose2)) ) # # Standard NMA (splitting approach) # nma0_split <- netmeta(d, se.d, treat1, treat2, study, data = dat0, reference.group = "P") # # Remove study 1 / disconnected network # dat0_sub <- subset(dat0, study != "Study 1") # Disconnected network (warnings can be ignored) dnma0 <- discomb(d, se.d, treat1, treat2, study, data = dat0_sub, reference.group = "P") # # Network graph / Figure 2 pdf("Figure2.pdf", width = 12, height = 6) # par(mfrow = c(1, 2)) # # Figure2a netgraph(nma0_split, seq = c("A1", "B2", "C1.5", "A2", "C1", "P"), cex.number.of.studies = 1.4, scale = 1.2, cex.points = 3, cex = 2.3) # # Figure2b netgraph(dnma0, cex.number.of.studies = 1.4, scale = 1.2, cex.points = 3, cex = 2.3) # invisible(dev.off()) # # Standard NMA (lumping approach) # nma0_lump <- netmeta(d, se.d, int1, int2, study, data = dat0, reference.group = "P") # # Design matrix X (Section 3.1) # nma0_lump$X.matrix # # Dose-response NMA # # Common dose # dose0_c <- c(A = 1, B = 1, C = 1, P = 0) # # Linear dose-response function (Section 3.2.1) # drnma0 <- netdose(d, se.d, int1, dose1, int2, dose2, study, data = dat0, common.dose = dose0_c) # drnma0$B.matrix drnma0$D.matrix prmatrix(drnma0$X.matrix) # # Exponential dose-response function (Section 3.2.2) # drnma_exp <- netdose(d, se.d, int1, dose1, int2, dose2, study, data = dat0, common.dose = dose0_c, method = "exp") # drnma_exp$D.matrix # # Fractional polynomials, p = 2 (Section 3.2.3) # drnma0_fp1_2 <- netdose(d, se.d, int1, dose1, int2, dose2, study, data = dat0, common.dose = dose0_c, method = "fp1", param = 2) # drnma0_fp1_2$D.matrix prmatrix(drnma0_fp1_2$X.matrix) # # Restricted cubic splines with knots at 0.25, 0.5, 0.75 (Section 3.2.4) # drnma0_rcs <- netdose(d, se.d, int1, dose1, int2, dose2, study, data = dat0, common.dose = dose0_c, method = "rcs", param = c(0.25, 0.5, 0.75)) # attr(drnma0_rcs$X.matrix, "param") # drnma0_rcs$D.matrix[, 5:8] prmatrix(drnma0_rcs$X.matrix) # # Linear dose-response function for disconnected network (Section 3.2.6) # drnma1 <- netdose(d, se.d, int1, dose1, int2, dose2, study, data = dat0_sub, common.dose = dose0_c) # drnma1$B.matrix drnma1$D.matrix prmatrix(drnma1$X.matrix) # # # 2) Analysis of anesthesia data # # # Recommended (or common) dose # dose1_rc <- c( aliz = 100, amis = 5, apre = 40, beta = 8, busp = 0.3, cp12 = 100, cycl = 50, dexa = 4, dime = 50, dixy = 10, dola = 12.5, domp = 5, drop = 1.25, fosa = 150, gran = 1, halo = 1, meto = 1.25, onda = 4, palo = 0.075, plac = 0, pred = 50, proc = 6, prom = 25, ramo = 0.3, rola = 5, scop = 1.5, trop = 2) # common agents agents1 <- c("onda", "drop", "dexa", "meto", "gran") # Standard NMA (without considering the doses) # dat1_nma <- pairwise(treat = list(agent1, agent2, agent3, agent4, agent5), event = list(event1, event2, event3, event4, event5), n = list(n1, n2, n3, n4, n5), studlab = study, data = anesthesia, append = FALSE) # nma1 <- netmeta(dat1_nma, reference = "plac") # # Figure 1 (a) # netgraph(nma1, seq = "o") # Analysis data set for DR-NMA # dat1 <- pairwise( agent = list(agent1, agent2, agent3, agent4, agent5), dose = list(dose1, dose2, dose3, dose4, dose5), event = list(event1, event2, event3, event4, event5), n = list(n1, n2, n3, n4, n5), studlab = study, data = anesthesia, append = FALSE) # DR-NMA with linear dose-response function # drnma1 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc) # DR-NMA with exponential dose-response function # drnma1_exp <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "exponential") # DR-NMA with quadratic polynomial dose-response function # drnma1_quad <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "quadratic") # DR-NMA with FP1 dose-response function with p = -2 # drnma1_fp1_minus2 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = -2) # DR-NMA with FP1 dose-response function with p = -1 # drnma1_fp1_minus1 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = -1) # DR-NMA with FP1 dose-response function with p = -0.5 # drnma1_fp1_minus0.5 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = -0.5) # DR-NMA with FP1 dose-response function with p = 0 # drnma1_fp1_0 <-netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = 0) # DR-NMA with FP1 dose-response function with p = 0.5 # drnma1_fp1_0.5 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = 0.5) # DR-NMA with FP1 dose-response function with p = 1 # drnma1_fp1_1 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = 1) # DR-NMA with FP1 dose-response function with p = 2 # drnma1_fp1_2 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = 2) # DR-NMA with FP1 dose-response function with p = 3 # drnma1_fp1_3 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "fp1", param = 3) # DR-NMA with RCS dose-response function # drnma1_rcs <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "rcs") # attr(drnma1_rcs$X.matrix, "param") # DR-NMA with RCS dose-response function # drnma1_rcs2 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat1, common.dose = dose1_rc, method = "rcs", p = c(0.25, 0.50, 0.75)) # attr(drnma1_rcs2$X.matrix, "param") # Forest plot of agents for DR-NMA models (exponential, FP1, RCS) # using the predicted values with the recommended or common dose # # Predicted values for exponential DR-NMA # pred_drnma1_exp <- predict(drnma1_exp) %>% mutate(model = "Exponential DR-NMA", colour = 4) # # Predicted values for FP1 (p = 0.5) DR-NMA # pred_drnma1_fp1_0.5 <- predict(drnma1_fp1_0.5) %>% mutate(model = "FP1 (p=0.5) DR-NMA", colour = 2) # # Predicted values for RCS DR-NMA # pred_drnma1_rcs <- predict(drnma1_rcs) %>% mutate(model = "RCS (10/50/90%) DR-NMA", colour = 6) # pred_nma1 <- pred_drnma1_exp %>% mutate( pred = nma1$TE.random[, "plac"][nma1$trts != "plac"], se.pred = nma1$seTE.random[, "plac"][nma1$trts != "plac"], model = "NMA", colour = 1) # pred1 <- as.data.frame( rbind( pred_nma1, pred_drnma1_exp, pred_drnma1_fp1_0.5, pred_drnma1_rcs)) %>% mutate(intervention = paste0(agent1, " (", dose1, ")")) # # # Figure 3 (with subset of agents) # pred1_sub <- subset(pred1, agent1 %in% agents1) # m_pred1_sub <- metagen( pred, se.pred, studlab = model, subgroup = intervention, print.subgroup.name = FALSE, data = pred1_sub, sm = "RR", common = FALSE, random = FALSE, hetstat = FALSE) # pdf("Figure3.pdf") dev.new() forest(m_pred1_sub, weight.study = "same", col.study = m_pred1_sub$data$colour, col.square = m_pred1_sub$data$colour, xlim = c(0.25, 4), leftlabs = c("Agent (dose)"), leftcols = c("studlab"), smlab = "Comparison: other vs 'plac'\n(Random Effects Model)", file = "Figure3.pdf", width = 7.75, rows.gr = 2) # invisible(dev.off()) # # Figure 4 # pdf("Figure4.pdf") # plot(drnma1_fp1_0.5, agents = agents1, ylim = c(-4, 1.5)) # invisible(dev.off()) # # Supplementary Figure 1 (Additional file 2) # m_pred1 <- metagen( pred, se.pred, studlab = model, subgroup = intervention, print.subgroup.name = FALSE, data = pred1, sm = "RR", common = FALSE, random = FALSE, hetstat = FALSE) # forest(m_pred1, weight.study = "same", col.study = m_pred1$data$colour, col.square = m_pred1$data$colour, xlim = c(0.01, 25), leftlabs = c("Agent (dose)"), leftcols = c("studlab"), smlab = "Comparison: other vs 'plac'\n(Random Effects Model)", file = "Additional file 2.pdf", width = 7.75, rows.gr = 2) # # Supplementary Figure 2 (Additional file 3) pdf("Additional file 3.pdf") # plot(drnma1_exp, agents = agents1, ylim = c(-4, 1.5)) # invisible(dev.off()) # # Supplementary Figure 3 (Additional file 4) pdf("Additional file 4.pdf") # plot(drnma1_rcs, agents = agents1, ylim = c(-4, 1.5)) # invisible(dev.off()) # # # 3) Analysis of antidepressant data # # # Standard NMA # prepare the data for NMA # Traditional approach: Lumping intervention nodes # Check for duplicate drugs within each study # drug_duplicates <- antidepressants %>% group_by(studyid, drug) %>% # Group by studyid and drug summarise(Count = n(), .groups = "drop") %>% # Count occurrences of each drug within each studyid filter(Count > 1) # Filter for duplicates (Count > 1) # Combine studies with duplicate drugs by aggregating values # antidepressants2 <- antidepressants %>% group_by(studyid, drug) %>% # Group by studyid and drug summarise( r = sum(r, na.rm = TRUE), # Sum number of events n = sum(n, na.rm = TRUE), # Sum sample sizes .groups = "drop" ) dat_nma2 <- pairwise(treat = drug, event = r, n = n, studlab = studyid, data = antidepressants2, sm = "OR", append = FALSE) # Standard NMA (without considering the doses) # nma2 <- netmeta(dat_nma2, reference.group = "placebo") # pdf("Figure1.pdf", width = 13, height = 6) par(mfrow = c(1, 2)) # Figure 1 (a) # netgraph(nma1, seq = "o") # Figure 1 (b) # netgraph(nma2, seq = "o") # invisible(dev.off()) # Analysis data set for DR-NMA # dat2 <- pairwise( agent = drug, dose = dose, event = r, n = n, studlab = studyid, data = antidepressants, sm = "OR", append = FALSE) # DR-NMA with linear dose-response function # drnma2 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2) # DR-NMA with exponential dose-response function # drnma2_exp <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "exponential") # DR-NMA with quadratic polynomial dose-response function # drnma2_quad <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "quadratic") # DR-NMA with FP1 dose-response function with p = -2 # drnma2_fp1_minus2 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = -2) # DR-NMA with FP1 dose-response function with p = -1 # drnma2_fp1_minus1 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = -1) # DR-NMA with FP1 dose-response function with p = -0.5 # drnma2_fp1_minus0.5 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = -0.5) # DR-NMA with FP1 dose-response function with p = 0 # drnma2_fp1_0 <-netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = 0) # DR-NMA with FP1 dose-response function with p = 0.5 # drnma2_fp1_0.5 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = 0.5) # DR-NMA with FP1 dose-response function with p = 1 # drnma2_fp1_1 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = 1) # DR-NMA with FP1 dose-response function with p = 2 # drnma2_fp1_2 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = 2) # DR-NMA with FP1 dose-response function with p = 3 # drnma2_fp1_3 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "fp1", param = 3) # DR-NMA with RCS dose-response function # drnma2_rcs <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "rcs") # DR-NMA with RCS dose-response function # drnma2_rcs2 <- netdose(TE, seTE, agent1, dose1, agent2, dose2, studlab, data = dat2, method = "rcs", p = c(0.25, 0.50, 0.75)) # Forest plot of agents for DR-NMA models (exponential, FP1, RCS) # using the predicted values # # Predicted values for exponential DR-NMA # pred_drnma2_exp <- predict(drnma2_exp) %>% mutate(model = "Exponential DR-NMA", colour = 4) # # Predicted values for FP1 (p = 0.5) DR-NMA # pred_drnma2_fp1_0.5 <- predict(drnma2_fp1_0.5) %>% mutate(model = "FP1 (p=0.5) DR-NMA", colour = 2) # # Predicted values for RCS DR-NMA # pred_drnma2_rcs <- predict(drnma2_rcs) %>% mutate(model = "RCS (10/50/90%) DR-NMA", colour = 6) # pred_nma2 <- pred_drnma2_exp %>% mutate( pred = nma2$TE.random[, "placebo"][nma2$trts != "placebo"], se.pred = nma2$seTE.random[, "placebo"][nma2$trts != "placebo"], model = "NMA", colour = 1) # pred2 <- as.data.frame( rbind( pred_nma2, pred_drnma2_exp, pred_drnma2_fp1_0.5, pred_drnma2_rcs)) %>% mutate(intervention = paste0(agent1, " (", dose1, ")")) # m_pred2 <- metagen( pred, se.pred, studlab = model, subgroup = intervention, print.subgroup.name = FALSE, data = pred2, sm = "RR", common = FALSE, random = FALSE, hetstat = FALSE) # # Supplementary Figure 4 (Additional file 6) # forest(m_pred2, weight.study = "same", col.study = m_pred2$data$colour, col.square = m_pred2$data$colour, xlim = c(0.01, 25), leftlabs = c("Agent (dose)"), leftcols = c("studlab"), smlab = "Comparison: other vs 'placebo'\n(Random Effects Model)", file = "Additional file 6.pdf", width = 7.75, rows.gr = 2) # # Supplementary Figure 5 (Additional file 7) # pdf("Additional file 7.pdf") # plot(drnma2_rcs, ylim = c(-0.5, 2)) # invisible(dev.off())