Background
Often more than one regimen is studied in dose-finding studies. If there are enough doses within each regimen, one may still utilize MCP-Mod. But specific assumptions are needed, and it depends on the situation, whether or not these are appropriate (and thus usage of MCP-Mod).
The first idea is to bring the doses for each regimen on a common scale (total dose per time unit). For example if once daily (od) dosing and twice daily (bid) dosing are used in a study, one might utilize the total daily dose.
It is usually not appropriate to then perform MCP-Mod on the total daily dose (ignoring from which regimen the doses originate): The study investigated more than one regimen, so assessing the difference between regimen (for example for the same total daily dose) is of interest. This would not be possible with a modelling approach that ignores the regimen.
The most general approach would be to perform MCP-Mod separately by regimen, and for example to adjust p-values originating from the MCP-part using a Bonferroni correction. This approach assumes that the regimen don’t share any similarity. Due to the double-blind nature of trials, all patients would receive two administrations per day (patients in the od group receive one placebo per day), so that there is no real od group and in particular no separate placebo od group. So it often makes more sense to assume that the placebo group is common to both the od and bid dose-response curve. For the MCP-step contrasts for both od and bid are taken with respect to the same placebo group and in the modelling step one would assume the intercept to be the same across regimen, but all other parameters separate.
One could also assume further parameters to be common across regimen (for example the Emax or the ED50 parameter for the Emax model), but in the following example no such assumption is made.
The motivation for the simulated data below is taken from a recently completed dose-finding study, where the dose-response of the drug Licogliflozin was assessed for the od and bid regimen (Bays et al. 2020), see also the corresponding page at clinicaltrials.gov.
Note that this study used MCP-Mod, but the analysis presented here has been modified and simplified (in terms of candidate models and dose-response modelling strategy).
For most of the following code it is useful to structure the first-stage estimates like this: \[ \hat\mu=(\hat\mu_{\mathrm{placebo}}, \hat\mu_{\mathrm{od}}, \hat\mu_{\mathrm{bid}}) \] The length of the sub-vectors \(\hat\mu_{\mathrm{od}}\) and \(\hat\mu_{\mathrm{bid}}\) correspond to the number of different doses in the two regimens. They can be different, but in our example both have 4 elements.
Also as discussed above everything is modeled on the total daily dose scale.
library(DoseFinding)
library(ggplot2)
## collect estimates and dosage information in one place
example_estimates <- function() {
## ANOVA mean estimates and ci bounds extracted from fig. 3 of Bays (2020).
## clinicaltrials.gov page already seems to contain values from the dose-response model fit
mn <- c(-0.55, -1.78, -1.95, -3.29, -4.43, -1.14, -2.74, -4.03, -4.47)
lb <- c(-1.56, -3.15, -3.36, -4.85, -5.40, -2.49, -4.10, -5.50, -5.50)
ub <- c( 0.40, -0.30, -0.54, -1.76, -3.48, 0.24, -1.38, -2.65, -3.44)
se <- (ub - lb)/(2*qnorm(0.975)) # approximate standard error
return(list(mu_hat = mn,
daily_dose = c(0, 2.5, 10, 50, 150, 5, 10, 50, 100),
S_hat = diag(se^2),
# keep track of which elements correspond to which regimen:
index = list(placebo = 1, od = 2:5, bid = 6:9)))
}
## restructure estimates for easy plotting with ggplot
tidy_estimates <- function(est) {
se <- sqrt(diag(est$S_hat))
tidy <- data.frame(daily_dose = est$daily_dose, mu_hat = est$mu_hat,
ub = est$mu_hat + qnorm(0.975) * se, lb = est$mu_hat - qnorm(0.975) * se)
tidy <- rbind(tidy[1, ], tidy) # duplicate placebo
tidy$regimen <- c("od", "bid", rep("od", length(est$index$od)), rep("bid", length(est$index$bid)))
return(tidy)
}
plot_estimates <- function(est) {
df <- tidy_estimates(est)
ggplot(df, aes(daily_dose, mu_hat)) + geom_point() +
geom_errorbar(aes(ymin = lb, ymax = ub)) +
facet_wrap(vars(regimen), labeller = label_both) +
xlab("daily dose") + ylab("percent body weight cange") +
labs(title = "ANOVA estimates with 95% confindence intervals")
}
est <- example_estimates()
plot_estimates(est)
