Marginal-ish effects
The following DAG includes surgical, antibiotic backbone, extended prophylactic and choice domains. Additionally, site of infection (knee/hip) and baseline clinician preferences are included where the latter are assumed to take primacy in selection of surgery type etc. The setup is intended to represent the approximate structure, not to replicate the exact data generating process.
Preference for surgery primarily impacts those assigned to revision but would also be applicable for a patient that ends up switching treatment, for example switch from dair to revision due to a irreparable loose joint. Preference for antibiotic duration comes into play for dair and two-stage patients since these units are never randomly assigned to duration. For example, if a unit received dair, then the preference would lead to some duration of backbone antibiotic therapy. Analogously, preference for extended prophylactic applies where the unit didn’t enter the extended prophylaxis domain even though they received a two-stage revision; here prophylaxis would be determined by the clinician assessment of what is best. For antibiotic choice (rifampicin), it is conceivable that patients are eligible for the domain but the clinician is not willing to risk them being allocated to the control arm, hence the possibility of clinical selection is included for that domain as well. Only a subset of the cohort will ever enter the choice domain anyway.
A joint-by-surgical type (dair, one, two-stage) interaction is explicitly included (denoted with \(S_A \times J\)). There may be others that are of interest. My guess is that infection site will also influence the approach adopted for revision, but that is not reflected here (I don’t think it would make any difference to the statistical considerations anyway).
Previously, we have talked about direct paths, but these are assumed not to exist. This is because you cannot define, for example surgical procedure, as a mediator because \(Y(a,m)\) would need to be defined for both \(a=0\) and \(a=1\) and it is only makes sense in the context of the latter. There are therefore no direct paths within the causal structure.
Mentally I am restricting myself to the late silo and specifically to the question of the surgical intervention and hence conditioning on that silo is implicit here and not mentioned within the DAG or models.
Given the above, some of the following may not align with a generic perspective applicable to the other silos.
Throughout, I will represent the nodes with:
- site of infection (knee, hip): \(J\)
- assigned surgery: \(A\)
- clinical selection for surgery: \(U_A\)
- received surgery: \(S_A\)
- assigned duration: \(D1\)
- clinical selection for duration: \(U_{D1}\)
- received duration: \(S_{D1}\)
- assigned prohphylaxis: \(D2\)
- clinical selection for prohphylaxis: \(U_{D2}\)
- received prohphylaxis: \(S_{D2}\)
- assigned choice: \(C\)
- clinical selection for choice: \(U_{C}\)
- received prohphylaxis: \(S_{C}\)
- outcome: \(Y\)
The total effect of surgery (\(A\)) on outcome (\(Y\)) is identifiable. However, the total effect is the effect of treatment through all paths, so it amounts to a comparison of surgery plus some possible antibiotic duration and extended prophylaxis, i.e. all the open paths below.
paths open
[1,] "a -> sa -> axj -> y" "TRUE"
[2,] "a -> sa -> axj <- j -> y" "FALSE"
[3,] "a -> sa -> sd1 -> y" "TRUE"
[4,] "a -> sa -> sd1 <- ud1 -> y" "FALSE"
[5,] "a -> sa -> sd2 -> y" "TRUE"
[6,] "a -> sa -> sd2 <- ud2 -> y" "FALSE"
[7,] "a -> sa -> y" "TRUE"
[8,] "a -> sa <- ua -> y" "FALSE"No adjustment is necessary to obtain the total effect but under the design this obfuscates attribution due to the various interventions.
Surgical domain model
Per the above, we could simply regress on assigned surgical treatment to obtain the total effect. However, for the surgical domain it is informative to use the following model specification:
\[ \begin{aligned} \mathbb{E}[Y | \pmb{Q}; \pmb{\zeta}] &= \beta_0 + \beta_1 (S_A = 1) + \beta_2 (S_A = 2) + \beta_3 U_A + \\ & \beta_4 S_{D1} + \beta_5 U_{D1} + \\ & \beta_6 S_{C} + \beta_7 U_{C} + \\ & \beta_8 J + \\ & \beta_9 S_{D2} + \beta_{10} U_{D2} + \\ & \gamma_1 (S_A = 1) J + \gamma_2 (S_A = 2) J \end{aligned} \]
where \(\pmb{Q}\) and \(\pmb{\zeta}\) represent the vector of conditional elements and with the parameters for \(D2\) (\(\beta_9\) and \(\beta_{10}\)) being relevant only for the cohort receiving two-stage revision.
To get back to a view on dair vs revision we can use the g-formula to average over all the other terms and then compute a weighted combination of the one-stage and two-stage effects based on the distribution of these in the sample. For example, say the only conditioning element we have is \(U_A\) and \(Y(a)\) denotes the potential outcomes with respect to the randomised surgery type:
\[ \begin{aligned} \mathbb{E}[Y(0)] &= \mathbb{E}_{U_A}[ \mathbb{E}[Y | S_A = 0, U_A = u_a] ] \\ &= \sum_{u_A} \mathbb{P}(U_A = u_a) \left(\sum_{y} y \mathbb{P}(Y = y | S_A = 0, U_A = u_a) \right) \end{aligned} \]
\[ \begin{aligned} \mathbb{E}[Y(1)] &= \mathbb{E}_{U_A}[ \mathbb{E}[Y | S_A \in \{ 1, 2\}, U_A = u_a] ] \\ &= \sum_{u_A} \mathbb{P}(U_A = u_a) \left(\sum_{y} y \mathbb{P}(Y = y | S_A = \in \{ 1, 2\}, U_A = u_a) \right) \end{aligned} \]
The model is implemented in stan as follows:
Logistic regression with bootstrap in stan
data {
int N;
array[N] int y;
array[N] int n;
// for one an two-stage indicators
// these are the main effects of interest
int P1;
matrix[N, P1] X1;
// other terms bar interactions and ext proph related
int P2;
matrix[N, P2] X2;
// extra fields for extended proph
int P3;
matrix[N, P3] X3;
// fields for interactions
int P4;
matrix[N, P4] X4;
// auxiliary terms
// observed probability of one-stage
real pr_one;
// number preferring 1/2 stage under rev
int N1;
int N2;
array[N1] int ix1;
array[N2] int ix2;
// priors
int prior_only;
real sd_a0;
real sd_b1;
real sd_b2;
real sd_b3;
real sd_b4;
}
transformed data {
}
parameters{
real a0;
// surgical effects
vector[P1] b1;
// other terms
vector[P2] b2;
// ext prophy
vector[P3] b3;
// interactions
vector[P4] b4;
}
transformed parameters{
vector[N] eta;
for(i in 1:N){
// X1[i,2] should be an indicator of whether two-stage has been done
if(X1[i,2] == 0){
eta[i] = a0 + X1[i,]*b1 + X2[i,]*b2 + X4[i, ]*b4 ;
}
if(X1[i,2] == 1){
eta[i] = a0 + X1[i,]*b1 + X2[i,]*b2 + X3[i,]*b3 + X4[i, ]*b4 ;
}
}
}
model{
target += normal_lpdf(a0 | 0, sd_a0);
target += normal_lpdf(b1 | 0, sd_b1);
target += normal_lpdf(b2 | 0, sd_b2);
target += normal_lpdf(b3 | 0, sd_b3);
target += normal_lpdf(b4 | 0, sd_b4);
if(!prior_only){
target += binomial_logit_lpmf(y | n, eta);
}
}
generated quantities{
vector[N] p;
real marg_p0;
real marg_p1;
real rd;
vector[N] cond_p0;
vector[N] cond_p1;
// to get to pate rather than sate
// Bayesian bootstrap (weights)
vector[N] w = dirichlet_rng(to_vector(n));
p = inv_logit(eta);
// drop out the terms relating to sa == 1 and sa == 2
// drop out last two terms that are only applicable for sa == 2
// dair - no X4 since sa = 0 hence zero contribution.
cond_p0 = inv_logit(a0 + X2*b2);
// one
cond_p1[ix1] = inv_logit(a0 + b1[1] + X2[ix1,]*b2 + X4[ix1,]*b4);
// two
cond_p1[ix2] = inv_logit(a0 + b1[2] + X2[ix2,]*b2 + X3[ix2,]*b3 + X4[ix2,]*b4);
// taking average over bayesian bootstrap weights
marg_p0 = w' * cond_p0;
marg_p1 = w' * cond_p1;
rd = marg_p1 - marg_p0;
}where the bootstrap functionality is used as a standardisation procedure and which could be extended to estimate strata level effects such as those relating to site of infection (knee/hip).
Below are scenarios with different combinations of effects that aims to illustrate the main limitation of the design. The assumed effect sizes are purely illustrative.
Example 1 - null case
Generic data generation function
n = 200
get_data <- function(
n = 200,
ff = function(j, a, ua, sa, d1, ud1, sd1, d2, ud2, sd2, c, uc, sc){
# provide a sensible linear predictor.
p <- rep(0.5, length(a))
p
})
{
# joint
j <- rbinom(n, 1, 0.5)
# choice
c <- rbinom(n, 1, 0.5)
# The clinical perspective overlay encapsulate entry into the domain.
# Entry implies the unit meets the conditions for rand to rif
uc <- rbinom(n, 1, 0.6)
# Type received
sc <- rep(NA, n)
sc[uc == 0] <- c[uc == 0]
# those that do not enter do not get rif
sc[uc == 1] <- 0
# assigned surgery - just assume that all enter into this domain
a <- rbinom(n, 1, 0.5)
# under receipt of rev, most clinicians would prefer to use two-stage
ua <- rbinom(n, 1, 0.75)
# dair = 0, one = 1, two = 2
sa <- rep(NA, n)
# assume full adherence in dair group
sa[a == 0] <- 0
# assigned to rev: gets either one or two stage
ix <- which(a == 1);
sa[ix] <- ua[ix] + 1
# random assignment of duration 6vs12
d1 <- sample(c(6,12), n, T)
# pref for duration of therapy after first op is towards longer durn
ud1 <- sample(seq(4, 12, by = 2), n, T, 1:5/sum(1:5))
# random assignment of prophylaxis
d2 <- rbinom(n, 1, 0.5)
# most prefer to use prophylaxis
ud2 <- rbinom(n, 1, 0.7)
# duration of antibiotic following first procedure
sd1 <- rep(NA, n)
# for those getting dair, duration after first opn is exactly per pref
sd1[sa == 0] <- ud1[sa == 0]
# assume full adherence in one-stage group
# 6 wk vs 12 wk
sd1[sa == 1] <- d1[sa == 1]
# for those getting two-stage, duration after first opn is per pref
sd1[sa == 2] <- ud1[sa == 2]
# prophylaxis
sd2 <- rep(NA, n)
# for those getting two-stage it is per the randomisation
sd2[sa == 2] <- d2[sa == 2]
# otherwise it isn't defined but assign it to what it would be set to
# if the unit were to get two-stage
# these units will not be used in the likelihood, but they would be used
# in the bootstrap step.
sd2[sa != 2] <- d2[sa != 2]
d <- data.table(j, a, ua, sa, d1, ud1, sd1, d2, ud2, sd2, c, uc, sc)
# outcome
d[, p_y := ff(j, a, ua, sa, d1, ud1, sd1, d2, ud2, sd2, c, uc, sc)]
d[, y := rbinom(.N, 1, prob = p_y)]
d
}Here all domain treatment effects are set to zero, however, the some baseline variables have been assumed to create heterogeneity in the response. Rather than repeat simulation, the model is fit to a large simulated dataset to get a sense of consistency.
Data simulation and parameter estimation
set.seed(987654321)
g_pars <- c(paste0("b1[",1:2,"]"),
paste0("b2[",1:6,"]"),
paste0("b3[",1:2,"]"),
paste0("b4[",1:2,"]")
)
names(g_pars) <- c(
"b1 (sa1)", "b2 (sa2)", "b3 (ua)",
"b4 (sd1)", "b5 (ud1)",
"b6 (sc)", "b7 (uc)",
"b8 (j)",
"b9 (sd2)", "b10 (ud2)",
"g1 (sa1 x j)", "g2 (sa2 x j)"
)
m1 <- cmdstanr::cmdstan_model("stan/logistic-demo-03.stan")
m2 <- cmdstanr::cmdstan_model("stan/logistic-demo-02.stan")
d <- get_data(
n = 1e5,
ff = function(j, a, ua, sa, d1, ud1, sd1, d2, ud2, sd2, c, uc, sc){
p <- rep(NA, length(a))
# dair
ix <- which(sa == 0)
# e.g.
# ua: pref for two-stage suggests less sev pt
# sa: one-stage does nothing
# sa: two-stage does nothing
# ud1: pref for longer duration suggests more sev pt
# d1: duration does nothing
# d2: is irrelevant for dair cohort
# uc: does nothing
# c: does nothing
p[ix] <- 0.6 +
0*(sa[ix] == 1) + 0*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
# one-stage
ix <- which(sa == 1)
p[ix] <- 0.6 +
0*(sa[ix] == 1) + 0*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
# two-stage
ix <- which(sa == 2)
p[ix] <- 0.6 +
0*(sa[ix] == 1) + 0*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
# last two parameters not included for dair and one-stage
0.0*sd2[ix] + 0*ud2[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
p
})
# d[, .(.N, mean(y)), keyby = a]
d[, `:=`(sa1 = as.numeric(d$sa == 1), sa2 = as.numeric(d$sa == 2))]
d_s <- d[, .(y = sum(y), n = .N, p_y = unique(p_y)), keyby = .(
sa1, sa2, ua, sd1, ud1, sc, uc, j, sd2, ud2, j
)]
# Model treats X[, 2] as being indicator of two-stage
X1 <- cbind(d_s$sa1, d_s$sa2)
# all others bar interactions
X2 <- cbind(d_s$ua, d_s$sd1, d_s$ud1, d_s$sc, d_s$uc, d_s$j)
# extended prophylaxis terms - only applies to two-stage pt
X3 <- cbind(d_s$sd2, d_s$ud2)
# interactions
X4 <- cbind(d_s$j*d_s$sa1, d_s$j*d_s$sa2)
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X1 = X1, X2 = X2, X3 = X3, X4 = X4,
P1 = ncol(X1), P2 = ncol(X2), P3 = ncol(X3), P4 = ncol(X4),
# number of recs preference for 1/2 under revision
N1 = d_s[ua == 0, .N], N2 = d_s[ua == 1, .N],
# idx for pref for 1/2
ix1 = d_s[ua == 0, which = T],
ix2 = d_s[ua == 1, which = T],
pr_one = d[, mean(ua == 0)],
prior_only = 0,
sd_a0 = 1.5,
sd_b1 = 3, sd_b2 = 3, sd_b3 = 3, sd_b4 = 3
)
f1 <- m1$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 1 finished in 11.5 seconds.
Chain 2 finished in 12.2 seconds.
Both chains finished successfully.
Mean chain execution time: 11.9 seconds.
Total execution time: 12.4 seconds.Data simulation and parameter estimation
# snk <- capture.output(
# f1 <- m1$pathfinder(ld, num_paths=20, single_path_draws=200,
# history_size=50, max_lbfgs_iters=100,
# refresh = 0, draws = 2000)
# )
# f_mode <- m1$optimize(data = ld, jacobian = TRUE, refresh =0)
# f1 <- m1$laplace(data = ld, mode = f_mode, refresh = 0)
d_s <- d[, .(y = sum(y), n = .N), keyby = .(
a
)]
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X = matrix(d_s$a),
prior_only = 0,
sd_a0 = 1.5,
sd_b = 3
)
f2 <- m2$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Both chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.1 seconds.The parameter estimates obviously do not align with the parameters used in the linear predictor as these are estimated on the log odds scale. But the magnitudes of the effects should be in line with the parameters on the risk scale. The parameters that had non-null effects in the data generation process have been shaded red.
Posterior distributions for model parameters
d_post <- data.table(f1$draws(variables = c("a0","b1","b2","b3","b4"), format = "matrix"))
names(d_post) <- c("a0", names(g_pars))
# parameters
cols <- c("a0", names(g_pars))
d_fig <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig[variable == "a0", value := plogis(value)]
d_fig[variable != "a0", value := exp(value)]
d_col <- unique(d_fig[, .(variable)])
d_col[, fill := factor(c(
1, 1, 1, 2,
1, 2,
1, 1,
2,
1, 1,
1, 1
))]
ggplot(d_fig, aes(x = value, group = variable)) +
geom_rect(data = d_col, aes(fill = fill),
xmin = -Inf,xmax = Inf,
ymin = -Inf,ymax = Inf, alpha = 0.1,
inherit.aes = F) +
scale_fill_manual("", values = c("white", "red")) +
geom_density() +
geom_vline(data = d_fig[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
scale_x_continuous(
labels = label_number(accuracy = 0.01)
) +
facet_wrap(~variable, scales = "free") +
theme(legend.position = "none")
Below are the recovered marginal probabilities of \(Y\) under each surgical treatment type (dair vs rev). The red overlay shows the marginal probabilities of \(Y\) based on using a truncated model that includes a single fixed effect for surgical intervention (dair vs rev) and no other terms.
Marginal probability of outcome by surgical type
cols <- c("marg_p0", "marg_p1")
d_post <- data.table(f1$draws(variables = cols, format = "matrix"))
names(d_post) <- cols
d_fig1 <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig2 <- data.table(f2$draws(variables = cols, format = "matrix"))
d_fig2 <- melt(d_fig2, measure.vars = cols)
ggplot(d_fig1, aes(x = value, group = variable)) +
geom_density() +
geom_density(data = d_fig2, aes(x = value, group = variable),
col = 2) +
geom_vline(data = d_fig1[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
facet_wrap(~variable, ncol = 1) +
scale_x_continuous(
"Proportion with successful outcome",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
)
Risk difference (rev - dair). The red dashed line shows the posterior distribution obtained from a reduced model that includes a single term for surgical treatment assignment (dair/rev). The result is centred around zero and aligns with what we expect.
Risk difference (dair vs revision)
d_fig1 <- data.table(f1$draws(variables = c("rd"), format = "matrix"))
d_fig2 <- data.table(f2$draws(variables = c("rd"), format = "matrix"))
ggplot(d_fig1, aes(x = rd)) +
geom_density() +
geom_density(data = d_fig2,
aes(x = rd), col = 2, lty = 2) +
scale_x_continuous(
"Risk difference",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
) +
geom_vline(data = d_fig1[, .(mu = mean(rd))],
aes(xintercept = mu)) +
geom_vline(data = d_fig2[, .(mu = mean(rd))],
aes(xintercept = mu), col = 2) 
Example 2 - revision effects
Here equal magnitude effects are used for revision for both one-stage and two-stage approach with worse outcomes for hip joints but no interactions. In the other domains, the treatment effects are set to zero.
Data simulation and parameter estimation
m1 <- cmdstanr::cmdstan_model("stan/logistic-demo-03.stan")
m2 <- cmdstanr::cmdstan_model("stan/logistic-demo-02.stan")
d <- get_data(
n = 1e5,
ff = function(j, a, ua, sa, d1, ud1, sd1, d2, ud2, sd2, c, uc, sc){
p <- rep(NA, length(a))
# dair
ix <- which(sa == 0)
p[ix] <- 0.6 +
0.1*(sa[ix] == 1) + 0.1*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
-0.00*j[ix]**(sa[ix] == 1) -0.00*j[ix]**(sa[ix] == 2)
# one-stage
ix <- which(sa == 1)
p[ix] <- 0.6 +
0.1*(sa[ix] == 1) + 0.1*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
-0.00*j[ix]**(sa[ix] == 1) -0.00*j[ix]**(sa[ix] == 2)
# two-stage
ix <- which(sa == 2)
p[ix] <- 0.6 +
0.1*(sa[ix] == 1) + 0.1*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
# last two parameters not included for dair and one-stage
0.0*sd2[ix] + 0*ud2[ix] +
-0.00*j[ix]**(sa[ix] == 1) -0.00*j[ix]**(sa[ix] == 2)
# p[ix] <- 0.6 +
# 0.1*(sa[ix] == 1) + 0.1*(sa[ix] == 2) + 0.05*ua[ix] +
# 0*sd1[ix] - 0.01*ud1[ix] +
# 0*sc[ix] + 0*uc[ix] +
# -0.1*j[ix] +
# # last two parameters not included for dair and one-stage
# 0.0*sd2[ix] + 0*ud2[ix] +
# -0.03*j[ix]**(sa[ix] == 1) -0.03*j[ix]**(sa[ix] == 2)
p
})
d[, `:=`(sa1 = as.numeric(d$sa == 1), sa2 = as.numeric(d$sa == 2))]
d_s <- d[, .(y = sum(y), n = .N, p_y = unique(p_y)), keyby = .(
sa1, sa2, ua, sd1, ud1, sc, uc, j, sd2, ud2, j
)]
# Model treats X[, 2] as being indicator of two-stage
X1 <- cbind(d_s$sa1, d_s$sa2)
# all others bar interactions
X2 <- cbind(d_s$ua, d_s$sd1, d_s$ud1, d_s$sc, d_s$uc, d_s$j)
# extended prophylaxis terms - only applies to two-stage pt
X3 <- cbind(d_s$sd2, d_s$ud2)
# interactions
X4 <- cbind(d_s$j*d_s$sa1, d_s$j*d_s$sa2)
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X1 = X1, X2 = X2, X3 = X3, X4 = X4,
P1 = ncol(X1), P2 = ncol(X2), P3 = ncol(X3), P4 = ncol(X4),
# number of recs preference for 1/2 under revision
N1 = d_s[ua == 0, .N], N2 = d_s[ua == 1, .N],
# idx for pref for 1/2
ix1 = d_s[ua == 0, which = T],
ix2 = d_s[ua == 1, which = T],
pr_one = d[, mean(ua == 0)],
prior_only = 0,
sd_a0 = 1.5,
sd_b1 = 3, sd_b2 = 3, sd_b3 = 3, sd_b4 = 3
)
f1 <- m1$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 1 finished in 11.0 seconds.
Chain 2 finished in 11.2 seconds.
Both chains finished successfully.
Mean chain execution time: 11.1 seconds.
Total execution time: 11.3 seconds.Data simulation and parameter estimation
# snk <- capture.output(
# f1 <- m1$pathfinder(ld, num_paths=20, single_path_draws=200,
# history_size=50, max_lbfgs_iters=100,
# refresh = 0, draws = 2000)
# )
# f_mode <- m1$optimize(data = ld, jacobian = TRUE, refresh =0)
# f1 <- m1$laplace(data = ld, mode = f_mode, refresh = 0)
d_s <- d[, .(y = sum(y), n = .N), keyby = .(
a
)]
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X = matrix(d_s$a),
prior_only = 0,
sd_a0 = 1.5,
sd_b = 10
)
f2 <- m2$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Both chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.1 seconds.Parameters estimated on the log-odds scale, but exponentiated below (intercept reported on risk scale).
Posterior distributions for model parameters
d_post <- data.table(f1$draws(variables = c("a0","b1","b2","b3","b4"), format = "matrix"))
names(d_post) <- c("a0", names(g_pars))
# parameters
cols <- c("a0", names(g_pars))
d_fig <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig[variable == "a0", value := plogis(value)]
d_fig[variable != "a0", value := exp(value)]
d_col <- unique(d_fig[, .(variable)])
d_col[, fill := factor(c(
1, 2, 2, 2,
1, 2,
1, 1,
2,
1, 1,
2, 2
))]
ggplot(d_fig, aes(x = value, group = variable)) +
geom_rect(data = d_col, aes(fill = fill),
xmin = -Inf,xmax = Inf,
ymin = -Inf,ymax = Inf, alpha = 0.1,
inherit.aes = F) +
scale_fill_manual("", values = c("white", "red")) +
geom_density() +
geom_vline(data = d_fig[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
scale_x_continuous(
labels = label_number(accuracy = 0.01)
) +
facet_wrap(~variable, scales = "free") +
theme(legend.position = "none")
The distribution of the marginal probability of \(Y\) under each treatment type.
Marginal probability of outcome by surgical type
cols <- c("marg_p0", "marg_p1")
d_post <- data.table(f1$draws(variables = cols, format = "matrix"))
names(d_post) <- cols
d_fig1 <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig2 <- data.table(f2$draws(variables = cols, format = "matrix"))
d_fig2 <- melt(d_fig2, measure.vars = cols)
ggplot(d_fig1, aes(x = value, group = variable)) +
geom_density() +
geom_density(data = d_fig2, aes(x = value, group = variable),
col = 2) +
geom_vline(data = d_fig1[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
facet_wrap(~variable, ncol = 1) +
scale_x_continuous(
"Proportion with successful outcome",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
)
Risk difference (rev - dair). The red dashed line shows the posterior distribution obtained from a reduced model that includes a single term for surgical treatment assignment (dair/rev). The result is suggests a positive effect of surgery as expected.
Risk difference (dair vs revision)
d_fig1 <- data.table(f1$draws(variables = c("rd"), format = "matrix"))
d_fig2 <- data.table(f2$draws(variables = c("rd"), format = "matrix"))
ggplot(d_fig1, aes(x = rd)) +
geom_density() +
scale_x_continuous(
"Risk difference",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
) +
geom_density(data = d_fig2,
aes(x = rd), col = 2, lty = 2) +
geom_vline(data = d_fig1[, .(mu = mean(rd))],
aes(xintercept = mu)) +
geom_vline(data = d_fig2[, .(mu = mean(rd))],
aes(xintercept = mu), col = 2) 
Example 3 - extended prophylaxis effects
Here revision effects are set to zero but a non-zero effect is set for the extended prophylaxis.
This is where things start to go awry.
Due to the fact that extended prophylaxis is only applicable for patients receiving a two-stage revision (the revision being the important part), any extended prophylaxis effect will manifest as being an effect of revision. Thus the observed effect has nothing to do with the surgery performed but rather the intervention received in the entangled domain.
Is there any way around this in the setting where we are referring to a generalised revision effect?
Data simulation and parameter estimation
m1 <- cmdstanr::cmdstan_model("stan/logistic-demo-03.stan")
m2 <- cmdstanr::cmdstan_model("stan/logistic-demo-02.stan")
d <- get_data(
n = 1e5,
ff = function(j, a, ua, sa, d1, ud1, sd1, d2, ud2, sd2, c, uc, sc){
p <- rep(NA, length(a))
# dair
ix <- which(sa == 0)
p[ix] <- 0.6 +
0*(sa[ix] == 1) + 0*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
# one-stage
ix <- which(sa == 1)
p[ix] <- 0.6 +
0*(sa[ix] == 1) + 0*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
# two-stage
ix <- which(sa == 2)
p[ix] <- 0.6 +
0*(sa[ix] == 1) + 0*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
# last two parameters not included for dair and one-stage
0.2*sd2[ix] + 0*ud2[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
p
})
d[, `:=`(sa1 = as.numeric(d$sa == 1), sa2 = as.numeric(d$sa == 2))]
d_s <- d[, .(y = sum(y), n = .N, p_y = unique(p_y)), keyby = .(
sa1, sa2, ua, sd1, ud1, sc, uc, j, sd2, ud2, j
)]
# Model treats X[, 2] as being indicator of two-stage
X1 <- cbind(d_s$sa1, d_s$sa2)
# all others bar interactions
X2 <- cbind(d_s$ua, d_s$sd1, d_s$ud1, d_s$sc, d_s$uc, d_s$j)
# extended prophylaxis terms - only applies to two-stage pt
X3 <- cbind(d_s$sd2, d_s$ud2)
# interactions
X4 <- cbind(d_s$j*d_s$sa1, d_s$j*d_s$sa2)
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X1 = X1, X2 = X2, X3 = X3, X4 = X4,
P1 = ncol(X1), P2 = ncol(X2), P3 = ncol(X3), P4 = ncol(X4),
# number of recs preference for 1/2 under revision
N1 = d_s[ua == 0, .N], N2 = d_s[ua == 1, .N],
# idx for pref for 1/2
ix1 = d_s[ua == 0, which = T],
ix2 = d_s[ua == 1, which = T],
pr_one = d[, mean(ua == 0)],
prior_only = 0,
sd_a0 = 1.5,
sd_b1 = 3, sd_b2 = 3, sd_b3 = 3, sd_b4 = 3
)
f1 <- m1$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 2 finished in 11.3 seconds.
Chain 1 finished in 12.3 seconds.
Both chains finished successfully.
Mean chain execution time: 11.8 seconds.
Total execution time: 12.4 seconds.Data simulation and parameter estimation
# snk <- capture.output(
# f1 <- m1$pathfinder(ld, num_paths=20, single_path_draws=200,
# history_size=50, max_lbfgs_iters=100,
# refresh = 0, draws = 2000)
# )
# f_mode <- m1$optimize(data = ld, jacobian = TRUE, refresh =0)
# f1 <- m1$laplace(data = ld, mode = f_mode, refresh = 0)
d_s <- d[, .(y = sum(y), n = .N), keyby = .(
a
)]
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X = matrix(d_s$a),
prior_only = 0,
sd_a0 = 1.5,
sd_b = 3
)
f2 <- m2$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Both chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.1 seconds.Parameters estimated on the log-odds scale, but exponentiated below (intercept reported on risk scale).
Posterior distributions for model parameters
d_post <- data.table(f1$draws(variables = c("a0","b1","b2","b3","b4"), format = "matrix"))
names(d_post) <- c("a0", names(g_pars))
# parameters
cols <- c("a0", names(g_pars))
d_fig <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig[variable == "a0", value := plogis(value)]
d_fig[variable != "a0", value := exp(value)]
d_col <- unique(d_fig[, .(variable)])
d_col[, fill := factor(c(
1, 1, 1, 2,
1, 2,
1, 1,
2,
2, 1,
1, 1
))]
ggplot(d_fig, aes(x = value, group = variable)) +
geom_rect(data = d_col, aes(fill = fill),
xmin = -Inf,xmax = Inf,
ymin = -Inf,ymax = Inf, alpha = 0.1,
inherit.aes = F) +
scale_fill_manual("", values = c("white", "red")) +
geom_density() +
geom_vline(data = d_fig[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
scale_x_continuous(
labels = label_number(accuracy = 0.01)
) +
facet_wrap(~variable, scales = "free") +
theme(legend.position = "none")
The distribution of the marginal probability of \(Y\) under each treatment type.
Marginal probability of outcome by surgical type
cols <- c("marg_p0", "marg_p1")
d_post <- data.table(f1$draws(variables = cols, format = "matrix"))
names(d_post) <- cols
d_fig1 <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig2 <- data.table(f2$draws(variables = cols, format = "matrix"))
d_fig2 <- melt(d_fig2, measure.vars = cols)
ggplot(d_fig1, aes(x = value, group = variable)) +
geom_density() +
geom_density(data = d_fig2, aes(x = value, group = variable),
col = 2) +
geom_vline(data = d_fig1[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
facet_wrap(~variable, ncol = 1) +
scale_x_continuous(
"Proportion with successful outcome",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
)
Risk difference (rev - dair). The red dashed line shows the posterior distribution obtained from a reduced model that includes a single term for surgical treatment assignment (dair/rev). Even though there is no effect of surgery, a positive risk difference for dair vs revision manifests due to the presence of an effect in the extended prophylaxis domain.
Risk difference (dair vs revision)
d_fig1 <- data.table(f1$draws(variables = c("rd"), format = "matrix"))
d_fig2 <- data.table(f2$draws(variables = c("rd"), format = "matrix"))
ggplot(d_fig1, aes(x = rd)) +
geom_density() +
scale_x_continuous(
"Risk difference",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
) +
geom_density(data = d_fig2,
aes(x = rd), col = 2, lty = 2) +
geom_vline(data = d_fig1[, .(mu = mean(rd))],
aes(xintercept = mu)) +
geom_vline(data = d_fig2[, .(mu = mean(rd))],
aes(xintercept = mu), col = 2) 
Example 4 - combined revision and extended prophylaxis effects
Here revision effects are set to negative values (i.e. people on revision do worse) and a positive effect is set for the extended prophylaxis.
For similar reasons as above, the negative effect of revision (people do worse under surgical removal and replacement of the joint) is cancelled out when evaluated at the margins due to the offsetting effect of extended prophylaxis.
Data simulation and parameter estimation
m1 <- cmdstanr::cmdstan_model("stan/logistic-demo-03.stan")
m2 <- cmdstanr::cmdstan_model("stan/logistic-demo-02.stan")
d <- get_data(
n = 1e5,
ff = function(j, a, ua, sa, d1, ud1, sd1, d2, ud2, sd2, c, uc, sc){
p <- rep(NA, length(a))
# dair
ix <- which(sa == 0)
p[ix] <- 0.6 +
-0.1*(sa[ix] == 1) + -0.1*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
# one-stage
ix <- which(sa == 1)
p[ix] <- 0.6 +
-0.07*(sa[ix] == 1) + -0.07*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
# two-stage
ix <- which(sa == 2)
p[ix] <- 0.6 +
-0.07*(sa[ix] == 1) + -0.07*(sa[ix] == 2) + 0.05*ua[ix] +
0*sd1[ix] - 0.01*ud1[ix] +
0*sc[ix] + 0*uc[ix] +
-0.1*j[ix] +
# last two parameters not included for dair and one-stage
0.2*sd2[ix] + 0*ud2[ix] +
0*j[ix]**(sa[ix] == 1) + 0*j[ix]**(sa[ix] == 2)
p
})
d[, `:=`(sa1 = as.numeric(d$sa == 1), sa2 = as.numeric(d$sa == 2))]
d_s <- d[, .(y = sum(y), n = .N, p_y = unique(p_y)), keyby = .(
sa1, sa2, ua, sd1, ud1, sc, uc, j, sd2, ud2, j
)]
# Model treats X[, 2] as being indicator of two-stage
X1 <- cbind(d_s$sa1, d_s$sa2)
# all others bar interactions
X2 <- cbind(d_s$ua, d_s$sd1, d_s$ud1, d_s$sc, d_s$uc, d_s$j)
# extended prophylaxis terms - only applies to two-stage pt
X3 <- cbind(d_s$sd2, d_s$ud2)
# interactions
X4 <- cbind(d_s$j*d_s$sa1, d_s$j*d_s$sa2)
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X1 = X1, X2 = X2, X3 = X3, X4 = X4,
P1 = ncol(X1), P2 = ncol(X2), P3 = ncol(X3), P4 = ncol(X4),
# number of recs preference for 1/2 under revision
N1 = d_s[ua == 0, .N], N2 = d_s[ua == 1, .N],
# idx for pref for 1/2
ix1 = d_s[ua == 0, which = T],
ix2 = d_s[ua == 1, which = T],
pr_one = d[, mean(ua == 0)],
prior_only = 0,
sd_a0 = 1.5,
sd_b1 = 3, sd_b2 = 3, sd_b3 = 3, sd_b4 = 3
)
f1 <- m1$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 1 finished in 11.2 seconds.
Chain 2 finished in 11.5 seconds.
Both chains finished successfully.
Mean chain execution time: 11.3 seconds.
Total execution time: 11.6 seconds.Data simulation and parameter estimation
# snk <- capture.output(
# f1 <- m1$pathfinder(ld, num_paths=20, single_path_draws=200,
# history_size=50, max_lbfgs_iters=100,
# refresh = 0, draws = 2000)
# )
# f_mode <- m1$optimize(data = ld, jacobian = TRUE, refresh =0)
# f1 <- m1$laplace(data = ld, mode = f_mode, refresh = 0)
d_s <- d[, .(y = sum(y), n = .N), keyby = .(
a
)]
ld <- list(
N = nrow(d_s),
y = d_s$y,
n = d_s$n,
X = matrix(d_s$a),
prior_only = 0,
sd_a0 = 1.5,
sd_b = 3
)
f2 <- m2$sample(
ld, iter_warmup = 1000, iter_sampling = 1000,
parallel_chains = 2, chains = 2, refresh = 0, show_exceptions = F,
max_treedepth = 10)Running MCMC with 2 parallel chains...
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Both chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.1 seconds.Parameters estimated on the log-odds scale, but exponentiated below (intercept reported on risk scale).
Posterior distributions for model parameters
d_post <- data.table(f1$draws(variables = c("a0","b1","b2","b3","b4"), format = "matrix"))
names(d_post) <- c("a0", names(g_pars))
# parameters
cols <- c("a0", names(g_pars))
d_fig <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig[variable == "a0", value := plogis(value)]
d_fig[variable != "a0", value := exp(value)]
d_col <- unique(d_fig[, .(variable)])
d_col[, fill := factor(c(
1, 2, 2, 2,
1, 2,
1, 1,
2,
2, 1,
1, 1
))]
ggplot(d_fig, aes(x = value, group = variable)) +
geom_rect(data = d_col, aes(fill = fill),
xmin = -Inf,xmax = Inf,
ymin = -Inf,ymax = Inf, alpha = 0.1,
inherit.aes = F) +
scale_fill_manual("", values = c("white", "red")) +
geom_density() +
geom_vline(data = d_fig[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
scale_x_continuous(
labels = label_number(accuracy = 0.01)
) +
facet_wrap(~variable, scales = "free") +
theme(legend.position = "none")
The distribution of the marginal probability of \(Y\) under each treatment type.
Marginal probability of outcome by surgical type
cols <- c("marg_p0", "marg_p1")
d_post <- data.table(f1$draws(variables = cols, format = "matrix"))
names(d_post) <- cols
d_fig1 <- melt(d_post[, cols, with = F], measure.vars = cols)
d_fig2 <- data.table(f2$draws(variables = cols, format = "matrix"))
d_fig2 <- melt(d_fig2, measure.vars = cols)
ggplot(d_fig1, aes(x = value, group = variable)) +
geom_density() +
geom_density(data = d_fig2, aes(x = value, group = variable),
col = 2) +
geom_vline(data = d_fig1[, .(mu = mean(value)), keyby = variable],
aes(xintercept = mu, group = variable)) +
facet_wrap(~variable, ncol = 1) +
scale_x_continuous(
"Proportion with successful outcome",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
)
Risk difference (rev - dair). The red dashed line shows the posterior distribution obtained from a reduced model that includes a single term for surgical treatment assignment (dair/rev).
Risk difference (dair vs revision)
d_fig1 <- data.table(f1$draws(variables = c("rd"), format = "matrix"))
d_fig2 <- data.table(f2$draws(variables = c("rd"), format = "matrix"))
ggplot(d_fig1, aes(x = rd)) +
geom_density() +
scale_x_continuous(
"Risk difference",
# labels = label_number(accuracy = 0.001),
labels = scales::percent_format(accuracy = 0.1)
) +
geom_density(data = d_fig2,
aes(x = rd), col = 2, lty = 2) +
geom_vline(data = d_fig1[, .(mu = mean(rd))],
aes(xintercept = mu)) +
geom_vline(data = d_fig2[, .(mu = mean(rd))],
aes(xintercept = mu), col = 2) 
Discussion
In a factorial design, the main effects could be isolated and have an interpretation conditional on all other terms in the model held at a constant level. When we are talking about an effect of \(P\) it relates specifically to intervening on \(P\). However, due to the roadmap design using the marginal perspective could lead us to identify positive or negative revision effects solely as a consequence of the non-zero effects in other domains.
The phenomena is an artefact of the design rather than the estimation process or the specification of the statistical model.
The only reasonable way to address this seems to be to define the concept of revision as a strategy that involves multiple domain elements. Probably the most interpretaable approach is to have dair vs revision corresponding to dair being a strict strategy of dair plus 12 weeks backbone therapy and where revision conrresponds to one-stage revision with 12 weeks backbone antibiotic or two-stage revision with 12 weeks backbone antibiotic and no extended prophylaxis. However, given that we do not have an intervention arm that is the absence of extended prophylaxis, we would have to settle on 7 days.
In no way does the above fix the design, but it provides a more coherent quantity to evaluate and report on than an arbitrary mix of strategies.