In this vignette we will illustrate the usage of the DoseFinding package for analyzing continuously distributed data. There is a separate vignette with details on sample size and power calculation.
We will use data from Verkindre et al. (2010), who actually use a cross-over design and utilize MCP-Mod in a supportive analysis. More information can be found at the corresponding clinicaltrials.gov page and on the R help page ?glycobrom.
The main purpose Verkindre et al. (2010) was to provide safety and efficacy data on Glycopyrronium Bromide (NVA237) in patients with stable Chronic Obstructive Pulmonary Disease (COPD). The primary endpoint in this study was the mean of two measurements of forced expiratory volume in 1 second (FEV1) at 23h 15min and 23h 45min post dosing, following 7 days of treatment.
In order to keep this exposition simple, we will ignore the active control and focus on the placebo group and the four dose groups (12.5, 25, 50, and 100μg).
For the purpose here, we recreate a dataset that mimicks a parallel group design, based on the published summary statistics. These can be found in the glycobrom dataset coming with the DoseFinding package. Here fev1 and sdev contain the mean and standard deviation of the mean (standard error) of the primary endpoint for each group, while n denotes the number of participants.
library(DoseFinding)
data(glycobrom)
print(glycobrom) dose fev1 sdev n
1 0.0 1.243 0.0156 49
2 12.5 1.317 0.0145 55
3 25.0 1.333 0.0151 51
4 50.0 1.374 0.0148 53
5 100.0 1.385 0.0148 53We want to create a dataset with 60 participants in each of the five groups. Noticing that the standard errors are essentially equal across all groups, we draw five vectors of measurement errors centered at 0 with identical variances 60 * 0.015^2 which we add to the observed means.
Note that here we use MASS::mvrnorm instead of rnorm because it lets us generate random numbers with the specified sample mean and sd.
set.seed(1, kind = "Mersenne-Twister", sample.kind = "Rejection", normal.kind = "Inversion")
rand <- rep(MASS::mvrnorm(60, 0, 60 * 0.015^2, empirical = TRUE), 5)
NVA <- data.frame(dose = rep(glycobrom$dose, each = 60),
FEV1 = rep(glycobrom$fev1, each = 60) + rand)
ggplot(NVA) + geom_jitter(aes(dose, FEV1), height = 0, width = 4) +
labs(title = "Simulated FEV1 by dose (jittered horizontally)") +
xlab("dose [μg]") + ylab("FEV1 [l]")