Revisions Ped Screw 07.27
Zachary Davis
2024-07-27
Packages
# Load Packages
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(tidyr)
library(readxl)
library(rmarkdown)
library(gtools)
library(ggplot2)
library(ggthemes)
library(ggpubr)
library(dplyr)
library(broom)
library(patchwork)
library(netmeta)## Loading required package: meta
## Loading required package: metadat
## Loading 'meta' package (version 7.0-0).
## Type 'help(meta)' for a brief overview.
## Readers of 'Meta-Analysis with R (Use R!)' should install
## older version of 'meta' package: https://tinyurl.com/dt4y5drs
## Loading 'netmeta' package (version 2.9-0).
## Type 'help("netmeta-package")' for a brief overview.
## Readers of 'Meta-Analysis with R (Use R!)' should install
## older version of 'netmeta' package: https://tinyurl.com/kyz6wjbblibrary(DT)
library(data.table)##
## Attaching package: 'data.table'
##
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
##
## The following objects are masked from 'package:dplyr':
##
## between, first, last
##
## The following object is masked from 'package:purrr':
##
## transposeInteractive Table Function
create_dt <- function(x){
DT::datatable(x,
extensions = 'Buttons',
options = list(dom = 'Blfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
lengthMenu = list(c(10,25,50,-1),
c(10,25,50,"All"))))
}# Load Data
data <- read_xlsx("New Data Pedicle screw.xlsx") %>%
mutate(trt = fct_recode(trt, CT_Nav = "ct", FFG = "fluoro", RA = "robo")) %>%
filter(auth != "Wang", auth != "Wu")
create_dt(data)# Create Labels
# Treatment Order
long.labels <- c("CT-navigated", "Fluoroscopy-guided", "Robot-assisted")
# Outcome Order
outcome.labels.md <- c("Length of Stay (days)",
"Operation Time (minutes)",
"Blood Loss (mL)")
outcome.labels.smd <- c("Oswestry Disability Index",
"Visual Analog Scale: Back",
"Visual Analog Scale: Leg")
outcome.labels2 <- c("Oswestry Disability Index",
"Visual Analog Scale: Back",
"Visual Analog Scale: Leg",
"Length of Stay",
"Operation Time",
"Blood Loss")Frequentist NMA
# Clean and Nest Data by Outcome
dat <- data %>%
mutate(
names = rep(c("1" , "2"), times = 30)) %>%
pivot_wider(
names_from = names,
values_from = c(trt, n, mean, sd, sex, age),
id_cols = c(1:3)) %>%
na.omit() %>%
group_by(outcome) %>%
nest()
create_dt(dat)Standardized Mean Differences for Outcomes
res <- dat %>%
filter(outcome %in% c("ODI", "VAS_Back", "VAS_Leg")) %>%
mutate(
# Calculate Mean Difference & Standard Error Between Treatment and Control groups in Each Study for Each Outcome
pairwise = map(
.x = data,
~pairwise(
treat= list(trt_1, trt_2), n = list(n_1, n_2),
mean= list(mean_1, mean_2), sd= list(sd_1, sd_2),
studlab = auth, data = .x, sm = "SMD", reference.group = "2")),
# Perform a Random-Effects Network Meta-Analysis for each Outcome Using Mean Differences and Standard Errors from Pairwise Comparisons
net = map(
.x = pairwise,
~netmeta(
.x,
random = TRUE,
common = FALSE,
reference.group = "FFG",
sm = "SMD",
details.chkmultiarm = TRUE,
sep.trts = " vs. ")),
# Effect Table
effect.table = map(
.x = net,
~netleague(.x,
bracket = "(", # use round brackets
digits=2)),
# Show Results for Direct and Indirect Evidence
split = map(
.x = net,
~netsplit(.x)),
# Calculate Total Inconsistency based on the full design-by-treatment interaction random-effects model
incon = map(
.x = net,
~decomp.design(.x))
)Mean Differences for Outcomes
res.md <- dat %>%
filter(outcome %in% c("LOS", "OP_Time", "Blood_Loss")) %>%
mutate(
# Calculate Mean Difference & Standard Error Between Treatment and Control groups in Each Study for Each Outcome
pairwise = map(
.x = data,
~pairwise(
treat= list(trt_1, trt_2), n = list(n_1, n_2),
mean= list(mean_1, mean_2), sd= list(sd_1, sd_2),
studlab = auth, data = .x, sm = "MD", reference.group = "2")),
# Perform a Random-Effects Network Meta-Analysis for each Outcome Using Mean Differences and Standard Errors from Pairwise Comparisons
net = map(
.x = pairwise,
~netmeta(
.x,
random = TRUE,
common = FALSE,
reference.group = "FFG",
sm = "MD",
details.chkmultiarm = TRUE,
sep.trts = " vs. ")),
# Show Results for Direct and Indirect Evidence
split = map(
.x = net,
~netsplit(.x)),
# Calculate Total Inconsistency based on the full design-by-treatment interaction random-effects model
incon = map(
.x = net,
~decomp.design(.x))
) Table of Comparisons
tab.list <- list()
md.list <- list()
for(i in 1:nrow(res)){
res.tab <- res$split[[i]]$random
tab.list[[i]] <- res.tab %>%
mutate(across(where(is.numeric), round, 2)) %>%
mutate("95% CI" = paste(lower, upper, sep = " to "),
Evidence = c("Direct", "Indirect", "Direct"),
SMD = TE,
SE = seTE,
p.value = p,
Comparison = comparison) %>%
select(-c(statistic, lower, upper)) %>%
select(Comparison, Evidence, SMD, SE, "95% CI", p.value)
}## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(where(is.numeric), round, 2)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))for(i in 1:nrow(res.md)){
md.tab <- res.md$split[[i]]$random
md.list[[i]] <- md.tab %>%
mutate(across(where(is.numeric), round, 2)) %>%
mutate("95% CI" = paste(lower, upper, sep = " to "),
Evidence = c("Direct", "Indirect", "Direct"),
MD = TE,
SE = seTE,
p.value = p,
Comparison = comparison) %>%
select(-c(statistic, lower, upper)) %>%
select(Comparison, Evidence, MD, SE, "95% CI", p.value)
}Network Graphs
Network Graph for Paper
netgraph(res$net[[1]], labels = long.labels)
Standardized Mean Difference
Oswestry Disability Index
create_dt(tab.list[[1]])Visual Analog Scale: Back
Operation Time
create_dt(tab.list[[2]])Visual Analog Scale: Leg
create_dt(tab.list[[3]])Mean Difference
Length of Stay
create_dt(md.list[[1]])Operation Time
create_dt(md.list[[2]])Blood Loss
create_dt(md.list[[3]])Forest Plots for SMD’s
for(i in 1:nrow(res)){
forest(res$net[[i]],
ref= c("RA", "CT_Nav"),
xlim= c(-1.5, 1.5),
baseline= FALSE,
drop= TRUE,
pooled = "random",
smlab = outcome.labels.smd[i],
labels = long.labels,
label.left = "Favors Treatment",
label.right = "Favors 'Other'")
}
# Length of Stay
forest(res.md$net[[1]],
ref= c("RA", "CT_Nav"),
xlim= c(-5, 5),
baseline= FALSE,
drop= TRUE,
pooled = "random",
smlab = outcome.labels.md[1],
labels = long.labels,
label.left = "Favors Treatment",
label.right = "Favors 'Other'")
# Operation Time
forest(res.md$net[[2]],
ref= c("RA", "CT_Nav"),
xlim= c(-130, 130),
baseline= FALSE,
drop= TRUE,
pooled = "random",
smlab = outcome.labels.md[2],
labels = long.labels,
label.left = "Favors Treatment",
label.right = "Favors 'Other'")
# Blood Loss
forest(res.md$net[[3]],
ref= c("RA", "CT_Nav"),
xlim= c(-250, 150),
baseline= FALSE,
drop= TRUE,
pooled = "random",
smlab = outcome.labels.md[3],
labels = long.labels,
label.left = "Favors Treatment",
label.right = "Favors 'Other'")
P.Scores
plot(
netrank(res$net[[1]]),
netrank(res$net[[2]]),
netrank(res$net[[3]]),
netrank(res.md$net[[1]]),
netrank(res.md$net[[2]]),
netrank(res.md$net[[3]]),
name = outcome.labels2,
digits = 2)
Network Graphs for Supplemental
The numbers on each bar represent the number of studies for each comparison
for(i in 1:nrow(res)){
print(netgraph(res$net[[i]], labels = long.labels, number.of.studies = TRUE), title(main = outcome.labels.smd[i]))
}
## $nodes
## trts labels seq srt xpos ypos xpos.labels
## CT_Nav CT_Nav CT-navigated CT_Nav 0 -0.2886751 0.5 -0.3109031
## FFG FFG Fluoroscopy-guided FFG 0 -0.2886751 -0.5 -0.3109031
## RA RA Robot-assisted RA 0 0.5773503 0.0 0.5995783
## ypos.labels offset.x offset.y cex col pch bg adj.x adj.y
## CT_Nav 0.52222799 0.02222799 0.02222799 1 red 20 red 1 0
## FFG -0.52222799 0.02222799 0.02222799 1 red 20 red 1 1
## RA -0.02222799 0.02222799 0.02222799 1 red 20 red 0 1
##
## $edges
## treat1 treat2 n.stud xpos ypos adj pos.number.of.studies
## CT_Nav vs. FFG CT_Nav FFG 1 -0.2886751 0.00 0.5 0.5
## FFG vs. RA FFG RA 2 0.1443376 -0.25 0.5 0.5
## col
## CT_Nav vs. FFG black
## FFG vs. RA black
## $nodes
## trts labels seq srt xpos ypos xpos.labels
## CT_Nav CT_Nav CT-navigated CT_Nav 0 -0.2886751 0.5 -0.3109031
## FFG FFG Fluoroscopy-guided FFG 0 -0.2886751 -0.5 -0.3109031
## RA RA Robot-assisted RA 0 0.5773503 0.0 0.5995783
## ypos.labels offset.x offset.y cex col pch bg adj.x adj.y
## CT_Nav 0.52222799 0.02222799 0.02222799 1 red 20 red 1 0
## FFG -0.52222799 0.02222799 0.02222799 1 red 20 red 1 1
## RA -0.02222799 0.02222799 0.02222799 1 red 20 red 0 1
##
## $edges
## treat1 treat2 n.stud xpos ypos adj pos.number.of.studies
## CT_Nav vs. FFG CT_Nav FFG 1 -0.2886751 0.00 0.5 0.5
## FFG vs. RA FFG RA 3 0.1443376 -0.25 0.5 0.5
## col
## CT_Nav vs. FFG black
## FFG vs. RA black
## $nodes
## trts labels seq srt xpos ypos xpos.labels
## CT_Nav CT_Nav CT-navigated CT_Nav 0 -0.2886751 0.5 -0.3109031
## FFG FFG Fluoroscopy-guided FFG 0 -0.2886751 -0.5 -0.3109031
## RA RA Robot-assisted RA 0 0.5773503 0.0 0.5995783
## ypos.labels offset.x offset.y cex col pch bg adj.x adj.y
## CT_Nav 0.52222799 0.02222799 0.02222799 1 red 20 red 1 0
## FFG -0.52222799 0.02222799 0.02222799 1 red 20 red 1 1
## RA -0.02222799 0.02222799 0.02222799 1 red 20 red 0 1
##
## $edges
## treat1 treat2 n.stud xpos ypos adj pos.number.of.studies
## CT_Nav vs. FFG CT_Nav FFG 1 -0.2886751 0.00 0.5 0.5
## FFG vs. RA FFG RA 2 0.1443376 -0.25 0.5 0.5
## col
## CT_Nav vs. FFG black
## FFG vs. RA blackfor(i in 1:nrow(res.md)){
print(netgraph(res.md$net[[i]], labels = long.labels, number.of.studies = TRUE), title(main = outcome.labels.md[i]))
}
## $nodes
## trts labels seq srt xpos ypos xpos.labels
## CT_Nav CT_Nav CT-navigated CT_Nav 0 -0.2886751 0.5 -0.3109031
## FFG FFG Fluoroscopy-guided FFG 0 -0.2886751 -0.5 -0.3109031
## RA RA Robot-assisted RA 0 0.5773503 0.0 0.5995783
## ypos.labels offset.x offset.y cex col pch bg adj.x adj.y
## CT_Nav 0.52222799 0.02222799 0.02222799 1 red 20 red 1 0
## FFG -0.52222799 0.02222799 0.02222799 1 red 20 red 1 1
## RA -0.02222799 0.02222799 0.02222799 1 red 20 red 0 1
##
## $edges
## treat1 treat2 n.stud xpos ypos adj pos.number.of.studies
## CT_Nav vs. FFG CT_Nav FFG 2 -0.2886751 0.00 0.5 0.5
## FFG vs. RA FFG RA 2 0.1443376 -0.25 0.5 0.5
## col
## CT_Nav vs. FFG black
## FFG vs. RA black
## $nodes
## trts labels seq srt xpos ypos xpos.labels
## CT_Nav CT_Nav CT-navigated CT_Nav 0 -0.2886751 0.5 -0.3109031
## FFG FFG Fluoroscopy-guided FFG 0 -0.2886751 -0.5 -0.3109031
## RA RA Robot-assisted RA 0 0.5773503 0.0 0.5995783
## ypos.labels offset.x offset.y cex col pch bg adj.x adj.y
## CT_Nav 0.52222799 0.02222799 0.02222799 1 red 20 red 1 0
## FFG -0.52222799 0.02222799 0.02222799 1 red 20 red 1 1
## RA -0.02222799 0.02222799 0.02222799 1 red 20 red 0 1
##
## $edges
## treat1 treat2 n.stud xpos ypos adj pos.number.of.studies
## CT_Nav vs. FFG CT_Nav FFG 2 -0.2886751 0.00 0.5 0.5
## FFG vs. RA FFG RA 3 0.1443376 -0.25 0.5 0.5
## col
## CT_Nav vs. FFG black
## FFG vs. RA black
## $nodes
## trts labels seq srt xpos ypos xpos.labels
## CT_Nav CT_Nav CT-navigated CT_Nav 0 -0.2886751 0.5 -0.3109031
## FFG FFG Fluoroscopy-guided FFG 0 -0.2886751 -0.5 -0.3109031
## RA RA Robot-assisted RA 0 0.5773503 0.0 0.5995783
## ypos.labels offset.x offset.y cex col pch bg adj.x adj.y
## CT_Nav 0.52222799 0.02222799 0.02222799 1 red 20 red 1 0
## FFG -0.52222799 0.02222799 0.02222799 1 red 20 red 1 1
## RA -0.02222799 0.02222799 0.02222799 1 red 20 red 0 1
##
## $edges
## treat1 treat2 n.stud xpos ypos adj pos.number.of.studies
## CT_Nav vs. FFG CT_Nav FFG 2 -0.2886751 0.00 0.5 0.5
## FFG vs. RA FFG RA 2 0.1443376 -0.25 0.5 0.5
## col
## CT_Nav vs. FFG black
## FFG vs. RA blackTable of Comparisons for Length of Stay (MD)
create_dt(md.list[[1]])Table of Comparisons for Operation Time (MD)
create_dt(md.list[[2]])Table of Comparisons for Blood Loss (MD)
create_dt(md.list[[3]])Summary Outputs for Outcomes with SMD’s
lapply(res$net, summary)## [[1]]
## Original data:
##
## treat1 treat2 TE seTE
## Chen CT_Nav FFG -0.3411 0.3010
## Cui FFG RA 0.3942 0.2917
## Feng FFG RA 0.0688 0.2237
##
## Number of treatment arms (by study):
## narms
## Chen 2
## Cui 2
## Feng 2
##
## Results (random effects model):
##
## treat1 treat2 SMD 95%-CI
## Chen CT_Nav FFG -0.3411 [-0.9310; 0.2488]
## Cui FFG RA 0.1892 [-0.1587; 0.5371]
## Feng FFG RA 0.1892 [-0.1587; 0.5371]
##
## Number of studies: k = 3
## Number of pairwise comparisons: m = 3
## Number of observations: o = 173
## Number of treatments: n = 3
## Number of designs: d = 2
##
## Random effects model
##
## Treatment estimate (sm = 'SMD', comparison: other treatments vs 'FFG'):
## SMD 95%-CI z p-value
## CT_Nav -0.3411 [-0.9310; 0.2488] -1.13 0.2571
## FFG . . . .
## RA -0.1892 [-0.5371; 0.1587] -1.07 0.2864
##
## Quantifying heterogeneity / inconsistency:
## tau^2 = 0; tau = 0; I^2 = 0%
##
## Tests of heterogeneity (within designs) and inconsistency (between designs):
## Q d.f. p-value
## Total 0.78 1 0.3760
## Within designs 0.78 1 0.3760
## Between designs 0.00 0 --
##
## [[2]]
## Original data:
##
## treat1 treat2 TE seTE
## Chen CT_Nav FFG -0.4939 0.3033
## Cui FFG RA 0.4617 0.2927
## Feng FFG RA 0.0403 0.2236
## Hyun FFG RA 0.6263 0.2645
##
## Number of treatment arms (by study):
## narms
## Chen 2
## Cui 2
## Feng 2
## Hyun 2
##
## Results (random effects model):
##
## treat1 treat2 SMD 95%-CI
## Chen CT_Nav FFG -0.4939 [-1.2007; 0.2129]
## Cui FFG RA 0.3481 [-0.0178; 0.7139]
## Feng FFG RA 0.3481 [-0.0178; 0.7139]
## Hyun FFG RA 0.3481 [-0.0178; 0.7139]
##
## Number of studies: k = 4
## Number of pairwise comparisons: m = 4
## Number of observations: o = 233
## Number of treatments: n = 3
## Number of designs: d = 2
##
## Random effects model
##
## Treatment estimate (sm = 'SMD', comparison: other treatments vs 'FFG'):
## SMD 95%-CI z p-value
## CT_Nav -0.4939 [-1.2007; 0.2129] -1.37 0.1708
## FFG . . . .
## RA -0.3481 [-0.7139; 0.0178] -1.86 0.0622
##
## Quantifying heterogeneity / inconsistency:
## tau^2 = 0.0380; tau = 0.1950; I^2 = 36.2% [0.0%; 79.6%]
##
## Tests of heterogeneity (within designs) and inconsistency (between designs):
## Q d.f. p-value
## Total 3.14 2 0.2085
## Within designs 3.14 2 0.2085
## Between designs 0.00 0 --
##
## [[3]]
## Original data:
##
## treat1 treat2 TE seTE
## Chen CT_Nav FFG -0.1548 0.2993
## Feng FFG RA -0.0370 0.2236
## Hyun FFG RA 0.2034 0.2589
##
## Number of treatment arms (by study):
## narms
## Chen 2
## Feng 2
## Hyun 2
##
## Results (random effects model):
##
## treat1 treat2 SMD 95%-CI
## Chen CT_Nav FFG -0.1548 [-0.7413; 0.4317]
## Feng FFG RA 0.0657 [-0.2660; 0.3974]
## Hyun FFG RA 0.0657 [-0.2660; 0.3974]
##
## Number of studies: k = 3
## Number of pairwise comparisons: m = 3
## Number of observations: o = 185
## Number of treatments: n = 3
## Number of designs: d = 2
##
## Random effects model
##
## Treatment estimate (sm = 'SMD', comparison: other treatments vs 'FFG'):
## SMD 95%-CI z p-value
## CT_Nav -0.1548 [-0.7413; 0.4317] -0.52 0.6049
## FFG . . . .
## RA -0.0657 [-0.3974; 0.2660] -0.39 0.6978
##
## Quantifying heterogeneity / inconsistency:
## tau^2 = 0; tau = 0; I^2 = 0%
##
## Tests of heterogeneity (within designs) and inconsistency (between designs):
## Q d.f. p-value
## Total 0.49 1 0.4822
## Within designs 0.49 1 0.4822
## Between designs 0.00 0 --Summary Outputs for Outcomes with MD’s
lapply(res.md$net, summary)## [[1]]
## Original data:
##
## treat1 treat2 TE seTE
## Chen CT_Nav FFG -3.2200 0.9659
## Cui FFG RA 2.7000 0.4932
## Hyun FFG RA 2.6000 1.0578
## Tian CT_Nav FFG -1.0500 0.3084
##
## Number of treatment arms (by study):
## narms
## Chen 2
## Cui 2
## Hyun 2
## Tian 2
##
## Results (random effects model):
##
## treat1 treat2 MD 95%-CI
## Chen CT_Nav FFG -1.7777 [-3.2537; -0.3018]
## Cui FFG RA 2.6652 [ 1.0814; 4.2490]
## Hyun FFG RA 2.6652 [ 1.0814; 4.2490]
## Tian CT_Nav FFG -1.7777 [-3.2537; -0.3018]
##
## Number of studies: k = 4
## Number of pairwise comparisons: m = 4
## Number of observations: o = 214
## Number of treatments: n = 3
## Number of designs: d = 2
##
## Random effects model
##
## Treatment estimate (sm = 'MD', comparison: other treatments vs 'FFG'):
## MD 95%-CI z p-value
## CT_Nav -1.7777 [-3.2537; -0.3018] -2.36 0.0182
## FFG . . . .
## RA -2.6652 [-4.2490; -1.0814] -3.30 0.0010
##
## Quantifying heterogeneity / inconsistency:
## tau^2 = 0.7581; tau = 0.8707; I^2 = 56.4% [0.0%; 87.6%]
##
## Tests of heterogeneity (within designs) and inconsistency (between designs):
## Q d.f. p-value
## Total 4.59 2 0.1009
## Within designs 4.59 2 0.1009
## Between designs 0.00 0 --
##
## [[2]]
## Original data:
##
## treat1 treat2 TE seTE
## Chen CT_Nav FFG 87.8300 8.0816
## Cui FFG RA -32.9000 2.7332
## Feng FFG RA 34.3800 13.2115
## Hyun FFG RA -0.0000 16.6885
## Tian CT_Nav FFG 46.1400 5.5535
##
## Number of treatment arms (by study):
## narms
## Chen 2
## Cui 2
## Feng 2
## Hyun 2
## Tian 2
##
## Results (random effects model):
##
## treat1 treat2 MD 95%-CI
## Chen CT_Nav FFG 66.6774 [ 19.3217; 114.0330]
## Cui FFG RA -1.1049 [-41.3228; 39.1129]
## Feng FFG RA -1.1049 [-41.3228; 39.1129]
## Hyun FFG RA -1.1049 [-41.3228; 39.1129]
## Tian CT_Nav FFG 66.6774 [ 19.3217; 114.0330]
##
## Number of studies: k = 5
## Number of pairwise comparisons: m = 5
## Number of observations: o = 294
## Number of treatments: n = 3
## Number of designs: d = 2
##
## Random effects model
##
## Treatment estimate (sm = 'MD', comparison: other treatments vs 'FFG'):
## MD 95%-CI z p-value
## CT_Nav 66.6774 [ 19.3217; 114.0330] 2.76 0.0058
## FFG . . . .
## RA 1.1049 [-39.1129; 41.3228] 0.05 0.9571
##
## Quantifying heterogeneity / inconsistency:
## tau^2 = 1119.7305; tau = 33.4624; I^2 = 93.5% [86.5%; 96.9%]
##
## Tests of heterogeneity (within designs) and inconsistency (between designs):
## Q d.f. p-value
## Total 46.13 3 < 0.0001
## Within designs 46.13 3 < 0.0001
## Between designs 0.00 0 --
##
## [[3]]
## Original data:
##
## treat1 treat2 TE seTE
## Chen CT_Nav FFG -176.9100 11.4643
## Cui FFG RA 158.5000 6.0017
## Feng FFG RA 72.5000 31.0066
## Tian CT_Nav FFG -89.0200 23.7073
##
## Number of treatment arms (by study):
## narms
## Chen 2
## Cui 2
## Feng 2
## Tian 2
##
## Results (random effects model):
##
## treat1 treat2 MD 95%-CI
## Chen CT_Nav FFG -135.4998 [-220.0307; -50.9690]
## Cui FFG RA 120.6217 [ 34.8582; 206.3853]
## Feng FFG RA 120.6217 [ 34.8582; 206.3853]
## Tian CT_Nav FFG -135.4998 [-220.0307; -50.9690]
##
## Number of studies: k = 4
## Number of pairwise comparisons: m = 4
## Number of observations: o = 234
## Number of treatments: n = 3
## Number of designs: d = 2
##
## Random effects model
##
## Treatment estimate (sm = 'MD', comparison: other treatments vs 'FFG'):
## MD 95%-CI z p-value
## CT_Nav -135.4998 [-220.0307; -50.9690] -3.14 0.0017
## FFG . . . .
## RA -120.6217 [-206.3853; -34.8582] -2.76 0.0058
##
## Quantifying heterogeneity / inconsistency:
## tau^2 = 3385.8700; tau = 58.1882; I^2 = 89.2% [70.7%; 96.0%]
##
## Tests of heterogeneity (within designs) and inconsistency (between designs):
## Q d.f. p-value
## Total 18.55 2 < 0.0001
## Within designs 18.55 2 < 0.0001
## Between designs 0.00 0 --Inconsistency and Heterogeneity
LOS
res.md$incon[[1]]## Q statistics to assess homogeneity / consistency
##
## Q df p-value
## Total 4.59 2 0.1009
## Within designs 4.59 2 0.1009
## Between designs 0.00 0 --
##
## Design-specific decomposition of within-designs Q statistic
##
## Design Q df p-value
## FFG vs. CT_Nav 4.58 1 0.0323
## FFG vs. RA 0.01 1 0.9317
##
## Q statistic to assess consistency under the assumption of
## a full design-by-treatment interaction random effects model
##
## Q df p-value tau.within tau2.within
## Between designs 0.00 0 -- 0.8707 0.7581OP
res.md$incon[[2]]## Q statistics to assess homogeneity / consistency
##
## Q df p-value
## Total 46.13 3 < 0.0001
## Within designs 46.13 3 < 0.0001
## Between designs 0.00 0 --
##
## Design-specific decomposition of within-designs Q statistic
##
## Design Q df p-value
## FFG vs. RA 28.05 2 < 0.0001
## FFG vs. CT_Nav 18.08 1 < 0.0001
##
## Q statistic to assess consistency under the assumption of
## a full design-by-treatment interaction random effects model
##
## Q df p-value tau.within tau2.within
## Between designs 0.00 0 -- 33.4624 1119.7305ODI
res$incon[[1]]## Q statistics to assess homogeneity / consistency
##
## Q df p-value
## Total 0.78 1 0.3760
## Within designs 0.78 1 0.3760
## Between designs 0.00 0 --
##
## Design-specific decomposition of within-designs Q statistic
##
## Design Q df p-value
## FFG vs. RA 0.78 1 0.3760
##
## Q statistic to assess consistency under the assumption of
## a full design-by-treatment interaction random effects model
##
## Q df p-value tau.within tau2.within
## Between designs 0.00 0 -- 0 0VAS: Back
res$incon[[2]]## Q statistics to assess homogeneity / consistency
##
## Q df p-value
## Total 3.14 2 0.2085
## Within designs 3.14 2 0.2085
## Between designs 0.00 0 --
##
## Design-specific decomposition of within-designs Q statistic
##
## Design Q df p-value
## FFG vs. RA 3.14 2 0.2085
##
## Q statistic to assess consistency under the assumption of
## a full design-by-treatment interaction random effects model
##
## Q df p-value tau.within tau2.within
## Between designs 0.00 0 -- 0.1950 0.0380Vas: Leg
res$incon[[3]]## Q statistics to assess homogeneity / consistency
##
## Q df p-value
## Total 0.49 1 0.4822
## Within designs 0.49 1 0.4822
## Between designs 0.00 0 --
##
## Design-specific decomposition of within-designs Q statistic
##
## Design Q df p-value
## FFG vs. RA 0.49 1 0.4822
##
## Q statistic to assess consistency under the assumption of
## a full design-by-treatment interaction random effects model
##
## Q df p-value tau.within tau2.within
## Between designs 0.00 0 -- 0 0Blood Loss
res.md$incon[[3]]## Q statistics to assess homogeneity / consistency
##
## Q df p-value
## Total 18.55 2 < 0.0001
## Within designs 18.55 2 < 0.0001
## Between designs 0.00 0 --
##
## Design-specific decomposition of within-designs Q statistic
##
## Design Q df p-value
## FFG vs. CT_Nav 11.14 1 0.0008
## FFG vs. RA 7.42 1 0.0065
##
## Q statistic to assess consistency under the assumption of
## a full design-by-treatment interaction random effects model
##
## Q df p-value tau.within tau2.within
## Between designs 0.00 0 -- 58.1882 3385.8700Bayesian Network Meta-Regression
Load Packages
library(gemtc)## Warning: package 'gemtc' was built under R version 4.3.3## Loading required package: coda##
## Attaching package: 'gemtc'## The following object is masked from 'package:meta':
##
## forestlibrary(rjags)## Linked to JAGS 4.3.1## Loaded modules: basemod,bugs# Data Wrangling for Bayesian format
## Outcomes with SMD Effect Sizes
out.smd <- c("ODI", "VAS_Back", "VAS_Leg")
treat.codes = tibble(id = c("CT_Nav", "FFG", "RA"),
description = c("CT-navigation", "Fluoroscopy-guided", "Robotic-assisted"))
bayes.df.clean <- function(df){
bayes.df <- df %>%
select(1:3) %>%
unnest() %>%
select(1:5, TE, seTE) %>%
pivot_longer(cols = 4:5,
names_to = "tx",
values_to = "treatment") %>%
mutate(diff = ifelse(treatment == "FFG", NA, TE),
std.err = ifelse(treatment == "FFG", NA, seTE),
study = auth,
randomized = ifelse(design == "rct", 1, 0))
data <- bayes.df %>%
select(1, 10, 8:9, 7) %>%
nest(!outcome)
study.info <- bayes.df %>%
select(1, 10:11) %>%
group_by(outcome) %>%
unique() %>%
nest(!outcome)
res.bayes <- merge(data, study.info, by = "outcome") %>%
mutate(data = data.x,
study.info = data.y) %>%
select(1, 4:5)
}
md.bayes.smd <- bayes.df.clean(res)## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(data, pairwise)`.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.md.bayes.md <- bayes.df.clean(res.md)## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(data, pairwise)`.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.md.bayes.smd <- bayes.df.clean(res)## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(data, pairwise)`.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.md.bayes.md <- bayes.df.clean(res.md)## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(data, pairwise)`.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = !outcome`?
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.Create Bayesian Meta-regression Model for SMD and MD Outcomes
Network Meta-regression Function
regressor <- list(coefficient = "shared",
variable = "randomized",
control = "FFG")
bayes.meta.regression.func <- function(dataframe){
bayes <- dataframe %>%
mutate(network.mr = map2(.x = data,
.y = study.info,
~mtc.network(data.re = .x,
studies = .y,
treatments = treat.codes)),
network.null = map(data, ~mtc.network(data.re = .x,
treatments = treat.codes)),
sum.net = map(network.mr, ~summary(.x)),
null.sum.net = map(network.null, ~summary(.x)),
# # Network Meta-regression Model Compilation
model.mr = map(network.mr, ~mtc.model(.x,
likelihood = "normal",
link = "identity",
type = "regression",
linearModel = "random",
regressor = regressor,
n.chain = 4)),
model.null = map(network.null,
~mtc.model(.x,
likelihood = "normal",
link = "identity",
linearModel = "random",
n.chain = 4)),
## Markov Chain Monte Carlo Sampling ##
# NMR with Randomization Covariate #
mcmc1 = map(model.mr,
~mtc.run(.x,
n.adapt = 50,
n.iter = 1000,
thin = 10)),
mcmc3 = map(model.mr, ~mtc.run(.x,
n.adapt = 5000,
n.iter = 1e6,
thin = 10)),
# NMR Null Model #
null.mcmc1 = map(model.null,
~mtc.run(.x,
n.adapt = 50,
n.iter = 1000,
thin = 10)),
null.mcmc2 = map(model.null,
~mtc.run(.x,
n.adapt = 5000,
n.iter = 1e6,
thin = 10)),
## Gelman-Rubin plot: Potential Scale Reduction Factor (PSRF) ##
cov.psrf.mcmc1 = map(mcmc1, ~gelman.diag(.x)$mpsrf),
cov.psrf.mcmc3 = map(mcmc3, ~gelman.diag(.x)$mpsrf),
null.psrf.mcmc1 = map(null.mcmc1, ~gelman.diag(.x)$mpsrf),
null.psrf.mcmc2 = map(null.mcmc2, ~gelman.diag(.x)$mpsrf))
for(i in 2:length(bayes)){
bayes[[i]] <- setNames(bayes[[i]], bayes[[1]])
}
return(bayes)
}Run Meta-regression Models
Checking Convergence of NMR
SMD Outcome
NMR Models - Covariate Models
rbind(Low_Itr = smd.nmr.final$cov.psrf.mcmc1, High_Itr = smd.nmr.final$cov.psrf.mcmc3)## ODI VAS_Back VAS_Leg
## Low_Itr 1.26748 1.452925 2.586471
## High_Itr 1.00258 1.000027 1.000137NMR Models - Null Models
rbind(Low_Itr = smd.nmr.final$null.psrf.mcmc1, High_Itr = smd.nmr.final$null.psrf.mcmc2)## ODI VAS_Back VAS_Leg
## Low_Itr 1.028621 1.036124 1.008958
## High_Itr 1.000036 1.000025 1.000053MD Outcomes
NMR Models - Covariate Models
rbind(Low_Itr = md.nmr.final$cov.psrf.mcmc1, High_Itr = md.nmr.final$cov.psrf.mcmc3)## Blood_Loss LOS OP_Time
## Low_Itr 10.31206 1.096186 4.108816
## High_Itr 1.027022 1.018334 1.032859SMD Outcome NMR Models - Null Models
rbind(Low_Itr = md.nmr.final$null.psrf.mcmc1, High_Itr = md.nmr.final$null.psrf.mcmc2)## Blood_Loss LOS OP_Time
## Low_Itr 2.130358 1.012832 1.045029
## High_Itr 1.000014 1.000011 1.000015DIC Tables for Moderator Analysis:
Model <- c("Null", "Covariate")
smd.dic.cov <- map_dbl(smd.nmr.final$mcmc3, ~summary(.x)$DIC["DIC"])
smd.dic.null <- map_dbl(smd.nmr.final$null.mcmc2, ~summary(.x)$DIC["DIC"])
smd.dic.comb <- as.data.frame(rbind(smd.dic.null, smd.dic.cov)) %>% round(., digits = 2)
smd.dic.df <- tibble(Model, smd.dic.comb)
md.dic.cov <- map_dbl(md.nmr.final$mcmc3, ~summary(.x)$DIC["DIC"])
md.dic.null <- map_dbl(md.nmr.final$null.mcmc2, ~summary(.x)$DIC["DIC"])
md.dic.comb <- as.data.frame(rbind(md.dic.null, md.dic.cov)) %>% round(., digits = 2)
md.dic.df <- tibble(Model, md.dic.comb)DIC Table for SMD Outcomes
create_dt(smd.dic.df)DIC Table for MD Outcomes
create_dt(md.dic.df)Summary data for table function
smd.test.list <- smd.nmr.final %>%
mutate(test.cov = map(mcmc3, ~summary(.x)),
test.null = map(null.mcmc2, ~summary(.x)))
md.test.list <- md.nmr.final %>%
mutate(test.cov = map(mcmc3, ~summary(.x)),
test.null = map(null.mcmc2, ~summary(.x)))
smd.list.cov <- smd.test.list$test.cov
smd.list.null <- smd.test.list$test.null
md.list.cov <- md.test.list$test.cov
md.list.null <- md.test.list$test.nullSummaries for Moderator Analyses
SMD Outcomes
Covariate Models
map(smd.nmr.final$mcmc3, ~summary(.x))## $ODI
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -0.30897 1.1160 0.0017645 0.013771
## d.FFG.RA -0.21837 0.5779 0.0009138 0.006996
## sd.d 0.18597 0.1124 0.0001778 0.000231
## B 0.04612 1.5841 0.0025047 0.020673
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -2.353735 -0.73476 -0.335043 0.06717 1.9612
## d.FFG.RA -1.379018 -0.43586 -0.204492 0.02260 0.8349
## sd.d 0.008631 0.08851 0.180273 0.28073 0.3820
## B -2.916705 -0.36416 0.004489 0.38020 3.3568
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 2.722435 2.346363 5.068798
##
## 3 data points, ratio 0.9075, I^2 = 27%
##
## -- Regression settings:
##
## Regression on "randomized", shared coefficients, "FFG" as control
## Input standardized: x' = (randomized - 0.6666667) / 1
## Estimates at the centering value: randomized = 0.6666667
##
##
## $VAS_Back
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -0.46456 1.8458 0.0029184 0.0272480
## d.FFG.RA -0.36267 0.6453 0.0010204 0.0091343
## sd.d 0.28423 0.1698 0.0002685 0.0003727
## B 0.03948 2.3927 0.0037831 0.0364579
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -4.04642 -1.0808 -0.491316 0.1068 3.3646
## d.FFG.RA -1.67484 -0.6006 -0.352356 -0.1086 0.8733
## sd.d 0.01526 0.1426 0.271548 0.4186 0.6009
## B -4.66280 -0.5855 0.004291 0.5975 5.0639
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 3.971923 3.172315 7.144238
##
## 4 data points, ratio 0.993, I^2 = 24%
##
## -- Regression settings:
##
## Regression on "randomized", shared coefficients, "FFG" as control
## Input standardized: x' = (randomized - 0.75) / 1
## Estimates at the centering value: randomized = 0.75
##
##
## $VAS_Leg
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -0.1533204 0.6204 9.810e-04 0.005485
## d.FFG.RA -0.0682466 0.3279 5.184e-04 0.002743
## sd.d 0.0989038 0.0583 9.218e-05 0.000122
## B 0.0002268 0.8074 1.277e-03 0.007932
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -1.326577 -0.44014 -0.1531796 0.13336 1.0195
## d.FFG.RA -0.688973 -0.22932 -0.0676445 0.09184 0.5486
## sd.d 0.004753 0.04838 0.0975540 0.14860 0.1978
## B -1.602883 -0.19095 0.0001728 0.19263 1.6029
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 2.502350 2.157040 4.659389
##
## 3 data points, ratio 0.8341, I^2 = 20%
##
## -- Regression settings:
##
## Regression on "randomized", shared coefficients, "FFG" as control
## Input standardized: x' = (randomized - 0.6666667) / 1
## Estimates at the centering value: randomized = 0.6666667Null Models
map(smd.nmr.final$null.mcmc2, ~summary(.x))## $ODI
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -0.3393 0.3702 0.0005854 0.0006913
## d.FFG.RA -0.2030 0.2362 0.0003734 0.0004210
## sd.d 0.1861 0.1122 0.0001773 0.0002326
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -1.073883 -0.58106 -0.3386 -0.09678 0.3922
## d.FFG.RA -0.678376 -0.35342 -0.2010 -0.05049 0.2609
## sd.d 0.008796 0.08914 0.1808 0.28060 0.3819
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 2.712409 2.359345 5.071754
##
## 3 data points, ratio 0.9041, I^2 = 26%
##
##
## $VAS_Back
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -0.4918 0.4481 0.0007085 0.0007702
## d.FFG.RA -0.3529 0.2421 0.0003827 0.0004037
## sd.d 0.2841 0.1695 0.0002681 0.0003698
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -1.3937 -0.7707 -0.4922 -0.2134 0.4126
## d.FFG.RA -0.8546 -0.4975 -0.3497 -0.2049 0.1311
## sd.d 0.0153 0.1424 0.2717 0.4178 0.6006
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 3.973884 3.187625 7.161509
##
## 4 data points, ratio 0.9935, I^2 = 25%
##
##
## $VAS_Leg
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -0.15198 0.31909 5.045e-04 0.0007759
## d.FFG.RA -0.06908 0.18763 2.967e-04 0.0004100
## sd.d 0.09885 0.05834 9.225e-05 0.0001223
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -0.779703 -0.36596 -0.15212 0.06225 0.4745
## d.FFG.RA -0.438422 -0.19424 -0.06890 0.05642 0.2990
## sd.d 0.004798 0.04836 0.09742 0.14869 0.1977
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 2.508730 2.161296 4.670026
##
## 3 data points, ratio 0.8362, I^2 = 20%MD Outcomes
Covariate Models
map(md.nmr.final$mcmc3, ~summary(.x))## $Blood_Loss
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -118.81 327.06 0.51713 24.3867
## d.FFG.RA -135.72 327.50 0.51782 24.3884
## sd.d 87.94 39.44 0.06236 0.0634
## B 32.22 640.26 1.01234 50.2270
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -686.56 -238.74 -137.980 -39.23 618.9
## d.FFG.RA -876.03 -215.34 -117.111 -15.73 432.7
## sd.d 27.97 56.07 81.811 116.49 168.7
## B -1079.63 -167.99 -5.057 152.92 1489.2
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 4.12005 -5.94135 -1.82130
##
## 4 data points, ratio 1.03, I^2 = 27%
##
## -- Regression settings:
##
## Regression on "randomized", shared coefficients, "FFG" as control
## Input standardized: x' = (randomized - 0.5) / 1
## Estimates at the centering value: randomized = 0.5
##
##
## $LOS
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -2.3063 7.5175 0.011886 0.32427
## d.FFG.RA -2.1822 7.5226 0.011894 0.32537
## sd.d 1.3524 0.8209 0.001298 0.00168
## B -0.9589 14.8428 0.023469 0.65438
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -19.39744 -3.7750 -1.84656 -0.04408 9.286
## d.FFG.RA -13.72760 -4.4546 -2.61441 -0.71098 14.949
## sd.d 0.09139 0.6945 1.23581 1.94414 3.039
## B -35.09991 -3.1783 -0.07685 2.91825 21.957
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 4.1698871 -0.1978215 3.9720656
##
## 4 data points, ratio 1.042, I^2 = 28%
##
## -- Regression settings:
##
## Regression on "randomized", shared coefficients, "FFG" as control
## Input standardized: x' = (randomized - 0.5) / 1
## Estimates at the centering value: randomized = 0.5
##
##
## $OP_Time
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav 96.31 184.84 0.29226 10.73701
## d.FFG.RA -18.83 124.44 0.19676 7.25866
## sd.d 44.20 17.66 0.02792 0.02803
## B 49.41 303.08 0.47921 17.96345
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -195.65 17.71 71.989 132.95 644.91
## d.FFG.RA -383.76 -46.09 -2.731 35.90 178.38
## sd.d 18.11 30.25 40.988 56.10 82.79
## B -428.78 -68.12 7.165 94.12 961.12
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 4.992124 -93.022198 -88.030074
##
## 5 data points, ratio 0.9984, I^2 = 20%
##
## -- Regression settings:
##
## Regression on "randomized", shared coefficients, "FFG" as control
## Input standardized: x' = (randomized - 0.6) / 1
## Estimates at the centering value: randomized = 0.6Null Models
map(md.nmr.final$null.mcmc2, ~summary(.x))## $Blood_Loss
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -135.13 69.45 0.10981 0.10970
## d.FFG.RA -119.54 69.90 0.11052 0.11017
## sd.d 88.01 39.49 0.06245 0.06338
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -279.44 -172.73 -136.16 -97.78 11.80
## d.FFG.RA -263.34 -158.09 -121.37 -81.57 29.31
## sd.d 28.07 56.03 81.71 116.78 168.80
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 4.120806 3.972875 8.093681
##
## 4 data points, ratio 1.03, I^2 = 27%
##
##
## $LOS
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav -1.827 1.2304 0.001945 0.002012
## d.FFG.RA -2.662 1.2379 0.001957 0.001963
## sd.d 1.353 0.8207 0.001298 0.001695
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -4.60747 -2.4572 -1.703 -1.135 0.54951
## d.FFG.RA -5.27503 -3.3166 -2.665 -2.009 -0.04322
## sd.d 0.08909 0.6959 1.238 1.943 3.04026
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 4.163037 3.611389 7.774426
##
## 4 data points, ratio 1.041, I^2 = 28%
##
##
## $OP_Time
##
## Results on the Mean Difference scale
##
## Iterations = 5010:1005000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.FFG.CT_Nav 66.5872 33.95 0.05368 0.05375
## d.FFG.RA 0.9265 28.34 0.04481 0.04482
## sd.d 44.2135 17.64 0.02790 0.02822
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.FFG.CT_Nav -3.773 47.47 66.56 85.65 137.20
## d.FFG.RA -58.617 -15.12 1.40 17.32 58.42
## sd.d 18.104 30.27 41.06 56.04 82.76
##
## -- Model fit (residual deviance):
##
## Dbar pD DIC
## 4.993603 4.852413 9.846016
##
## 5 data points, ratio 0.9987, I^2 = 20%Forest Plots for the Moderator Analyses
SMD Outcomes
Randomization = 1 Models
map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))
## $ODI
## NULL
##
## $VAS_Back
## NULL
##
## $VAS_Leg
## NULLRandomization = 0
map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))
## $ODI
## NULL
##
## $VAS_Back
## NULL
##
## $VAS_Leg
## NULLMD Outcomes
Randomization = 1 Models
map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))
## $Blood_Loss
## NULL
##
## $LOS
## NULL
##
## $OP_Time
## NULLRandomization = 0
map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))
## $Blood_Loss
## NULL
##
## $LOS
## NULL
##
## $OP_Time
## NULLFunction for creating Tables for Meta-regression
tables.func <- function(tab, mod){
list.names <- names(tab)
list_tabs <- list()
if(mod == "COV"){
for(i in 1:length(tab)){
zz <- tab[[i]]$summaries$statistics %>%
round(., digits = 2) %>%
as.data.frame(., row.names = row.names(.)) %>%
rownames_to_column(., var = "rownames") %>%
as_tibble()
xx <- tab[[i]]$summaries$quantiles %>%
round(., digits = 2) %>%
as.data.frame(., row.names = row.names(.)) %>%
rownames_to_column(., var = "rownames") %>%
as_tibble()
yy <- cbind(zz, xx) %>%
select(1:3, 7, 11) %>%
filter(rownames != "sd.d") %>%
mutate("95% CrI" = c(paste(.[[4]], .[[5]], sep = ", ")),
Estimates = case_when(
rownames == "d.FFG.CT_Nav" ~ "CT-Nav v.s. FFG",
rownames == "d.FFG.RA" ~ "RA v.s. FFG",
rownames == "B" ~ "Randomized",
TRUE ~ NA
)) %>%
select(Estimates, Mean, "95% CrI")
list_tabs[[i]] <- yy
}
} else{
for(i in 1:length(tab)){
zz <- tab[[i]]$summaries$statistics %>%
round(., digits = 2) %>%
as.data.frame(., row.names = row.names(.)) %>%
rownames_to_column(., var = "rownames") %>%
as_tibble()
xx <- tab[[i]]$summaries$quantiles %>%
round(., digits = 2) %>%
as.data.frame(., row.names = row.names(.)) %>%
rownames_to_column(., var = "rownames") %>%
as_tibble()
yy <- cbind(zz, xx) %>%
select(1:3, 7, 11) %>%
filter(rownames != "sd.d") %>%
mutate("95% CrI" = c(paste(.[[4]], .[[5]], sep = ", ")),
Estimates = case_when(
rownames == "d.FFG.CT_Nav" ~ "CT-Nav v.s. FFG",
rownames == "d.FFG.RA" ~ "RA v.s. FFG",
TRUE ~ NA
)) %>%
select(Estimates, Mean, "95% CrI")
list_tabs[[i]] <- yy
}
}
names(list_tabs) <- list.names
return(list_tabs)
}Tables of Regression Outcomes for SMD Outcomes
smd.Final.Tables.null <- tables.func(smd.list.null, "NULL")
smd.Final.Tables.cov <- tables.func(smd.list.cov, "COV")ODI
Null
create_dt(smd.Final.Tables.null$ODI)Covariate Model
create_dt(smd.Final.Tables.cov$ODI)VAS_Back
Null
create_dt(smd.Final.Tables.null$VAS_Back)Covariate Model
create_dt(smd.Final.Tables.cov$VAS_Back)VAS_Leg
Null
create_dt(smd.Final.Tables.null$VAS_Leg)Covariate Model
create_dt(smd.Final.Tables.cov$VAS_Leg)Tables of Regression Outcomes for MD Outcomes
md.Final.Tables.null <- tables.func(md.list.null, "NULL")
md.Final.Tables.cov <- tables.func(md.list.cov, "COV")Blood_Loss
Null
create_dt(md.Final.Tables.null$Blood_Loss)Covariate Model
create_dt(md.Final.Tables.cov$Blood_Loss)LOS
Null
create_dt(md.Final.Tables.null$LOS)Covariate Model
create_dt(md.Final.Tables.cov$LOS)OP_Time
Null
create_dt(md.Final.Tables.null$OP_Time)Covariate Model
create_dt(md.Final.Tables.cov$OP_Time)Comparison of Study Characteristics by Treatments
Since the assessment of transitivity is unachievable given the star network (a network where all treatments have been compared with a common treatment but not between themselves) nature of our NMA, distributions of study characteristics related to each intervention will be assessed graphically.
Boxplots of Sample Size
data %>%
ggplot(aes(x = trt, y = n, fill = trt)) +
geom_boxplot() +
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 2) +
coord_cartesian(ylim=c(10,50)) +
scale_fill_manual(values = c("#0099f8", "#e74c3c", "#2ecc71"))+
facet_wrap(~outcome)## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
Age
data %>%
ggplot(aes(x = trt, y = age, fill = trt)) +
geom_boxplot() +
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 2) +
coord_cartesian(ylim=c(40,80)) +
scale_fill_manual(values = c("#0099f8", "#e74c3c", "#2ecc71"))+
facet_wrap(~outcome)## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
Baseline VAS: Back and Leg
data %>%
filter(outcome %in% c("VAS_Back", "VAS_Leg")) %>%
na.omit() %>%
ggplot(aes(x = trt, y = baseline, fill = trt)) +
geom_boxplot() +
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 2) +
coord_cartesian(ylim=c(3,10)) +
scale_fill_manual(values = c("#0099f8", "#e74c3c", "#2ecc71"))+
facet_wrap(~outcome)## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
Baseline ODI
data %>%
filter(outcome == "ODI") %>%
na.omit() %>%
ggplot(aes(x = trt, y = baseline, fill = trt)) +
geom_boxplot() +
coord_cartesian(ylim=c(20,80)) +
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 2) +
scale_fill_manual(values = c("#0099f8", "#e74c3c", "#2ecc71"))## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.