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 errors
library(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/kyz6wjbb
library(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':
## 
##     transpose

Interactive 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     black
for(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     black

Table 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.7581

OP

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.7305

ODI

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           0

VAS: 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.0380

Vas: 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           0

Blood 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.8700

Bayesian 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':
## 
##     forest
library(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 = 2e5, 
                                  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 = 12e5, 
                              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.059176 1.048178 1.019067
## High_Itr 1.001879 1.008846 1.001244

NMR 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.035057 1.005784 1.034434
## High_Itr 1.000029 1.000014 1.000059

MD 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  12.36024   3.672858 5.920863
## High_Itr 1.013081   1.004675 1.011114

SMD 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  1.82975    1.014192 1.379404
## High_Itr 1.000022   1.000017 1.000009

DIC 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)

Summaries for Moderator Analyses

SMD Outcomes

Covariate Models

map(smd.nmr.final$mcmc3, ~summary(.x))
## $ODI
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.3272 1.0255 0.0036256      0.0227220
## d.FFG.RA     -0.2095 0.5339 0.0018876      0.0115603
## sd.d          0.1863 0.1122 0.0003968      0.0005162
## B             0.0192 1.4417 0.0050972      0.0339714
## 
## 2. Quantiles for each variable:
## 
##                   2.5%     25%        50%     75%  97.5%
## d.FFG.CT_Nav -2.280603 -0.7354 -3.383e-01 0.05619 1.6762
## d.FFG.RA     -1.250426 -0.4312 -2.004e-01 0.02461 0.7984
## sd.d          0.009231  0.0890  1.811e-01 0.28091 0.3818
## B            -2.815142 -0.3637  2.836e-05 0.36134 2.9223
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.727472 2.366326 5.093798 
## 
## 3 data points, ratio 0.9092, 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.47775 2.1995 0.0077765      0.0910089
## d.FFG.RA     -0.35788 0.7581 0.0026802      0.0301643
## sd.d          0.28157 0.1697 0.0005999      0.0008392
## B             0.02045 2.8780 0.0101754      0.1209809
## 
## 2. Quantiles for each variable:
## 
##                  2.5%     25%       50%     75%  97.5%
## d.FFG.CT_Nav -5.08506 -1.0955 -0.496884  0.1035 3.8816
## d.FFG.RA     -1.83000 -0.6010 -0.349015 -0.1035 1.1960
## sd.d          0.01473  0.1401  0.268934  0.4153 0.5995
## B            -6.06251 -0.6044 -0.001689  0.5991 5.7147
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 3.984558 3.171019 7.155577 
## 
## 4 data points, ratio 0.9961, I^2 = 25%
## 
## -- 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.13560 0.61362 0.0021695      0.0116643
## d.FFG.RA     -0.07629 0.32436 0.0011468      0.0058630
## sd.d          0.09872 0.05844 0.0002066      0.0002703
## B             0.02558 0.79685 0.0028173      0.0170848
## 
## 2. Quantiles for each variable:
## 
##                   2.5%      25%      50%    75%  97.5%
## d.FFG.CT_Nav -1.239307 -0.43546 -0.14733 0.1426 1.0900
## d.FFG.RA     -0.715765 -0.23130 -0.07033 0.0911 0.5090
## sd.d          0.004609  0.04816  0.09714 0.1486 0.1978
## B            -1.455115 -0.18628  0.00232 0.1965 1.7281
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.516263 2.159614 4.675877 
## 
## 3 data points, ratio 0.8388, I^2 = 21%
## 
## -- 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

Null Models

map(smd.nmr.final$null.mcmc2, ~summary(.x))
## $ODI
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.3397 0.3710 0.0005355      0.0006318
## d.FFG.RA     -0.2030 0.2362 0.0003410      0.0003806
## sd.d          0.1862 0.1122 0.0001620      0.0002131
## 
## 2. Quantiles for each variable:
## 
##                   2.5%      25%     50%      75%  97.5%
## d.FFG.CT_Nav -1.072461 -0.58294 -0.3396 -0.09698 0.3952
## d.FFG.RA     -0.679161 -0.35371 -0.2013 -0.05059 0.2631
## sd.d          0.008902  0.08907  0.1806  0.28037 0.3821
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.720692 2.366518 5.087209 
## 
## 3 data points, ratio 0.9069, I^2 = 26%
## 
## 
## $VAS_Back
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.4932 0.4493 0.0006486      0.0007031
## d.FFG.RA     -0.3524 0.2427 0.0003503      0.0003699
## sd.d          0.2843 0.1698 0.0002451      0.0003381
## 
## 2. Quantiles for each variable:
## 
##                  2.5%     25%     50%     75%  97.5%
## d.FFG.CT_Nav -1.39260 -0.7735 -0.4933 -0.2137 0.4124
## d.FFG.RA     -0.85596 -0.4975 -0.3487 -0.2037 0.1296
## sd.d          0.01527  0.1426  0.2716  0.4185 0.6009
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 3.979676 3.190921 7.170597 
## 
## 4 data points, ratio 0.9949, I^2 = 25%
## 
## 
## $VAS_Leg
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.15295 0.31875 4.601e-04      0.0007049
## d.FFG.RA     -0.06904 0.18775 2.710e-04      0.0003687
## sd.d          0.09901 0.05824 8.407e-05      0.0001104
## 
## 2. Quantiles for each variable:
## 
##                   2.5%      25%      50%     75%  97.5%
## d.FFG.CT_Nav -0.778450 -0.36686 -0.15279 0.06107 0.4725
## d.FFG.RA     -0.439601 -0.19404 -0.06880 0.05649 0.2996
## sd.d          0.004942  0.04858  0.09765 0.14873 0.1977
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.506104 2.157577 4.663681 
## 
## 3 data points, ratio 0.8354, I^2 = 20%

MD Outcomes

Covariate Models

map(md.nmr.final$mcmc3, ~summary(.x))
## $Blood_Loss
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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 -123.35 211.52   0.7478        15.7230
## d.FFG.RA     -131.86 211.65   0.7483        15.6667
## sd.d           87.94  39.57   0.1399         0.1415
## B              23.80 400.00   1.4142        30.4874
## 
## 2. Quantiles for each variable:
## 
##                 2.5%     25%      50%    75% 97.5%
## d.FFG.CT_Nav -548.69 -230.25 -130.388 -29.56 344.2
## d.FFG.RA     -600.24 -225.96 -125.634 -24.48 292.0
## sd.d           27.96   55.96   81.642 116.73 168.8
## B            -798.87 -151.88    9.121 170.67 921.5
## 
## -- Model fit (residual deviance):
## 
##      Dbar        pD       DIC 
##  4.118089 -1.425021  2.693068 
## 
## 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.6685  5.5340 0.019566       0.311619
## d.FFG.RA     -2.8177  5.5304 0.019553       0.309232
## sd.d          1.3549  0.8206 0.002901       0.003705
## B             0.3188 10.7919 0.038155       0.640319
## 
## 2. Quantiles for each variable:
## 
##                   2.5%     25%      50%      75%  97.5%
## d.FFG.CT_Nav -11.90919 -3.6242 -1.75412  0.08546 10.781
## d.FFG.RA     -15.27536 -4.5789 -2.70554 -0.86358  7.487
## sd.d           0.08883  0.6984  1.24518  1.94588  3.041
## B            -19.71105 -2.8999  0.06112  3.14163 24.948
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.169900 3.195648 7.365548 
## 
## 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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 73.789 116.01  0.41015         8.0951
## d.FFG.RA     -3.766  79.00  0.27929         5.2062
## sd.d         44.209  17.67  0.06248         0.0633
## B            11.553 184.97  0.65396        13.2125
## 
## 2. Quantiles for each variable:
## 
##                 2.5%    25%    50%    75%  97.5%
## d.FFG.CT_Nav -157.11  16.89 70.381 124.86 350.69
## d.FFG.RA     -191.08 -40.20 -1.391  36.49 154.13
## sd.d           18.13  30.23 40.976  56.10  82.77
## B            -357.79 -69.30  5.240  82.99 469.62
## 
## -- Model fit (residual deviance):
## 
##       Dbar         pD        DIC 
##  4.9840864 -0.5360103  4.4480761 
## 
## 5 data points, ratio 0.9968, 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.6

Null Models

map(md.nmr.final$null.mcmc2, ~summary(.x))
## $Blood_Loss
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.1 69.49  0.10030        0.09982
## d.FFG.RA     -119.6 70.06  0.10112        0.10112
## sd.d           88.0 39.45  0.05695        0.05774
## 
## 2. Quantiles for each variable:
## 
##                 2.5%     25%    50%    75%  97.5%
## d.FFG.CT_Nav -279.19 -172.89 -136.1 -97.70  12.79
## d.FFG.RA     -263.76 -157.99 -121.4 -81.78  30.39
## sd.d           28.08   56.12   81.8 116.54 168.88
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.118020 3.971143 8.089163 
## 
## 4 data points, ratio 1.03, I^2 = 27%
## 
## 
## $LOS
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.824 1.2285 0.001773       0.001834
## d.FFG.RA     -2.661 1.2375 0.001786       0.001801
## sd.d          1.352 0.8204 0.001184       0.001542
## 
## 2. Quantiles for each variable:
## 
##                  2.5%    25%    50%    75%   97.5%
## d.FFG.CT_Nav -4.59738 -2.454 -1.700 -1.133  0.5565
## d.FFG.RA     -5.26327 -3.316 -2.666 -2.011 -0.0406
## sd.d          0.08937  0.696  1.238  1.940  3.0407
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.164142 3.607736 7.771879 
## 
## 4 data points, ratio 1.041, I^2 = 28%
## 
## 
## $OP_Time
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.6132 34.01  0.04908        0.04912
## d.FFG.RA      0.9715 28.38  0.04096        0.04108
## sd.d         44.2341 17.66  0.02550        0.02560
## 
## 2. Quantiles for each variable:
## 
##                 2.5%    25%    50%   75%  97.5%
## d.FFG.CT_Nav  -3.548  47.57 66.488 85.65 137.39
## d.FFG.RA     -58.408 -15.17  1.486 17.31  58.52
## sd.d          18.117  30.25 41.044 56.15  82.79
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.985538 4.846598 9.832136 
## 
## 5 data points, ratio 0.9971, 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
## NULL

Randomization = 0

map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))

## $ODI
## NULL
## 
## $VAS_Back
## NULL
## 
## $VAS_Leg
## NULL

MD Outcomes

Randomization = 1 Models

map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))

## $Blood_Loss
## NULL
## 
## $LOS
## NULL
## 
## $OP_Time
## NULL

Randomization = 0

map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))

## $Blood_Loss
## NULL
## 
## $LOS
## NULL
## 
## $OP_Time
## NULL

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.null

Summaries for Moderator Analyses

SMD Outcomes

Covariate Models

map(smd.nmr.final$mcmc3, ~summary(.x))
## $ODI
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.3272 1.0255 0.0036256      0.0227220
## d.FFG.RA     -0.2095 0.5339 0.0018876      0.0115603
## sd.d          0.1863 0.1122 0.0003968      0.0005162
## B             0.0192 1.4417 0.0050972      0.0339714
## 
## 2. Quantiles for each variable:
## 
##                   2.5%     25%        50%     75%  97.5%
## d.FFG.CT_Nav -2.280603 -0.7354 -3.383e-01 0.05619 1.6762
## d.FFG.RA     -1.250426 -0.4312 -2.004e-01 0.02461 0.7984
## sd.d          0.009231  0.0890  1.811e-01 0.28091 0.3818
## B            -2.815142 -0.3637  2.836e-05 0.36134 2.9223
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.727472 2.366326 5.093798 
## 
## 3 data points, ratio 0.9092, 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.47775 2.1995 0.0077765      0.0910089
## d.FFG.RA     -0.35788 0.7581 0.0026802      0.0301643
## sd.d          0.28157 0.1697 0.0005999      0.0008392
## B             0.02045 2.8780 0.0101754      0.1209809
## 
## 2. Quantiles for each variable:
## 
##                  2.5%     25%       50%     75%  97.5%
## d.FFG.CT_Nav -5.08506 -1.0955 -0.496884  0.1035 3.8816
## d.FFG.RA     -1.83000 -0.6010 -0.349015 -0.1035 1.1960
## sd.d          0.01473  0.1401  0.268934  0.4153 0.5995
## B            -6.06251 -0.6044 -0.001689  0.5991 5.7147
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 3.984558 3.171019 7.155577 
## 
## 4 data points, ratio 0.9961, I^2 = 25%
## 
## -- 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.13560 0.61362 0.0021695      0.0116643
## d.FFG.RA     -0.07629 0.32436 0.0011468      0.0058630
## sd.d          0.09872 0.05844 0.0002066      0.0002703
## B             0.02558 0.79685 0.0028173      0.0170848
## 
## 2. Quantiles for each variable:
## 
##                   2.5%      25%      50%    75%  97.5%
## d.FFG.CT_Nav -1.239307 -0.43546 -0.14733 0.1426 1.0900
## d.FFG.RA     -0.715765 -0.23130 -0.07033 0.0911 0.5090
## sd.d          0.004609  0.04816  0.09714 0.1486 0.1978
## B            -1.455115 -0.18628  0.00232 0.1965 1.7281
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.516263 2.159614 4.675877 
## 
## 3 data points, ratio 0.8388, I^2 = 21%
## 
## -- 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

Null Models

map(smd.nmr.final$null.mcmc2, ~summary(.x))
## $ODI
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.3397 0.3710 0.0005355      0.0006318
## d.FFG.RA     -0.2030 0.2362 0.0003410      0.0003806
## sd.d          0.1862 0.1122 0.0001620      0.0002131
## 
## 2. Quantiles for each variable:
## 
##                   2.5%      25%     50%      75%  97.5%
## d.FFG.CT_Nav -1.072461 -0.58294 -0.3396 -0.09698 0.3952
## d.FFG.RA     -0.679161 -0.35371 -0.2013 -0.05059 0.2631
## sd.d          0.008902  0.08907  0.1806  0.28037 0.3821
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.720692 2.366518 5.087209 
## 
## 3 data points, ratio 0.9069, I^2 = 26%
## 
## 
## $VAS_Back
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.4932 0.4493 0.0006486      0.0007031
## d.FFG.RA     -0.3524 0.2427 0.0003503      0.0003699
## sd.d          0.2843 0.1698 0.0002451      0.0003381
## 
## 2. Quantiles for each variable:
## 
##                  2.5%     25%     50%     75%  97.5%
## d.FFG.CT_Nav -1.39260 -0.7735 -0.4933 -0.2137 0.4124
## d.FFG.RA     -0.85596 -0.4975 -0.3487 -0.2037 0.1296
## sd.d          0.01527  0.1426  0.2716  0.4185 0.6009
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 3.979676 3.190921 7.170597 
## 
## 4 data points, ratio 0.9949, I^2 = 25%
## 
## 
## $VAS_Leg
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.15295 0.31875 4.601e-04      0.0007049
## d.FFG.RA     -0.06904 0.18775 2.710e-04      0.0003687
## sd.d          0.09901 0.05824 8.407e-05      0.0001104
## 
## 2. Quantiles for each variable:
## 
##                   2.5%      25%      50%     75%  97.5%
## d.FFG.CT_Nav -0.778450 -0.36686 -0.15279 0.06107 0.4725
## d.FFG.RA     -0.439601 -0.19404 -0.06880 0.05649 0.2996
## sd.d          0.004942  0.04858  0.09765 0.14873 0.1977
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 2.506104 2.157577 4.663681 
## 
## 3 data points, ratio 0.8354, I^2 = 20%

MD Outcomes

Covariate Models

map(md.nmr.final$mcmc3, ~summary(.x))
## $Blood_Loss
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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 -123.35 211.52   0.7478        15.7230
## d.FFG.RA     -131.86 211.65   0.7483        15.6667
## sd.d           87.94  39.57   0.1399         0.1415
## B              23.80 400.00   1.4142        30.4874
## 
## 2. Quantiles for each variable:
## 
##                 2.5%     25%      50%    75% 97.5%
## d.FFG.CT_Nav -548.69 -230.25 -130.388 -29.56 344.2
## d.FFG.RA     -600.24 -225.96 -125.634 -24.48 292.0
## sd.d           27.96   55.96   81.642 116.73 168.8
## B            -798.87 -151.88    9.121 170.67 921.5
## 
## -- Model fit (residual deviance):
## 
##      Dbar        pD       DIC 
##  4.118089 -1.425021  2.693068 
## 
## 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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.6685  5.5340 0.019566       0.311619
## d.FFG.RA     -2.8177  5.5304 0.019553       0.309232
## sd.d          1.3549  0.8206 0.002901       0.003705
## B             0.3188 10.7919 0.038155       0.640319
## 
## 2. Quantiles for each variable:
## 
##                   2.5%     25%      50%      75%  97.5%
## d.FFG.CT_Nav -11.90919 -3.6242 -1.75412  0.08546 10.781
## d.FFG.RA     -15.27536 -4.5789 -2.70554 -0.86358  7.487
## sd.d           0.08883  0.6984  1.24518  1.94588  3.041
## B            -19.71105 -2.8999  0.06112  3.14163 24.948
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.169900 3.195648 7.365548 
## 
## 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:205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 20000 
## 
## 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 73.789 116.01  0.41015         8.0951
## d.FFG.RA     -3.766  79.00  0.27929         5.2062
## sd.d         44.209  17.67  0.06248         0.0633
## B            11.553 184.97  0.65396        13.2125
## 
## 2. Quantiles for each variable:
## 
##                 2.5%    25%    50%    75%  97.5%
## d.FFG.CT_Nav -157.11  16.89 70.381 124.86 350.69
## d.FFG.RA     -191.08 -40.20 -1.391  36.49 154.13
## sd.d           18.13  30.23 40.976  56.10  82.77
## B            -357.79 -69.30  5.240  82.99 469.62
## 
## -- Model fit (residual deviance):
## 
##       Dbar         pD        DIC 
##  4.9840864 -0.5360103  4.4480761 
## 
## 5 data points, ratio 0.9968, 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.6

Null Models

map(md.nmr.final$null.mcmc2, ~summary(.x))
## $Blood_Loss
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.1 69.49  0.10030        0.09982
## d.FFG.RA     -119.6 70.06  0.10112        0.10112
## sd.d           88.0 39.45  0.05695        0.05774
## 
## 2. Quantiles for each variable:
## 
##                 2.5%     25%    50%    75%  97.5%
## d.FFG.CT_Nav -279.19 -172.89 -136.1 -97.70  12.79
## d.FFG.RA     -263.76 -157.99 -121.4 -81.78  30.39
## sd.d           28.08   56.12   81.8 116.54 168.88
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.118020 3.971143 8.089163 
## 
## 4 data points, ratio 1.03, I^2 = 27%
## 
## 
## $LOS
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.824 1.2285 0.001773       0.001834
## d.FFG.RA     -2.661 1.2375 0.001786       0.001801
## sd.d          1.352 0.8204 0.001184       0.001542
## 
## 2. Quantiles for each variable:
## 
##                  2.5%    25%    50%    75%   97.5%
## d.FFG.CT_Nav -4.59738 -2.454 -1.700 -1.133  0.5565
## d.FFG.RA     -5.26327 -3.316 -2.666 -2.011 -0.0406
## sd.d          0.08937  0.696  1.238  1.940  3.0407
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.164142 3.607736 7.771879 
## 
## 4 data points, ratio 1.041, I^2 = 28%
## 
## 
## $OP_Time
## 
## Results on the Mean Difference scale
## 
## Iterations = 5010:1205000
## Thinning interval = 10 
## Number of chains = 4 
## Sample size per chain = 120000 
## 
## 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.6132 34.01  0.04908        0.04912
## d.FFG.RA      0.9715 28.38  0.04096        0.04108
## sd.d         44.2341 17.66  0.02550        0.02560
## 
## 2. Quantiles for each variable:
## 
##                 2.5%    25%    50%   75%  97.5%
## d.FFG.CT_Nav  -3.548  47.57 66.488 85.65 137.39
## d.FFG.RA     -58.408 -15.17  1.486 17.31  58.52
## sd.d          18.117  30.25 41.044 56.15  82.79
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 4.985538 4.846598 9.832136 
## 
## 5 data points, ratio 0.9971, 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
## NULL

Randomization = 0

map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))

## $ODI
## NULL
## 
## $VAS_Back
## NULL
## 
## $VAS_Leg
## NULL

MD Outcomes

Randomization = 1 Models

map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))

## $Blood_Loss
## NULL
## 
## $LOS
## NULL
## 
## $OP_Time
## NULL

Randomization = 0

map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))

## $Blood_Loss
## NULL
## 
## $LOS
## NULL
## 
## $OP_Time
## NULL

Function 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`.