Warning: package 'cmdstanr' was built under R version 4.4.3Simulation results 5
Sequential design with early stopping (restricted action set) - risk based
Load simulation results
# Each input file corresponds to the results from a single simulation
# scenario/configuration.
# Load all the files into a single list.
# files of interest
sim_lab <- "sim05-15"
flist <- list.files(paste0("data/", sim_lab), pattern = "sim05")
toks <- list()
l <- list()
i <- 1
for(i in 1:length(flist)){
l[[i]] <- qs::qread(file.path(paste0("data/", sim_lab), flist[i]))
toks[[i]] <- unlist(tstrsplit(flist[i], "[-.]"))
}Introduction
In the present design, the model is used to compute unit level risk (probability) and assesses decisions based on risk difference for the intervention comparisons by domain. Log-odds scale had limited interpretability for clinical users and is inconsistent in absolute terms across strata with varying baseline rates. It is also useful to explore what an absolute perspective on effectiveness translates to in terms of operating characteristics. Simulation parameters continue to be expressed as log-odds-ratios. This makes the simulation process simpler. Results are presented on both the odds and risk scale. Aim is to give us a better intuition and transparency on the magnitude of effects we are contemplating and determine if these are reasonable assumptions.
We have a single multivariable logistic regression model. One of the sets of parameters accounts for clinician preference for revision type in order to achieve conditional exchangeability across groups. At each interim, we assess the posterior and if a decision threshold is met, we act. For example, if a superiority decision is reached in one of the domains for which this decision type is relevant, then we consider that domain dealt with and all subsequent participants are assigned to receive the superior intervention. We continue to update the posterior for the stopped comparisons in subsequent interim analyses until we get to the point where all questions have been answered in all domains. Once all questions have been answered, the trial recruitment will stop.
The priors were as follows:
- Reference log-odds of response: logistic distribution, mean 0.7 and scale 0.7
- Silo effects: normal distribution, mean 0 and scale 1
- Joint effects: normal distribution, mean 0 and scale 1
- Preference effects: normal distribution, mean 0 and scale 1
- Treatment effects: normal distribution, mean 0 and scale 1
Superiority and non-inferiority are applicable to some domains and not others, however, we define reference and threshold values for all domains, just in case.
For the superiority decision, a reference value of 0 was used and the probability thresholds were:
- 0.92 for surgical domain
- 0.95 for antibiotic duration domain
- 0.95 for extended prophylaxis domain
- 0.995 for antibiotic choice domain
For the futility decision (in relation to superiority) a reference value of 0.05 was used and the probability thresholds were:
- 0.3 for surgical domain
- 0.25 for antibiotic duration domain
- 0.25 for extended prophylaxis domain
- 0.25 for antibiotic choice domain
The above means, for example, that if the probability that the risk difference is greater than 0.05 is less than 0.3 in the surgical domain comparison, then we say the superiority goal is futile.
For the ni decision, a reference value of -0.05 was used and the probability thresholds were:
- 0.975 for surgical domain
- 0.925 for antibiotic duration domain
- 0.975 for extended prophylaxis domain
- 0.975 for antibiotic choice domain
The above means, for example, that if the probability that the risk difference is greater than -0.05 is greater than 0.925, in the antibiotic duration domain, then we will say the intervention is non-inferior.
The futility decision (in relation to non-inferiority) has a reference value of 0 and the probability thresholds were:
- 0.25 for surgical domain
- 0.1 for antibiotic duration domain
- 0.25 for extended prophylaxis domain
- 0.25 for antibiotic choice domain
This means, for example, that if the probability that the risk difference is greater than 0 is less than 0.1, in the antibiotic duration domain, then we say the non-inferiority goal is futile.
Figure 1 attempts to put the superiority rules into pictures based on possible scenarios for assessment of the posterior risk difference for an arbitrary domain where superiority is being assessed. The approach, reference values and thresholds apply to all domains where superiority is assessed.
Analogously, Figure 2 puts the non-inferiority rule into pictures based on possible scenarios for assessment of the posterior risk difference for an arbitrary domain where non-inferiority is being assessed. The approach, reference values and thresholds apply to all domains.
For this set of simulations, the number of simulated trials per scenario was 1000, the simulation label is sim05-15.
Simulation results
Cumulative probability of each decision type
# Cumulative probability of decisions:
# Traverse the list of simulation results and for each one summarise the
# cumulative probability of each decision type.
i <- 1
d_tbl_1 <- data.table()
# For each scenario that was simulated
for(i in 1:length(l)){
# extract the decision matrix - sim, analysis, quantity, domain level decision
d_dec_1 <- copy(l[[i]]$d_decision)
# config for scenario
l_cfg <- copy(l[[i]]$cfg)
# number of enrolments at each interim (interim sample size sequence)
d_N <- data.table(analys = seq_along(l_cfg$N_pt), N = l_cfg$N_pt)
# long version of decisions
d_dec_2 <- melt(d_dec_1,
id.vars = c("sim", "analys", "quant"),
variable.name = "domain")
# Should be right, but just in case...
if(any(is.na(d_dec_2$value))){
message("Some of the decision values are NA in index ", i, " file ", flist[i])
d_dec_2[is.na(value), value := FALSE]
}
d_dec_2[, domain := as.numeric(gsub("d", "", domain))]
# Domains 1, 3 and 4 will stop for superiority or futility for superiority.
# Domaain 2 will stop for NI or futility for NI.
# No other stopping rules apply and so we only evaluate the operating
# characteristics on these, i.e. we do not care about the results for the
# cumualative probability of ni for domain 1, 3 and 4 because we would never
# stop for this.
d_dec_2 <- rbind(
d_dec_2[domain %in% c(1, 3, 4) & quant %in% c("sup", "fut_sup")],
d_dec_2[domain %in% c(2) & quant %in% c("ni", "fut_ni")]
)
# compute the cumulative instances of a decision being made by sim, each
# decision type and by parameter
d_dec_2[, value := as.logical(cumsum(value)>0), keyby = .(sim, quant, domain)]
d_dec_2 <- merge(d_dec_2, d_N, by = "analys")
# cumulative proportion for which each decision quantity has been met by
# analysis and domain
d_dec_cumprob <- d_dec_2[, .(pr_val = mean(value)), keyby = .(analys, N, quant, domain)]
d_tbl_1 <- rbind(
d_tbl_1,
cbind(scenario = i, desc = l_cfg$desc, d_dec_cumprob)
)
}Expected sample size
# Expected sample size
# Similar to above but focus on expected sample size
# Traverse the list of simulation results and for each one summarise the
# sample sizes at which stopping for a domain occurs for any reason.
# All we are trying to get to is the expected sample size by domain and
# are not really interested in what decision was made. The cumulative prob
# of each decision type is computed previously.
i <- 2
d_tbl_2 <- data.table()
for(i in 1:length(l)){
# extract the decision matrix - sim, analysis, quantity, domain level decision
d_dec_1 <- copy(l[[i]]$d_decision)
# config for scenario
l_cfg <- copy(l[[i]]$cfg)
# interim looks
d_N <- data.table(analys = seq_along(l_cfg$N_pt), N = l_cfg$N_pt)
# long version of decision
d_dec_2 <- melt(d_dec_1,
id.vars = c("sim", "analys", "quant"),
variable.name = "domain")
# Should be right, but just in case...
if(any(is.na(d_dec_2$value))){
message("Some of the decision values are NA in index ", i, " file ", flist[i])
d_dec_2[is.na(value), value := FALSE]
}
d_dec_2[, domain := as.numeric(gsub("d", "", domain))]
d_dec_2 <- merge(d_dec_2, d_N, by = "analys")
# First instance of any decision rule being hit by sim and domain.
# Domains 1, 3 and 4 will stop for superiority or futility for superiority.
# Domaain 2 will stop for NI or futility for NI.
# Sometimes a decision rule will not be hit for a domain and it will continue
# to the max sample size. We will deal with that in a minute.
d_dec_stop <- rbind(
d_dec_2[
domain %in% c(1, 3, 4) & value == T & quant %in% c("sup", "fut_sup"),
.SD[1], keyby = .(sim, domain)],
d_dec_2[
domain %in% c(2) & value == T & quant %in% c("ni", "fut_ni"),
.SD[1], keyby = .(sim, domain)]
)
setnames(d_dec_stop, "N", "N_stopped")
setnames(d_dec_stop, "value", "stopped_early")
# Add in any rows for which no early stopping happened
d_dec_stop <- merge(d_dec_stop[, .SD, .SDcols = !c("analys")],
# all combinations of sim and domain
unique(d_dec_2[, .(sim, domain)]),
by = c("sim", "domain"), all.y = T)
# If domain or trial not stopped then record as having run to the
# maximum sample size with no decision made.
d_dec_stop[is.na(N_stopped), N_stopped := max(d_N$N)]
d_dec_stop[is.na(stopped_early), stopped_early := F]
# So now we know where every domain was stopped and the reason for stopping
# within each sim. Great.
d_dec_stop[is.na(quant), quant := "null"]
d_tbl_2 <- rbind(
d_tbl_2,
cbind(
scenario = i,
desc = l_cfg$desc,
d_dec_stop
)
)
}Distributions of posterior means (unconditional)
# Distribution of posterior means for parameters of interest.
# Some simulated trials will have stopped prior to the maximum sample size and
# these will have NA for their posterior means. If you were to summarise these
# posterior means, they would be conditional on the trial having 'survived'
# until the relevant interim. This means that you have missing data at later
# interims, which creates a selection bias in that your selection of sims at any
# given interim are not a random sample, but rather a sample conditioned on the
# stopping rules.
# If you do not account for this in some way then a summary can be either
# optimistic or pessimistic depending on how the stopping rules interact
# with the data. Here we try to account for this missingness by imputing the
# missing posterior means with locf within each simulation.
# Note that this is really only a partial 'fix' to get a sense of whether
# our estimates is representative of the parameter values we used to simulate
# the data.
i <- 2
d_tbl_3 <- data.table()
for(i in 1:length(l)){
# config for scenario
l_cfg <- copy(l[[i]]$cfg)
# params
d_pars <- copy(l[[i]]$d_post_smry_1)
d_pars <- d_pars[par %in% c("lor", "rd")]
# interim looks
d_N <- data.table(analys = seq_along(l_cfg$N_pt), N = l_cfg$N_pt)
d_pars <- dcast(d_pars, sim + id_analys + domain ~ par, value.var = c("mu", "se"))
# locf
d_pars[, `:=`(mu_lor = nafill(mu_lor, type = "locf"),
mu_rd = nafill(mu_rd, type = "locf"),
se_lor = nafill(se_lor, type = "locf"),
se_rd = nafill(se_rd, type = "locf")
),
keyby = .(sim, domain)]
setnames(d_pars, "id_analys", "analys")
#
d_pars <- merge(d_pars, d_N, by = "analys")
d_tbl_3 <- rbind(
d_tbl_3,
cbind(
scenario = i, desc = l_cfg$desc,
d_pars[, .(analys, sim, domain, N, mu_lor, mu_rd, se_lor, se_rd)]
)
)
}Distributions of posterior means (conditional)
# Distribution of posterior means for parameters of interest.
# By way of contrast to the above, here we summarise based on the stopping
# rules. That is, the posterior means that are included in the summary are only
# those that 'survived' until the relevant interim. This is simply to give you
# a sense of how extreme the selection effects associated with stopping can be.
i <- 6
d_tbl_4 <- data.table()
for(i in 1:length(l)){
# extract the decision matrix - sim, analysis, quantity, domain level decision
d_dec_1 <- copy(l[[i]]$d_decision)
# config for scenario
l_cfg <- copy(l[[i]]$cfg)
bd1 <- unlist(l_cfg$d1)
bd2 <- unlist(l_cfg$d2)
bd3 <- unlist(l_cfg$d3)
bd4 <- unlist(l_cfg$d4)
d_tru <- copy(l[[i]]$d_all)
setnames(d_tru, "id_analys", "analys")
d_tru <- dcast(
d_tru[, .(N = sum(N)), keyby = .(sim, analys, pref_rev)],
sim + analys ~ pref_rev, value.var = "N"
)
colnames(d_tru) <- c("sim", "analys", "N1", "N2")
d_tru[, `:=`(prop_p1 = N1 / (N1 + N2),
prop_p2 = N2 / (N1 + N2))]
d_tru[, db1 := as.numeric(
(d_tru$prop_p1 * bd1[2] + d_tru$prop_p2 * bd1[3]) - bd1[1])]
d_tru[, db2 := as.numeric(bd2[3] - bd2[2])]
d_tru[, db3 := as.numeric(bd3[3] - bd3[2])]
d_tru[, db4 := as.numeric(bd4[3] - bd4[2])]
d_tru[, `:=`(
N1 = NULL, N2 = NULL, prop_p1 = NULL, prop_p2 = NULL
)]
d_tru <- melt(
d_tru,
measure.vars = c("db1", "db2", "db3", "db4"),
value.name = "lor_tru", variable.name = "domain"
)
d_tru[, domain := as.numeric(gsub("db", "", domain))]
# params
d_pars <- copy(l[[i]]$d_post_smry_1)
d_pars <- d_pars[par %in% c("lor", "rd")]
# interim looks
d_N <- data.table(analys = seq_along(l_cfg$N_pt), N = l_cfg$N_pt)
# long version of decision - assessment of each quantity by sim, domain and analysis
d_dec_2 <- melt(d_dec_1,
id.vars = c("sim", "analys", "quant"),
variable.name = "domain")
# Should be right, but just in case...
if(any(is.na(d_dec_2$value))){
message("Some of the decision values are NA in index ", i, " file ", flist[i])
d_dec_2[is.na(value), value := FALSE]
}
d_dec_2[, domain := as.numeric(gsub("d", "", domain))]
d_dec_2 <- merge(d_dec_2, d_N, by = "analys")
# First instance of any decision rule being hit by sim and domain.
# Domains 1, 3 and 4 will stop for superiority or futility for superiority.
# Domaain 2 will stop for NI or futility for NI.
# Sometimes a decision rule will not be hit for a domain and it will continue
# to the max sample size. We will deal with that in a minute.
d_dec_stop <- rbind(
# value being set to TRUE indicates we hit a stopping run
d_dec_2[
domain %in% c(1, 3, 4) & value == T & quant %in% c("sup", "fut_sup"),
.SD[1], keyby = .(sim, domain)],
d_dec_2[
domain %in% c(2) & value == T & quant %in% c("ni", "fut_ni"),
.SD[1], keyby = .(sim, domain)]
)
setnames(d_dec_stop, "N", "N_stopped")
setnames(d_dec_stop, "value", "stopped_early")
# Add in any the domains where we did not stop early
d_dec_stop <- merge(d_dec_stop[, .SD, .SDcols = !c("analys")],
# all combinations of sim and domain
unique(d_dec_2[, .(sim, domain)]),
by = c("sim", "domain"), all.y = T)
# If domain or trial not stopped then record as having run to the
# maximum sample size with no decision made.
d_dec_stop[is.na(N_stopped), N_stopped := max(d_N$N)]
# get rid of quant
d_dec_stop[is.na(stopped_early), stopped_early := F]
# So now we know where every domain was stopped for any reason
# within each sim. Great.
# Add in the other looks at the data so that we have a complete set
d_dec_stop <- merge(
d_dec_stop,
unique(d_dec_2[, .(analys, sim, domain, N)]),
by = c("sim", "domain"),
all.y = T
)
# Merge in the parameter estimates (posterior means)
d_dec_stop <- merge(
d_dec_stop,
dcast(d_pars, sim + id_analys + domain ~ par, value.var = "mu"),
by.x = c("analys", "sim", "domain"),
by.y = c("id_analys", "sim", "domain")
)
# so now drop any records for which we had stopped the trial
d_dec_stop <- d_dec_stop[!is.na(lor)]
d_dec_stop <- merge(
d_dec_stop,
d_tru,
by = c("sim", "analys", "domain")
)
# d_dec_stop[quant == "sup", mean(type_m_ratio), keyby = .(analys, domain, N)]
# out of all the decision that were made, i.e. those where quant is not NA
# what proportion
d_tbl_4 <- rbind(
d_tbl_4,
cbind(
scenario = i, desc = l_cfg$desc,
d_dec_stop[, .(sim, analys, domain, quant, N, lor, rd, stopped_early, N_stopped)]
)
)
}Number of participants for each randomised comparison over time
#
i <- 1
d_N_by_arm <- data.table()
for(i in 1:length(l)){
# extract the decision matrix - sim, analysis, quantity, domain level decision
d_dec_1 <- copy(l[[i]]$d_decision)
# long version of decisions
d_dec_2 <- melt(
d_dec_1, id.vars = c("sim", "analys", "quant"), variable.name = "domain")
# Should be right, but just in case...
if(any(is.na(d_dec_2$value))){
message("Some of the decision values are NA in index ", i, " file ", flist[i])
d_dec_2[is.na(value), value := FALSE]
}
d_dec_2[, domain := as.numeric(gsub("d", "", domain))]
# Domains 1, 3 and 4 will stop for superiority or futility for superiority.
# Domaain 2 will stop for NI or futility for NI.
# No other stopping rules apply and so we only evaluate the operating
# characteristics on these, i.e. we do not care about the results for the
# cumualative probability of ni for domain 1, 3 and 4 because we would never
# stop for this.
d_dec_2 <- rbind(
d_dec_2[domain %in% c(1, 3, 4) & quant %in% c("sup", "fut_sup")],
d_dec_2[domain %in% c(2) & quant %in% c("ni", "fut_ni")]
)
# config for scenario
l_cfg <- copy(l[[i]]$cfg)
# interim looks
d_enrolment <- data.table(analys = seq_along(l_cfg$N_pt), N_enrol = l_cfg$N_pt)
# observed trial data sets
d_tru <- copy(l[[i]]$d_all)
setnames(d_tru, "id_analys", "analys")
# late acute, surgical arms
# *** silo isn't included in the keyby because we want to average across
# the entire trial population, not silo specific sample sizes.
d_1_tmp <- d_tru[
silo == 2, .(N = sum(N), domain = 1), keyby = .(sim, analys, d1)]
d_1_tmp[, N := cumsum(N), keyby = .(sim, d1)]
setnames(d_1_tmp, old = "d1", "arm")
d_1_tmp <- merge(
d_1_tmp,
CJ(sim = 1:max(d_tru$sim), analys = 1:5, arm = 1:3, domain = 1),
by = c("sim", "analys", "arm", "domain"),
all.y = T)
d_1_tmp[, N := nafill(N, type = "locf"), keyby = .(sim, arm)]
setkey(d_1_tmp, sim, analys, arm)
# one-stage revision, ab duration randomised arms
d_2_tmp <- d_tru[
d1 == 2 & d2 %in% 2:3,
.(N = sum(N), domain = 2), keyby = .(sim, analys, d2)]
d_2_tmp[, N := cumsum(N), keyby = .(sim, d2)]
setnames(d_2_tmp, old = "d2", "arm")
d_2_tmp <- merge(
d_2_tmp,
CJ(sim = 1:max(d_tru$sim), analys = 1:5, arm = 2:3, domain = 2),
by = c("sim", "analys", "arm", "domain"),
all.y = T)
d_2_tmp[analys == 1 & is.na(N), N := 0]
d_2_tmp[, N := nafill(N, type = "locf"), keyby = .(sim, arm)]
setkey(d_2_tmp, sim, analys, arm)
# two-stage revision, extended prophy randomised arms
d_3_tmp <- d_tru[
d1 == 3 & d3 %in% 2:3,
.(N = sum(N), domain = 3), keyby = .(sim, analys, d3)]
d_3_tmp[, N := cumsum(N), keyby = .(sim, d3)]
setnames(d_3_tmp, old = "d3", "arm")
d_3_tmp <- merge(
d_3_tmp,
CJ(sim = 1:max(d_tru$sim), analys = 1:5, arm = 2:3, domain = 3),
by = c("sim", "analys", "arm", "domain"),
all.y = T)
d_3_tmp[analys == 1 & is.na(N), N := 0]
d_3_tmp[, N := nafill(N, type = "locf"), keyby = .(sim, arm)]
setkey(d_3_tmp, sim, analys, arm)
# ab choice randomised arms
d_4_tmp <- d_tru[
d4 %in% 2:3, .(N = sum(N), domain = 4), keyby = .(sim, analys, d4)]
d_4_tmp[, N := cumsum(N), keyby = .(sim, d4)]
setnames(d_4_tmp, old = "d4", "arm")
d_4_tmp <- merge(
d_4_tmp,
CJ(sim = 1:max(d_tru$sim), analys = 1:5, arm = 2:3, domain = 4),
by = c("sim", "analys", "arm", "domain"),
all.y = T)
d_4_tmp[analys == 1 & is.na(N), N := 0]
d_4_tmp[, N := nafill(N, type = "locf"), keyby = .(sim, arm)]
setkey(d_4_tmp, sim, analys, arm)
d_arm_tmp <- rbind(d_1_tmp, d_2_tmp, d_3_tmp, d_4_tmp)
# d_arm_tmp
# d_dec_2[value == T, .SD[1], keyby = .(sim, analys, domain, quant)]
# dcast(
# d_dec_1[
# sim == 5 & quant %in% c("sup", "fut_sup"), .(sim, analys, quant, d1)],
# sim + analys ~ quant, value.var = "d1"
# )
d_arm_tmp[, `:=`(
scenario = i, desc = l_cfg$desc
)]
d_arm_tmp <- merge(
d_arm_tmp,
d_enrolment,
by = "analys")
d_N_by_arm <- rbind(d_N_by_arm, d_arm_tmp)
}Probability of triggering decision
Table 1 shows the cumulative probability of a superiority decision across each of the scenarios simulated (the same information is shown in Figure 3). Operating characteristics are shown only for the relevant domains and the futility of a superiority decision is included in parentheses.
In the first eight scenarios, the surgical effects are simulated assuming equal sizes in all silos. In the second set of scenarios, silo specific heterogeneity is introduced. Under silo-specific heterogeneity there is inflation of the type-i assertion probability in the surgical domain.
Code
d_tbl_1_cur <- d_tbl_1[quant %in% c("sup", "fut_sup") & domain %in% c(1, 3, 4), .SD]
setorderv(d_tbl_1_cur, cols = c("scenario", "domain", "analys", "quant"),
order = c(1, 1, 1, -1))
d_tbl_1_cur <- dcast(d_tbl_1_cur, scenario + desc + domain ~ quant + N, value.var = "pr_val")
d_tbl_1_cur <- d_tbl_1_cur[, .SD, .SDcols = !c("scenario")]
d_tbl_1_cur[, domain := factor(domain,
levels = c(1, 3, 4),
labels = c("Surgical", "AB Ext-proph", "AB Choice"))]
g_tbl <- d_tbl_1_cur |>
gt(groupname_col = "desc") |>
gt::text_transform(
locations = cells_row_groups(),
fn = function(x) {
lapply(x, function(x) {
gt::md(paste0("*", x, "*"))
})
}
) |>
cols_align(
columns = 1,
align = "left"
) |>
cols_align(
columns = 2:ncol(d_tbl_1_cur),
align = "center"
) |>
cols_merge(
columns = c("sup_500", "fut_sup_500"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("sup_1000", "fut_sup_1000"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("sup_1500", "fut_sup_1500"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("sup_2000", "fut_sup_2000"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("sup_2500", "fut_sup_2500"
),
pattern = "<<{1}>><< ({2})>>"
) |>
tab_spanner(
label = md("Cumulative probability of **superiority** (futility) decision"),
columns = 3:ncol(d_tbl_1_cur)
) |>
cols_label(
domain = "Domain",
sup_500 = html("500"),
sup_1000 = html("1000"),
sup_1500 = html("1500"),
sup_2000 = html("2000"),
sup_2500 = html("2500")
) |>
tab_options(
table.font.size = "70%"
) |>
fmt_number(decimals = 3, drop_trailing_zeros = TRUE)
g_tbl| Domain |
Cumulative probability of superiority (futility) decision
|
||||
|---|---|---|---|---|---|
| 500 | 1000 | 1500 | 2000 | 2500 | |
| Null effect in all domains | |||||
| Surgical | 0.045 (0.521) | 0.075 (0.667) | 0.099 (0.709) | 0.107 (0.737) | 0.108 (0.743) |
| AB Ext-proph | 0.052 (0.437) | 0.081 (0.597) | 0.098 (0.682) | 0.108 (0.729) | 0.118 (0.756) |
| AB Choice | 0.006 (0.596) | 0.011 (0.787) | 0.013 (0.868) | 0.015 (0.916) | 0.017 (0.944) |
| Moderate (OR 1.75) surgical revision effect (one and two-stage) | |||||
| Surgical | 0.568 (0.045) | 0.739 (0.052) | 0.794 (0.054) | 0.818 (0.055) | 0.838 (0.056) |
| AB Ext-proph | 0.052 (0.506) | 0.082 (0.721) | 0.103 (0.801) | 0.11 (0.842) | 0.118 (0.863) |
| AB Choice | 0.01 (0.602) | 0.016 (0.786) | 0.02 (0.878) | 0.023 (0.931) | 0.024 (0.945) |
| Moderate (OR 1.75) surgical revision effect (one-stage) | |||||
| Surgical | 0.203 (0.26) | 0.297 (0.336) | 0.346 (0.359) | 0.373 (0.38) | 0.393 (0.391) |
| AB Ext-proph | 0.05 (0.49) | 0.087 (0.653) | 0.102 (0.734) | 0.113 (0.79) | 0.121 (0.816) |
| AB Choice | 0.003 (0.612) | 0.007 (0.814) | 0.008 (0.898) | 0.01 (0.942) | 0.01 (0.962) |
| Moderate (OR 1.75) surgical revision effect (two-stage) | |||||
| Surgical | 0.154 (0.258) | 0.257 (0.342) | 0.303 (0.378) | 0.34 (0.403) | 0.358 (0.411) |
| AB Ext-proph | 0.048 (0.479) | 0.077 (0.633) | 0.104 (0.728) | 0.115 (0.776) | 0.122 (0.806) |
| AB Choice | 0.005 (0.594) | 0.009 (0.778) | 0.012 (0.868) | 0.015 (0.923) | 0.017 (0.953) |
| Moderate (OR 1.75) antibiotic duration 6wk effect | |||||
| Surgical | 0.052 (0.494) | 0.087 (0.622) | 0.104 (0.667) | 0.111 (0.681) | 0.114 (0.689) |
| AB Ext-proph | 0.043 (0.493) | 0.073 (0.665) | 0.086 (0.722) | 0.096 (0.761) | 0.102 (0.801) |
| AB Choice | 0 (0.557) | 0.002 (0.787) | 0.006 (0.881) | 0.006 (0.921) | 0.007 (0.952) |
| Moderate (OR 1.75) antibiotic ext-proph 12wk effect | |||||
| Surgical | 0.059 (0.501) | 0.092 (0.643) | 0.107 (0.704) | 0.113 (0.721) | 0.12 (0.729) |
| AB Ext-proph | 0.414 (0.056) | 0.676 (0.069) | 0.784 (0.077) | 0.844 (0.078) | 0.869 (0.08) |
| AB Choice | 0.009 (0.58) | 0.011 (0.792) | 0.013 (0.871) | 0.015 (0.923) | 0.015 (0.948) |
| Moderate (OR 1.75) antibiotic choice rifampacin effect | |||||
| Surgical | 0.047 (0.519) | 0.078 (0.657) | 0.091 (0.718) | 0.096 (0.741) | 0.1 (0.75) |
| AB Ext-proph | 0.055 (0.456) | 0.074 (0.608) | 0.086 (0.702) | 0.102 (0.747) | 0.108 (0.781) |
| AB Choice | 0.337 (0.024) | 0.751 (0.028) | 0.896 (0.028) | 0.945 (0.028) | 0.966 (0.028) |
| Moderate (OR 1.75) effects in all domains | |||||
| Surgical | 0.649 (0.028) | 0.801 (0.037) | 0.853 (0.04) | 0.87 (0.041) | 0.881 (0.041) |
| AB Ext-proph | 0.379 (0.09) | 0.669 (0.116) | 0.815 (0.12) | 0.863 (0.123) | 0.88 (0.124) |
| AB Choice | 0.307 (0.044) | 0.642 (0.054) | 0.836 (0.057) | 0.904 (0.057) | 0.927 (0.058) |
| Null effect in all domains +silo specific effects | |||||
| Surgical | 0.115 (0.312) | 0.197 (0.424) | 0.241 (0.471) | 0.26 (0.5) | 0.282 (0.513) |
| AB Ext-proph | 0.041 (0.501) | 0.063 (0.672) | 0.078 (0.746) | 0.092 (0.797) | 0.097 (0.838) |
| AB Choice | 0.005 (0.578) | 0.009 (0.76) | 0.015 (0.854) | 0.021 (0.899) | 0.026 (0.931) |
| Moderate (OR 1.75) surgical revision effect (one and two-stage) +silo specific effects | |||||
| Surgical | 0.616 (0.033) | 0.778 (0.04) | 0.831 (0.042) | 0.854 (0.043) | 0.872 (0.045) |
| AB Ext-proph | 0.05 (0.51) | 0.076 (0.712) | 0.093 (0.802) | 0.108 (0.844) | 0.121 (0.864) |
| AB Choice | 0.009 (0.543) | 0.011 (0.774) | 0.015 (0.868) | 0.017 (0.921) | 0.017 (0.954) |
| Moderate (OR 1.75) surgical revision effect (one-stage) +silo specific effects | |||||
| Surgical | 0.286 (0.139) | 0.428 (0.181) | 0.503 (0.199) | 0.543 (0.213) | 0.566 (0.218) |
| AB Ext-proph | 0.045 (0.479) | 0.071 (0.672) | 0.083 (0.767) | 0.1 (0.819) | 0.113 (0.854) |
| AB Choice | 0.005 (0.606) | 0.006 (0.799) | 0.008 (0.872) | 0.011 (0.924) | 0.012 (0.944) |
| Moderate (OR 1.75) surgical revision effect (two-stage) +silo specific effects | |||||
| Surgical | 0.439 (0.067) | 0.596 (0.084) | 0.655 (0.093) | 0.674 (0.095) | 0.692 (0.097) |
| AB Ext-proph | 0.048 (0.478) | 0.07 (0.676) | 0.09 (0.791) | 0.1 (0.841) | 0.105 (0.864) |
| AB Choice | 0.004 (0.597) | 0.006 (0.8) | 0.008 (0.895) | 0.009 (0.938) | 0.009 (0.962) |
| Moderate (OR 1.75) antibiotic duration 6wk effect +silo specific effects | |||||
| Surgical | 0.17 (0.239) | 0.261 (0.328) | 0.31 (0.361) | 0.334 (0.376) | 0.351 (0.385) |
| AB Ext-proph | 0.06 (0.456) | 0.087 (0.654) | 0.108 (0.737) | 0.121 (0.79) | 0.124 (0.818) |
| AB Choice | 0.006 (0.592) | 0.01 (0.783) | 0.01 (0.874) | 0.011 (0.924) | 0.014 (0.955) |
| Moderate (OR 1.75) antibiotic ext-proph 12wk effect +silo specific effects | |||||
| Surgical | 0.202 (0.221) | 0.312 (0.305) | 0.364 (0.34) | 0.401 (0.36) | 0.425 (0.372) |
| AB Ext-proph | 0.405 (0.066) | 0.661 (0.078) | 0.794 (0.086) | 0.847 (0.089) | 0.881 (0.091) |
| AB Choice | 0.005 (0.603) | 0.01 (0.806) | 0.013 (0.904) | 0.013 (0.94) | 0.015 (0.96) |
| Moderate (OR 1.75) antibiotic choice rifampacin effect +silo specific effects | |||||
| Surgical | 0.148 (0.286) | 0.234 (0.366) | 0.285 (0.401) | 0.312 (0.427) | 0.337 (0.438) |
| AB Ext-proph | 0.043 (0.468) | 0.076 (0.651) | 0.096 (0.754) | 0.111 (0.793) | 0.114 (0.828) |
| AB Choice | 0.339 (0.034) | 0.709 (0.038) | 0.876 (0.038) | 0.943 (0.038) | 0.962 (0.038) |
| Moderate (OR 1.75) effects in all domains +silo specific effects | |||||
| Surgical | 0.653 (0.025) | 0.797 (0.032) | 0.839 (0.033) | 0.86 (0.035) | 0.875 (0.035) |
| AB Ext-proph | 0.367 (0.114) | 0.669 (0.141) | 0.801 (0.148) | 0.836 (0.153) | 0.85 (0.158) |
| AB Choice | 0.322 (0.033) | 0.673 (0.043) | 0.844 (0.045) | 0.919 (0.046) | 0.945 (0.047) |
Code
d_fig <- d_tbl_1[quant %in% c("sup", "fut_sup") & domain %in% c(1, 3, 4), .SD]
setorderv(d_fig, cols = c("scenario", "domain", "analys", "quant"),
order = c(1, 1, 1, -1))
d_fig[, domain := factor(domain,
levels = c(1, 3, 4),
labels = c("Surgical", "AB Ext-proph", "AB Choice"))]
d_fig[, quant := factor(quant,
levels = c("sup", "fut_sup"),
labels = c("Superiority", "Futility"))]
d_fig[, desc := factor(desc,
levels = unique(d_fig$desc),
labels = unique(d_fig$desc))]
ggplot(d_fig, aes(x = N, y = pr_val, group = quant, col = quant)) +
geom_line(lwd = 0.25) +
scale_y_continuous("") +
scale_color_discrete("") +
# facet_grid(desc ~ domain) +
# facet_grid2(desc ~ domain, render_empty = FALSE)
ggh4x::facet_grid2(desc ~ domain ,
labeller = labeller(desc = label_wrap_gen(35)),
scales = "free",
axes = "y",
independent = "y") +
# theme_minimal() +
theme(
legend.position = "bottom",
strip.text.y.right = element_text(angle = 0,
hjust = 0,
vjust = 0.2,
size = 4),
axis.ticks = element_blank(),
strip.text.x.top = element_text(size = 4),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey",
linewidth = 0.1,
linetype = 1),
axis.text.y = element_text(size = 4),
axis.text.x = element_text(size = 4),
axis.title.x = element_text(size = 5),
legend.text = element_text(size = 4)
)
Table 2 shows the cumulative probability of a non-inferiority decision with futility shown in parentheses (the same information is shown in Figure 4). The results are only shown for the domains for which non-inferiority is evaluated.
The “Null effect in all domains” is actually a bit of a misnomer as the true null with regards to the NI decision would be at the NI margin, whereas the scenario refers to the setting where all effects are set to zero. Thus an inflation over the usual type-i assertion probability is to be expected.
Code
d_tbl_1_cur <- d_tbl_1[quant %in% c("ni", "fut_ni") & domain %in% c(2), .SD]
setorderv(d_tbl_1_cur, cols = c("scenario", "domain", "analys", "quant"),
order = c(1, 1, 1, -1))
d_tbl_1_cur <- dcast(d_tbl_1_cur, scenario + desc + domain ~ quant + N, value.var = "pr_val")
d_tbl_1_cur <- d_tbl_1_cur[, .SD, .SDcols = !c("scenario")]
d_tbl_1_cur[, domain := factor(domain,
levels = c(2),
labels = c("AB Duration"))]
g_tbl <- d_tbl_1_cur |>
gt(groupname_col = "desc") |>
gt::text_transform(
locations = cells_row_groups(),
fn = function(x) {
lapply(x, function(x) {
gt::md(paste0("*", x, "*"))
})
}
) |>
cols_align(
columns = 1,
align = "left"
) |>
cols_align(
columns = 2:ncol(d_tbl_1_cur),
align = "center"
) |>
cols_merge(
columns = c("ni_500", "fut_ni_500"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("ni_1000", "fut_ni_1000"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("ni_1500", "fut_ni_1500"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("ni_2000", "fut_ni_2000"
),
pattern = "<<{1}>><< ({2})>>"
) |>
cols_merge(
columns = c("ni_2500", "fut_ni_2500"
),
pattern = "<<{1}>><< ({2})>>"
) |>
tab_spanner(
label = md("Cumulative probability of **NI** (futility) decision"),
columns = 3:ncol(d_tbl_1_cur)
) |>
cols_label(
domain = "Domain",
ni_500 = html("500"),
ni_1000 = html("1000"),
ni_1500 = html("1500"),
ni_2000 = html("2000"),
ni_2500 = html("2500")
) |>
tab_options(
table.font.size = "70%"
) |>
fmt_number(decimals = 3, drop_trailing_zeros = TRUE)
g_tbl| Domain |
Cumulative probability of NI (futility) decision
|
||||
|---|---|---|---|---|---|
| 500 | 1000 | 1500 | 2000 | 2500 | |
| Null effect in all domains | |||||
| AB Duration | 0.149 (0.095) | 0.248 (0.133) | 0.298 (0.161) | 0.342 (0.184) | 0.371 (0.197) |
| Moderate (OR 1.75) surgical revision effect (one and two-stage) | |||||
| AB Duration | 0.159 (0.098) | 0.284 (0.145) | 0.366 (0.185) | 0.44 (0.221) | 0.497 (0.235) |
| Moderate (OR 1.75) surgical revision effect (one-stage) | |||||
| AB Duration | 0.155 (0.097) | 0.256 (0.134) | 0.324 (0.173) | 0.387 (0.197) | 0.438 (0.213) |
| Moderate (OR 1.75) surgical revision effect (two-stage) | |||||
| AB Duration | 0.156 (0.084) | 0.254 (0.142) | 0.33 (0.173) | 0.389 (0.196) | 0.437 (0.215) |
| Moderate (OR 1.75) antibiotic duration 6wk effect | |||||
| AB Duration | 0.451 (0.006) | 0.65 (0.007) | 0.783 (0.008) | 0.855 (0.008) | 0.896 (0.008) |
| Moderate (OR 1.75) antibiotic ext-proph 12wk effect | |||||
| AB Duration | 0.148 (0.096) | 0.226 (0.141) | 0.276 (0.178) | 0.314 (0.199) | 0.351 (0.22) |
| Moderate (OR 1.75) antibiotic choice rifampacin effect | |||||
| AB Duration | 0.134 (0.084) | 0.218 (0.135) | 0.279 (0.167) | 0.34 (0.183) | 0.373 (0.198) |
| Moderate (OR 1.75) effects in all domains | |||||
| AB Duration | 0.436 (0.012) | 0.732 (0.014) | 0.875 (0.015) | 0.928 (0.016) | 0.968 (0.016) |
| Null effect in all domains +silo specific effects | |||||
| AB Duration | 0.139 (0.085) | 0.231 (0.135) | 0.314 (0.156) | 0.355 (0.18) | 0.402 (0.192) |
| Moderate (OR 1.75) surgical revision effect (one and two-stage) +silo specific effects | |||||
| AB Duration | 0.163 (0.075) | 0.299 (0.134) | 0.379 (0.174) | 0.454 (0.197) | 0.509 (0.216) |
| Moderate (OR 1.75) surgical revision effect (one-stage) +silo specific effects | |||||
| AB Duration | 0.159 (0.086) | 0.277 (0.14) | 0.356 (0.169) | 0.422 (0.19) | 0.469 (0.214) |
| Moderate (OR 1.75) surgical revision effect (two-stage) +silo specific effects | |||||
| AB Duration | 0.646 (0.001) | 0.898 (0.001) | 0.98 (0.001) | 0.993 (0.001) | 0.995 (0.001) |
| Moderate (OR 1.75) antibiotic duration 6wk effect +silo specific effects | |||||
| AB Duration | 0.401 (0.021) | 0.668 (0.025) | 0.81 (0.027) | 0.866 (0.028) | 0.9 (0.028) |
| Moderate (OR 1.75) antibiotic ext-proph 12wk effect +silo specific effects | |||||
| AB Duration | 0.145 (0.084) | 0.243 (0.137) | 0.312 (0.171) | 0.364 (0.195) | 0.41 (0.211) |
| Moderate (OR 1.75) antibiotic choice rifampacin effect +silo specific effects | |||||
| AB Duration | 0.138 (0.088) | 0.232 (0.125) | 0.312 (0.167) | 0.374 (0.186) | 0.419 (0.204) |
| Moderate (OR 1.75) effects in all domains +silo specific effects | |||||
| AB Duration | 0.427 (0.014) | 0.713 (0.019) | 0.877 (0.02) | 0.944 (0.02) | 0.961 (0.021) |
Code
d_fig <- d_tbl_1[quant %in% c("ni", "fut_ni") & domain %in% c(2), .SD]
setorderv(d_fig, cols = c("scenario", "domain", "analys", "quant"),
order = c(1, 1, 1, -1))
# d_fig <- rbind(
# d_fig,
# data.table(
# scenario = 1, desc = "Null effect in all domains",
# analys = 1, N = 500, quant = "ni", domain = 1, pr_val = NA
# ),
# data.table(
# scenario = 1, desc = "Null effect in all domains",
# analys = 1, N = 500, quant = "ni", domain = 3, pr_val = NA
# )
# )
# d_fig[, domain := factor(domain,
# levels = c(1, 2, 3),
# labels = c("1", "AB Duration", "3"))]
d_fig[, domain := factor(domain,
levels = c(2),
labels = c("AB Duration"))]
d_fig[, quant := factor(quant,
levels = c("ni", "fut_ni"),
labels = c("NI", "Futility for NI"))]
d_fig[, desc := factor(desc,
levels = unique(d_fig$desc),
labels = unique(d_fig$desc))]
ggplot(d_fig, aes(x = N, y = pr_val, group = quant, col = quant)) +
geom_line(lwd = 0.25) +
scale_y_continuous("", breaks = seq(0, 1, by = 0.2)) +
scale_color_discrete("") +
# facet_grid(desc ~ domain, space = "free") +
ggh4x::facet_grid2(desc ~ domain ,
labeller = labeller(desc = label_wrap_gen(35)),
scales = "free",
axes = "y",
independent = "y") +
ggh4x::force_panelsizes(cols = unit(c(4), "cm"), TRUE) +
# facet_manual(
# . ~ desc, design = matrix(1:17, ncol = 1),
# widths = unit(3, "cm")
# ) +
# theme_minimal() +
theme(
legend.position = "bottom",
strip.text.y.right = element_text(angle = 0,
hjust = 0,
vjust = 0.2,
size = 4),
axis.ticks = element_blank(),
strip.text.x.top = element_text(size = 4),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey",
linewidth = 0.1,
linetype = 1),
axis.text.y = element_text(size = 4),
axis.text.x = element_text(size = 4),
axis.title.x = element_text(size = 5),
legend.text = element_text(size = 4)
)
Sample size of randomised comparisons
Table 3 shows the number of enrolments when any stopping decision is made (or for reaching the maximum 2500 sample size). This is the total number of enrolments that are expected to have occurred, not how many are contributing to the inference in a given domain.
Figure 5 shows the number of participants entering into each of the randomised comparisons by domain and scenario.
The expected cumulative sample sizes are calculated by extracting the number of participants entering into each randomised comparison (i.e. restricted to the relevant strata) and taking the cumulative sum of these over time. For example, the domain 1 (surgical) expected values are based on the participants in the late acute silo that receive randomised surgical intervention. Similarly, the domain 2 (antibiotic duration) expected values are based on the participants across all silos that received one-stage revision and then receive one of the randomised treatment allocations (i.e. are not in the non-randomised treatment group for any reason).
If a decision was made in a domain, then subsequent enrolments would be assigned to the relevant arm and so the allocation would be expected to deviate away from 1:1. For example, if a superiority decision was made for domain 1 (and the trial was still ongoing) then all subsequent participants are assigned to the superior intervention. Finally, if a decision was made for all research questions, the trial stops early. To avoid any bias in the calculation, we use LOCF to propagate the expectations forward through to the maximum number of enrolments.
Code
d_tbl_2_cur <- d_tbl_2[, .(N_mu = mean(N_stopped)), keyby = .(scenario, desc, domain)]
d_tbl_2_cur <- dcast(d_tbl_2_cur, scenario + desc ~ domain, value.var = "N_mu")
d_tbl_2_cur <- d_tbl_2_cur[, .SD, .SDcols = !c("scenario")]
g_tbl <- d_tbl_2_cur |>
gt(groupname_col = "desc") |>
gt::text_transform(
locations = cells_row_groups(),
fn = function(x) {
lapply(x, function(x) {
gt::md(paste0("*", x, "*"))
})
}
) |>
cols_align(
columns = 1,
align = "left"
) |>
cols_align(
columns = 2:ncol(d_tbl_2_cur),
align = "center"
) |>
tab_spanner(
label = html("Expected number of total enrolments to hit stopping rule by domain"),
columns = 2:ncol(d_tbl_2_cur)
) |>
cols_label(
`1` = "Surgical",
`2` = "AB Duration",
`3` = "AB Ext-proph",
`4` = "AB choice"
) |>
tab_options(
table.font.size = "70%"
) |>
fmt_number(decimals = 0, drop_trailing_zeros = TRUE)
g_tbl
Expected number of total enrolments to hit stopping rule by domain
|
|||
|---|---|---|---|
| Surgical | AB Duration | AB Ext-proph | AB choice |
| Null effect in all domains | |||
| 1,026 | 1,695 | 1,108 | 894 |
| Moderate (OR 1.75) surgical revision effect (one and two-stage) | |||
| 943 | 1,553 | 898 | 867 |
| Moderate (OR 1.75) surgical revision effect (one-stage) | |||
| 1,232 | 1,638 | 992 | 853 |
| Moderate (OR 1.75) surgical revision effect (two-stage) | |||
| 1,293 | 1,638 | 1,023 | 898 |
| Moderate (OR 1.75) antibiotic duration 6wk effect | |||
| 1,104 | 1,116 | 1,034 | 920 |
| Moderate (OR 1.75) antibiotic ext-proph 12wk effect | |||
| 1,044 | 1,711 | 1,006 | 893 |
| Moderate (OR 1.75) antibiotic choice rifampacin effect | |||
| 1,031 | 1,730 | 1,090 | 982 |
| Moderate (OR 1.75) effects in all domains | |||
| 848 | 986 | 925 | 1,052 |
| Null effect in all domains +silo specific effects | |||
| 1,250 | 1,704 | 1,008 | 930 |
| Moderate (OR 1.75) surgical revision effect (one and two-stage) +silo specific effects | |||
| 890 | 1,564 | 906 | 921 |
| Moderate (OR 1.75) surgical revision effect (one-stage) +silo specific effects | |||
| 1,266 | 1,602 | 984 | 884 |
| Moderate (OR 1.75) surgical revision effect (two-stage) +silo specific effects | |||
| 1,152 | 740 | 954 | 872 |
| Moderate (OR 1.75) antibiotic duration 6wk effect +silo specific effects | |||
| 1,322 | 1,077 | 994 | 895 |
| Moderate (OR 1.75) antibiotic ext-proph 12wk effect +silo specific effects | |||
| 1,258 | 1,674 | 994 | 853 |
| Moderate (OR 1.75) antibiotic choice rifampacin effect +silo specific effects | |||
| 1,282 | 1,690 | 1,009 | 996 |
| Moderate (OR 1.75) effects in all domains +silo specific effects | |||
| 864 | 984 | 901 | 1,040 |
Code
d_fig_1 <- d_N_by_arm[, .(mu_N = mean(N)), keyby = .(desc, domain, arm, N_enrol)]
d_fig_1[, desc := factor(
desc,
levels = unique(d_N_by_arm$desc))]
d_fig_1[
, domain := factor(
domain,
levels = 1:4,
labels = c("Surgical", "AB Duration", "AB Ext-proph", "AB Choice"))]
d_fig_2 <- d_N_by_arm[
domain == 1 & arm %in% 2:3,
.(N = sum(N), arm = 4),
keyby = .(analys, sim, domain, scenario, desc, N_enrol)]
d_fig_2 <- d_fig_2[, .(mu_N = mean(N)), keyby = .(desc, domain, arm, N_enrol)]
d_fig_2[, desc := factor(
desc,
levels = unique(d_N_by_arm$desc))]
d_fig_2[
, domain := factor(
domain,
levels = 1:4,
labels = c("Surgical", "AB Duration", "AB Ext-proph", "AB Choice"))]
d_fig_1[, arm := factor(arm)]
ggplot(d_fig_1, aes(x = N_enrol, y = mu_N, col = arm)) +
geom_line(lwd = 0.2) +
geom_point(size = 0.4) +
geom_line(
data = d_fig_2,
aes(x = N_enrol, y = mu_N), lwd = 0.2, lty = 2, inherit.aes = F) +
scale_x_continuous("") +
scale_y_continuous("Expected number of participants") +
scale_color_discrete("Treatment arm: ") +
ggh4x::facet_grid2(
desc ~ domain,
labeller = labeller(desc = label_wrap_gen(35)),
scales = "free_y", independent = "y") +
# theme_minimal() +
theme(
legend.position = "bottom",
strip.text.y.right = element_text(angle = 0,
hjust = 0,
vjust = 0.2,
size = 4),
axis.ticks = element_blank(),
strip.text.x.top = element_text(size = 4),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey",
linewidth = 0.1,
linetype = 1),
axis.title.y=element_text(size = 5),
axis.text.y = element_text(size = 4),
axis.text.x = element_text(size = 4),
axis.title.x = element_text(size = 5),
legend.text = element_text(size = 4)
)
Parameter estimation
Figure 6 shows the median value and 95% quantiles from the posterior mean odds ratios obtained from the simulations by each domain and scenario. The estimates are unconditional, in that they ignore the stopping rule and propagate estimates forward with LOCF so that we are working from a random sample (albeit with static imputation) rather than a dependent sample of simulations.
Figure 7 shows the median value and 95% quantiles for the posterior mean risk differences obtained from the simulations for each domain and scenario. The transformation to a risk difference suggests that the odds ratios amount to effects in the order of 10-20% on the risk scale, dependent on the silo, domain etc.
The AB duration domain has high variance in the distribution of posterior means that we anticipate to observe. That is, when there is no true effect, we might still see posterior means as large as \(\pm 0.2\) on the absolute risk scale.
Code
# ggplot2::theme_update(text = element_text(size = 8))
# ggplot2::theme_update(legend.position = "bottom")
# # ggplot2::theme_update(legend.title = element_blank())
# ggplot2::theme_update(axis.text.x = element_text(size = 8))
# ggplot2::theme_update(axis.text.y = element_text(size = 8))
d_fig <- d_tbl_3[,
.(or = median(exp(mu_lor)),
q_025 = quantile(exp(mu_lor), prob = 0.025),
q_975 = quantile(exp(mu_lor), prob = 0.975)),
keyby = .(scenario, desc, analys, domain, N)]
# setorderv(d_fig, cols = "scenario", order = -1L)
d_fig[, desc := factor(desc, levels = unique(d_fig$desc))]
d_fig[, domain := factor(domain,
levels = 1:4,
labels = c("Surgical", "AB Duration", "AB Ext-proph", "AB Choice"))]
# d_tbl_3[scenario == 8, range(lor)]
ggplot(data = d_fig,
aes(x = N, y = or)) +
geom_line(lwd = 0.25) +
geom_errorbar(aes(ymin = q_025, ymax = q_975), lwd = 0.25, width = 50) +
scale_x_continuous("") +
scale_y_continuous("OR") +
ggh4x::facet_grid2(desc ~ domain ,
labeller = labeller(desc = label_wrap_gen(35)),
scales = "free",
axes = "y",
independent = "y") +
# theme_minimal() +
theme(text = element_text(size = 6),
strip.text.y.right = element_text(angle = 0,
hjust = 0,
vjust = 0.2,
size = 4),
strip.text.x = element_text(angle = 0, size = 5),
axis.ticks = element_blank(),
axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size = 5),
axis.text.y = element_text(size = 5),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey",
linewidth = 0.1,
linetype = 1))
Code
d_fig <- d_tbl_3[,
.(rd = median(mu_rd),
q_025 = quantile(mu_rd, prob = 0.025),
q_975 = quantile(mu_rd, prob = 0.975)),
keyby = .(scenario, desc, analys, domain, N)]
# setorderv(d_fig, cols = "scenario", order = -1L)
d_fig[, desc := factor(desc, levels = unique(d_fig$desc))]
d_fig[, domain := factor(domain,
levels = 1:4,
labels = c("Surgical", "AB Duration", "AB Ext-proph", "AB Choice"))]
# d_tbl_3[scenario == 8, range(lor)]
ggplot(data = d_fig,
aes(x = N, y = rd)) +
geom_line(lwd = 0.25) +
geom_errorbar(aes(ymin = q_025, ymax = q_975), lwd = 0.25, width = 50) +
scale_x_continuous("") +
scale_y_continuous("Risk difference") +
ggh4x::facet_grid2(desc ~ domain ,
labeller = labeller(desc = label_wrap_gen(35)),
scales = "free",
axes = "y",
independent = "y") +
# theme_minimal() +
theme(text = element_text(size = 6),
strip.text.y.right = element_text(angle = 0,
hjust = 0,
vjust = 0.2,
size = 4),
strip.text.x = element_text(angle = 0, size = 5),
axis.ticks = element_blank(),
axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size = 5),
axis.text.y = element_text(size = 5),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey",
linewidth = 0.1,
linetype = 1))
Fraction of uncertainty resolved
Figure 8 shows the median value and 95% quantiles for the fraction of uncertainty resolved based on
\[ \begin{aligned} \text{Fraction \ Resolved} = 1 - \frac{Var(\beta_{post})}{Var(\beta_{pri})} \end{aligned} \]
where \(Var(\beta_{post})\) and \(Var(\beta_{pri})\) represent the variance associated with the prior and posterior belief for the relevant log-odds ratio for the treatment effects. This is basically just a way to compare the prior and posterior variance. When the posterior is based on negligible data, the variance will be similar to that of the prior and the fraction resolved will be very small. A low fraction resolved (e.g. less than 0.5) suggests that any decision that was made was done so with a substantial amount of uncertainty remaining (you didn’t move far from your prior belief) whereas values close to unity suggest that a lot of the uncertainty has been resolved.
- What is obvious from the above plots is also obvious here, the decision made in the AB duration domain are subject to a substantial amount of uncertainty.
Code
l_cfg <- copy(l[[1]]$cfg)
# initial uncertainty
v0 <- (l_cfg$pri$b_trt[2])^2
d_fig <- copy(d_tbl_3)
d_fig[, fr_unc := 1 - (se_lor/v0)]
d_fig <- d_fig[,
.(fr_unc = median(fr_unc),
q_025 = quantile(fr_unc, prob = 0.025),
q_975 = quantile(fr_unc, prob = 0.975)),
keyby = .(scenario, desc, analys, domain, N)]
# setorderv(d_fig, cols = "scenario", order = -1L)
d_fig[, desc := factor(desc, levels = unique(d_fig$desc))]
d_fig[, domain := factor(domain,
levels = 1:4,
labels = c("Surgical", "AB Duration", "AB Ext-proph", "AB Choice"))]
# d_tbl_3[scenario == 8, range(lor)]
ggplot(data = d_fig,
aes(x = N, y = fr_unc)) +
geom_line(lwd = 0.25) +
geom_errorbar(aes(ymin = q_025, ymax = q_975), lwd = 0.25, width = 50) +
scale_x_continuous("") +
scale_y_continuous("Fraction of uncertainty resolved (1-(V_post/V_pri))") +
ggh4x::facet_grid2(desc ~ domain ,
labeller = labeller(desc = label_wrap_gen(35)),
scales = "free",
axes = "y",
independent = "y") +
# theme_clean() +
# theme_minimal() +
theme(text = element_text(size = 6),
strip.text.y.right = element_text(angle = 0,
hjust = 0,
vjust = 0.2,
size = 4),
strip.text.x = element_text(angle = 0, size = 5),
axis.ticks = element_blank(),
axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size = 5),
axis.text.y = element_text(size = 5),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey",
linewidth = 0.1,
linetype = 1))