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")
# 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)
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))
)
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))
)
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)
}
netgraph(res$net[[1]], labels = long.labels)
create_dt(tab.list[[1]])
create_dt(tab.list[[2]])
create_dt(tab.list[[3]])
create_dt(md.list[[1]])
create_dt(md.list[[2]])
create_dt(md.list[[3]])
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'")
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)
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
create_dt(md.list[[1]])
create_dt(md.list[[2]])
create_dt(md.list[[3]])
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 --
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 --
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
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
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
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
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
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
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.
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)
}
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
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
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
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
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)
create_dt(smd.dic.df)
create_dt(md.dic.df)
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
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%
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
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%
map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))
## $ODI
## NULL
##
## $VAS_Back
## NULL
##
## $VAS_Leg
## NULL
map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))
## $ODI
## NULL
##
## $VAS_Back
## NULL
##
## $VAS_Leg
## NULL
map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))
## $Blood_Loss
## NULL
##
## $LOS
## NULL
##
## $OP_Time
## NULL
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
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
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%
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
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%
map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))
## $ODI
## NULL
##
## $VAS_Back
## NULL
##
## $VAS_Leg
## NULL
map(smd.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 0)))
## $ODI
## NULL
##
## $VAS_Back
## NULL
##
## $VAS_Leg
## NULL
map(md.nmr.final$mcmc3, ~forest(relative.effect(.x, t1 = "FFG", covariate = 1)))
## $Blood_Loss
## NULL
##
## $LOS
## NULL
##
## $OP_Time
## NULL
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)
}
smd.Final.Tables.null <- tables.func(smd.list.null, "NULL")
smd.Final.Tables.cov <- tables.func(smd.list.cov, "COV")
create_dt(smd.Final.Tables.null$ODI)
create_dt(smd.Final.Tables.cov$ODI)
create_dt(smd.Final.Tables.null$VAS_Back)
create_dt(smd.Final.Tables.cov$VAS_Back)
create_dt(smd.Final.Tables.null$VAS_Leg)
create_dt(smd.Final.Tables.cov$VAS_Leg)
md.Final.Tables.null <- tables.func(md.list.null, "NULL")
md.Final.Tables.cov <- tables.func(md.list.cov, "COV")
create_dt(md.Final.Tables.null$Blood_Loss)
create_dt(md.Final.Tables.cov$Blood_Loss)
create_dt(md.Final.Tables.null$LOS)
create_dt(md.Final.Tables.cov$LOS)
create_dt(md.Final.Tables.null$OP_Time)
create_dt(md.Final.Tables.cov$OP_Time)
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.
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`.
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`.
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`.
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`.