In this vignette we illustrate how to use the
DoseFinding package with longitudinal, continuous
observations by fitting a mixed model for repeated measures (MMRM) and
applying the generalized MCP-Mod methodology to the resulting estimates.
We also show how to apply futility analyses during the trial, following
Bornkamp et al. (2024).
Data Simulation Code
In order to illustrate the methodology, we are going to consider data as simulated in Section 4.1 of Bornkamp et al. (2024). We include the corresponding functions here for completeness. This can be helpful for users who want to simulate their own data during trial design.
Calculate Mean Response from Emax Model
The calcMean function calculates the mean response at
each dose and each time point under the Emax model with \(\text{E}_0 = 0\) and \(\text{ED}_{50} = 1\) with the \(\text{E}_{max}\) parameter gradually
increasing so that a pre-specified maximum effect is achieved at the
last visit. The function takes three arguments: dose is a
vector of doses, times is a vector of time points, and
maxEff is the maximum effect at the final analysis (i.e.,
for the highest dose at the last visit). It returns a matrix of mean
responses, where the rows correspond to doses and the columns correspond
to time points.
calcMean <- function(dose, times, maxEff) {
emaxLastVisit <- Mods(emax = 1, doses = dose, maxEff = maxEff)$emax["eMax"]
maxTime <- max(times)
outer(dose, times, function(doses, times) {
Emax <- emaxLastVisit * (1 - exp(-0.5 * times)) / (1 - exp(-0.5 * maxTime))
Emax * doses / (doses + 1)
})
}
# Let's try it out:
calcMean(
dose = c(0, 0.5, 1, 2, 4),
times = 0:10,
maxEff = 0.1
) [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[2,] 0 0.01650577 0.02651703 0.03258916 0.03627210 0.03850591 0.03986079
[3,] 0 0.02475866 0.03977554 0.04888374 0.05440814 0.05775886 0.05979118
[4,] 0 0.03301154 0.05303405 0.06517832 0.07254419 0.07701182 0.07972157
[5,] 0 0.03961385 0.06364086 0.07821399 0.08705303 0.09241418 0.09566588
[,8] [,9] [,10] [,11]
[1,] 0.00000000 0.00000000 0.00000000 0.00000000
[2,] 0.04068256 0.04118099 0.04148330 0.04166667
[3,] 0.06102384 0.06177149 0.06222496 0.06250000
[4,] 0.08136512 0.08236198 0.08296661 0.08333333
[5,] 0.09763814 0.09883438 0.09955993 0.10000000Generate Normal Mean Vector and Covariance Matrix
The parGen function returns the mean and
(compound-symmetry) covariance matrix for the continuous outcomes, given
the same input as for calcMean above, plus the baseline
mean bslMean, and the standard deviation errSD
and constant correlation rho of the residuals:
parGen <- function(times, doses, maxEff, bslMean = 0, errSD = 0.55, rho = 0.9) {
times <- sort(times)
ndim <- length(times)
## matrix of outcome means for each dose and visit
MeanMat <- calcMean(doses, times, maxEff) + bslMean
rownames(MeanMat) <- paste0("Dose", doses)
colnames(MeanMat) <- paste0("Week", times)
## cov mat. for err
sdV <- if (length(errSD) == 1) {
rep(errSD, ndim)
} else if (length(errSD) == ndim) {
errSD
} else {
stop("Length of err sd should either be 1 (homogeneous) or equal to the number of visits in times")
}
## CS covariance structure
R <- diag(1, ndim)
R[lower.tri(R)] <- R[upper.tri(R)] <- rho
Sigm <- diag(sdV) %*% R %*% diag(sdV)
list(visit = times, Sigm = Sigm, MeanMat = MeanMat, bslMean = bslMean)
}
# Let's try it out:
pars <- parGen(
times = 0:10,
doses = c(0, 0.5, 1, 2, 4),
maxEff = 0.1,
bslMean = 1.5
)
pars$visit
[1] 0 1 2 3 4 5 6 7 8 9 10
$Sigm
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.30250 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225
[2,] 0.27225 0.30250 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225
[3,] 0.27225 0.27225 0.30250 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225
[4,] 0.27225 0.27225 0.27225 0.30250 0.27225 0.27225 0.27225 0.27225 0.27225
[5,] 0.27225 0.27225 0.27225 0.27225 0.30250 0.27225 0.27225 0.27225 0.27225
[6,] 0.27225 0.27225 0.27225 0.27225 0.27225 0.30250 0.27225 0.27225 0.27225
[7,] 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.30250 0.27225 0.27225
[8,] 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.30250 0.27225
[9,] 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.30250
[10,] 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225
[11,] 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225 0.27225
[,10] [,11]
[1,] 0.27225 0.27225
[2,] 0.27225 0.27225
[3,] 0.27225 0.27225
[4,] 0.27225 0.27225
[5,] 0.27225 0.27225
[6,] 0.27225 0.27225
[7,] 0.27225 0.27225
[8,] 0.27225 0.27225
[9,] 0.27225 0.27225
[10,] 0.30250 0.27225
[11,] 0.27225 0.30250
$MeanMat
Week0 Week1 Week2 Week3 Week4 Week5 Week6 Week7
Dose0 1.5 1.500000 1.500000 1.500000 1.500000 1.500000 1.500000 1.500000
Dose0.5 1.5 1.516506 1.526517 1.532589 1.536272 1.538506 1.539861 1.540683
Dose1 1.5 1.524759 1.539776 1.548884 1.554408 1.557759 1.559791 1.561024
Dose2 1.5 1.533012 1.553034 1.565178 1.572544 1.577012 1.579722 1.581365
Dose4 1.5 1.539614 1.563641 1.578214 1.587053 1.592414 1.595666 1.597638
Week8 Week9 Week10
Dose0 1.500000 1.500000 1.500000
Dose0.5 1.541181 1.541483 1.541667
Dose1 1.561771 1.562225 1.562500
Dose2 1.582362 1.582967 1.583333
Dose4 1.598834 1.599560 1.600000
$bslMean
[1] 1.5Note that we assume here that the time is in weeks, but that could of course easily be adapted as needed for your project.
Simulate Longitudinal Outcomes
The datGen function simulates the longitudinal outcome
data, given a sample size of patients N, a list of
parameters pars (as generated by parGen), a
vector of doses and corresponding allocation ratios
alRatio, the enrollment time of the last patient
LPFV.t. The distribution of the recruitment time can be
specified by RecrDist and RandRecr, with the
following options:
-
RecrDist = "Quad": quadratic, \(a t^2\), where \(a = N / LPFV.t^2\).-
RandRecr = TRUE: \(t \sim \sqrt{N \cdot \text{Unif}(0,1) / a}\). -
RandRecr = FALSE: \(t = \sqrt{i / a}\), where \(i = 1, \dotsc, N\).
-
-
Exp: exponential, with rate \(\lambda = -\log(0.9 / N) / LPFV.t\).-
RandRecr = TRUE: \(t \sim \text{Exp}(\lambda)\). -
RandRecr = FALSE: \(t = F_{\text{Exp}(\lambda)}^{-1}(i/N)\) for \(i = 1, \dotsc, N-1\); \(t = LPFV.t\) for \(i = N\).
-
-
Unif: uniform, \(t \sim \text{Unif}(0, LPFV.t)\).-
RandRecr = TRUE: \(t \sim \text{Unif}(0, LPFV.t)\). -
RandRecr = FALSE: \(t = i \cdot LPFV.t / N\) for \(i = 1, \dotsc, N\).
-
datGen <- function(N = 300, pars, doses, alRatio, LPFV.t, RecrDist = "Quad", RandRecr = TRUE) {
K <- length(doses)
nt <- length(pars$visit)
## generate list of dose arms
nblk <- ceiling(N / sum(alRatio)) # number of blocks
d <- unlist(lapply(1:nblk, function(i) sample(rep(1:K, alRatio), sum(alRatio))))
d <- d[1:N]
err <- mvtnorm::rmvnorm(N, mean = rep(0, nt), sigma = pars$Sigm) # error term
Y <- pars$MeanMat[d, ] + err # complete outcomes
## enrollment time (week)
randT <- if (RecrDist == "Quad") {
a <- N / LPFV.t^2
if (RandRecr) {
sqrt(N * runif(N) / a)
} else {
sqrt(c(1:N) / a)
}
} else if (RecrDist == "Exp") {
lamb <- -log(0.9 / N) / LPFV.t
if (RandRecr) {
rexp(N, rate = lamb)
} else {
c(qexp((1:(N - 1)) / N, lamb), LPFV.t)
}
} else if (RecrDist == "Unif") {
if (RandRecr) {
runif(N, 0, LPFV.t)
} else {
(1:N) * LPFV.t / N
}
}
trdat <- cbind(1:N, doses[d], randT, Y)
colnames(trdat)[1:3] <- c("SUBJID", "Dose", "enroll.time")
out <- tibble::as_tibble(trdat) |>
tidyr::gather("Visit", "response", -c("SUBJID", "Dose", "enroll.time")) |>
dplyr::arrange(SUBJID)
out$week <- rep(pars$visit, N)
out$cal.time <- out$enroll.time + out$week
bsl <- out |>
dplyr::filter(week <= 0) |>
dplyr::select("SUBJID", "response") |>
dplyr::rename(base = "response")
out <- bsl |>
dplyr::left_join(out, by = "SUBJID", multiple = "all") |>
dplyr::mutate(chg = response - base) |>
dplyr::select(- dplyr::all_of("response")) |>
dplyr::rename(response = "chg")
out |>
dplyr::mutate(USUBJID = paste0("PAT", SUBJID)) |>
dplyr::mutate_at("Visit", function(x) factor(x, levels = paste0("Week", pars$visit))) |>
dplyr::mutate_at("Dose", function(x) factor(x, levels = doses)) |>
dplyr::arrange(SUBJID, week)
}
# Let's try it out:
dat <- datGen(
N = 300,
pars = pars,
doses = c(0, 0.5, 1, 2, 4),
alRatio = c(1, 1, 1, 1, 1),
LPFV.t = 10,
RecrDist = "Quad",
RandRecr = TRUE
)
head(dat, 15)
[38;5;246m# A tibble: 15 × 9
[39m
SUBJID base Dose enroll.time Visit week cal.time response USUBJID
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[3m
[38;5;246m<int>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[38;5;250m 1
[39m 1 1.83 4 9.75 Week0 0 9.75 0 PAT1
[38;5;250m 2
[39m 1 1.83 4 9.75 Week1 1 10.8 0.198 PAT1
[38;5;250m 3
[39m 1 1.83 4 9.75 Week2 2 11.8 -
[31m0
[39m
[31m.
[39m
[31m103
[39m PAT1
[38;5;250m 4
[39m 1 1.83 4 9.75 Week3 3 12.8 0.270 PAT1
[38;5;250m 5
[39m 1 1.83 4 9.75 Week4 4 13.8 -
[31m0
[39m
[31m.
[39m
[31m150
[39m PAT1
[38;5;250m 6
[39m 1 1.83 4 9.75 Week5 5 14.8 0.183 PAT1
[38;5;250m 7
[39m 1 1.83 4 9.75 Week6 6 15.8 0.158 PAT1
[38;5;250m 8
[39m 1 1.83 4 9.75 Week7 7 16.8 0.454 PAT1
[38;5;250m 9
[39m 1 1.83 4 9.75 Week8 8 17.8 -
[31m0
[39m
[31m.
[39m
[31m0
[39m
[31m28
[4m1
[24m
[39m PAT1
[38;5;250m10
[39m 1 1.83 4 9.75 Week9 9 18.8 0.474 PAT1
[38;5;250m11
[39m 1 1.83 4 9.75 Week10 10 19.8 0.354 PAT1
[38;5;250m12
[39m 2 0.462 1 1.52 Week0 0 1.52 0 PAT2
[38;5;250m13
[39m 2 0.462 1 1.52 Week1 1 2.52 0.405 PAT2
[38;5;250m14
[39m 2 0.462 1 1.52 Week2 2 3.52 0.271 PAT2
[38;5;250m15
[39m 2 0.462 1 1.52 Week3 3 4.52 0.382 PAT2 So we see that the function generates a data.frame with
the subject ID, baseline value, dose, enrollment time, visit and week
plus calendar time for each longitudinal observation.
Cut Interim Data
The intDat function cuts the interim data when a certain
proportion pct of patients have their final outcome. It
returns a list with the interim analysis time int.time, the
cut data dat and the distribution of the last visit
dist (i.e., the number of patients with their last visit at
each time point).
intDat <- function(data, pct = 0.5) {
N <- max(data$SUBJID)
tmax <- max(data$week)
sorted_final_cal_time <- data |>
dplyr::filter(week == tmax) |>
dplyr::pull(cal.time) |>
sort()
N_required <- ceiling(pct * N)
t.cut <- sorted_final_cal_time[N_required]
out <- list(int.time = t.cut)
out$dat <- subset(data, cal.time <= t.cut)
out$dist <- out$dat |>
dplyr::group_by(SUBJID) |>
dplyr::summarize(lastvis = max(week)) |>
dplyr::group_by(lastvis) |>
dplyr::summarize(
n = dplyr::n(),
pct = (100 * dplyr::n() / N)
)
out
}
# Let's try it out:
int_dat <- intDat(dat, pct = 0.5)
int_dat$int.time
[1] 17.06187
$dat
[38;5;246m# A tibble: 2,998 × 9
[39m
SUBJID base Dose enroll.time Visit week cal.time response USUBJID
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<fct>
[39m
[23m
[3m
[38;5;246m<int>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[3m
[38;5;246m<chr>
[39m
[23m
[38;5;250m 1
[39m 1 1.83 4 9.75 Week0 0 9.75 0 PAT1
[38;5;250m 2
[39m 1 1.83 4 9.75 Week1 1 10.8 0.198 PAT1
[38;5;250m 3
[39m 1 1.83 4 9.75 Week2 2 11.8 -
[31m0
[39m
[31m.
[39m
[31m103
[39m PAT1
[38;5;250m 4
[39m 1 1.83 4 9.75 Week3 3 12.8 0.270 PAT1
[38;5;250m 5
[39m 1 1.83 4 9.75 Week4 4 13.8 -
[31m0
[39m
[31m.
[39m
[31m150
[39m PAT1
[38;5;250m 6
[39m 1 1.83 4 9.75 Week5 5 14.8 0.183 PAT1
[38;5;250m 7
[39m 1 1.83 4 9.75 Week6 6 15.8 0.158 PAT1
[38;5;250m 8
[39m 1 1.83 4 9.75 Week7 7 16.8 0.454 PAT1
[38;5;250m 9
[39m 2 0.462 1 1.52 Week0 0 1.52 0 PAT2
[38;5;250m10
[39m 2 0.462 1 1.52 Week1 1 2.52 0.405 PAT2
[38;5;246m# ℹ 2,988 more rows
[39m
$dist
[38;5;246m# A tibble: 4 × 3
[39m
lastvis n pct
[3m
[38;5;246m<int>
[39m
[23m
[3m
[38;5;246m<int>
[39m
[23m
[3m
[38;5;246m<dbl>
[39m
[23m
[38;5;250m1
[39m 7 48 16
[38;5;250m2
[39m 8 56 18.7
[38;5;250m3
[39m 9 46 15.3
[38;5;250m4
[39m 10 150 50 Longitudinal Data Analysis
Given a longitudinal dataset (at an interim analysis), we can perform two analyses:
- Completers analysis: Subset the data to patients who have reached the last visit, and then apply a simple linear model.
- Repeated measures analysis: Fit a mixed model for repeated measures (MMRM) to the data.
In both cases, we then calculate the least square means \(\widehat{\mu}_{0t}\) and the covariance matrix \(S_{0t}\) of the least square means, as well as the residual standard deviation \(\sigma\) of the model.
Let’s code these in two corresponding functions:
Completers Analysis
analyzeCompleters <- function(dat, eval.time = "Week12") {
complDat <- dat |>
dplyr::filter(Visit == eval.time)
fit <- lm(response ~ Dose + base, data = complDat)
lsm <- emmeans::emmeans(fit, ~ Dose)
summ <- subset(summary(lsm))
rowIDs <- as.numeric(rownames(summ))
list(
mu0t = summ$emmean,
sigma = sigma(fit),
S0t = vcov(lsm)
)
}
# Let's try it out:
resultCompleters <- analyzeCompleters(
dat = int_dat$dat,
eval.time = "Week10"
)
resultCompleters$mu0t
[1] -0.06011853 0.04853842 0.06463132 0.11936622 0.14242645
$sigma
[1] 0.2499916
$S0t
0 0.5 1 2 4
0 1.838760e-03 3.272396e-06 -4.632441e-06 4.028387e-07 6.215931e-07
0.5 3.272396e-06 2.248530e-03 -2.341013e-05 2.035753e-06 3.141233e-06
1 -4.632441e-06 -2.341013e-05 2.049133e-03 -2.881835e-06 -4.446765e-06
2 4.028387e-07 2.035753e-06 -2.881835e-06 2.500083e-03 3.866923e-07
4 6.215931e-07 3.141233e-06 -4.446765e-06 3.866923e-07 1.953591e-03Repeated Measures Analysis
analyzeRepeated <- function(dat, eval.time = "Week12") {
form <- "response ~ Visit + Dose + base + Dose*Visit + base*Visit"
postBaselineDat <- dat |>
dplyr::filter(week > 0) |>
droplevels()
fit <- mmrm::mmrm(as.formula(paste0(form, " + us(Visit | USUBJID)")), data = postBaselineDat)
lsm <- emmeans::emmeans(fit, ~ Dose + Visit)
summ <- subset(summary(lsm), Visit == eval.time)
rowIDs <- as.numeric(rownames(summ))
list(
mu0t = summ$emmean,
sigma = sqrt(diag(mmrm::VarCorr(fit)))[eval.time],
S0t = vcov(lsm)[rowIDs, rowIDs]
)
}
# Let's try it out:
resultRepeated <- analyzeRepeated(
dat = int_dat$dat,
eval.time = "Week10"
)
resultRepeated$mu0t
[1] -0.02818037 0.05291721 0.09861362 0.13468919 0.14456095
$sigma
Week10
0.2513171
$S0t
0 Week10 0.5 Week10 1 Week10 2 Week10
0 Week10 1.430501e-03 -1.818752e-06 1.529028e-06 -5.639547e-07
0.5 Week10 -1.818752e-06 1.626728e-03 -1.358336e-05 1.101611e-06
1 Week10 1.529028e-06 -1.358336e-05 1.539021e-03 -8.100511e-08
2 Week10 -5.639547e-07 1.101611e-06 -8.100511e-08 1.743068e-03
4 Week10 1.596990e-07 -6.294172e-07 7.022021e-07 -1.325705e-07
4 Week10
0 Week10 1.596990e-07
0.5 Week10 -6.294172e-07
1 Week10 7.022021e-07
2 Week10 -1.325705e-07
4 Week10 1.484844e-03Generalized MCP-Mod
Final Analysis
At the final analysis, we would like to use the generalized MCP-Mod
methodology to test the hypothesis of a dose-response relationship. This
is easily done by providing the least square means via resp
and their covariance matrix via S to the
MCTtest function, which will then perform the MCP-Mod
analysis. We just need to also specify the dose levels and the models we
want to test: Here we just use three models, it could be more in
practice.
doses <- c(0, 0.5, 1, 2, 4)
models <- Mods(
emax = 2,
sigEmax = c(0.5, 3),
quadratic = -0.2,
doses = doses
)
resultFinal <- analyzeRepeated(
dat = dat,
eval.time = "Week10"
)
testFinal <- MCTtest(
dose = doses,
resp = resultFinal$mu0t,
S = resultFinal$S0t,
models = models,
alternative = "one.sided",
type = "general",
critV = TRUE,
pVal = TRUE,
alpha = 0.025
)
testFinalMultiple Contrast Test
Contrasts:
emax sigEmax quadratic
0 -0.657 -0.788 -0.722
0.5 -0.271 -0.203 -0.223
1 -0.013 0.251 0.167
2 0.309 0.363 0.611
4 0.632 0.378 0.167
Contrast Correlation:
emax sigEmax quadratic
emax 1.000 0.921 0.827
sigEmax 0.921 1.000 0.941
quadratic 0.827 0.941 1.000
Multiple Contrast Test:
t-Stat adj-p
emax 3.987 <0.001
sigEmax 3.588 <0.001
quadratic 3.574 <0.001
Critical value: 2.18 (alpha = 0.025, one-sided)
testFinal$critV[1] 2.180088
attr(,"Calc")
[1] TRUE
testFinal$tStat emax sigEmax quadratic
3.986915 3.588199 3.574414
attr(,"pVal")
[1] 3.940677e-05 2.383234e-04 2.215375e-04
attr(testFinal$tStat, "pVal")[1] 3.940677e-05 2.383234e-04 2.215375e-04Futility Interim Analysis
Here we show how to apply the futility analysis during the trial, as
described in Bornkamp et al. (2024). The idea is that based on the
available interim data, we calculate the predictive or conditional power
to detect a significant dose-response relationship at the final
analysis. If this value is below a certain threshold, we may want to
stop the trial for futility. We can use the powMCTInterim
function to calculate these values. To do so, we need two additional
ingredients:
- The contrast matrix for the models we want to test, which is
provided via the argument
contMat. - The covariance matrix at the final analysis, denoted by \(S_{[0,1]}\) in Bornkamp et al. (2024) and provided via the argument
S_01here. Here we assume that this matrix is diagonal, with constant entries \(\sigma^2 / n\), where \(\sigma\) is the residual standard deviation of the model and \(n\) is the number of patients in the final analysis. - For the conditional power, we also need to specify the assumed dose
group means at the primary analysis time-point. This is provided via the
argument
mu_assumed. Typically, this is based on adding the observed placebo response rate to the originally planned treatment effect (before study start).
Let’s calculate these now for the example interim data from above:
w <- rep(1, length(doses))
contMat <- optContr(models = models, w = w)$contMat
n_final <- rep(60, length(doses))
S01 <- diag(resultRepeated$sigma^2 / n_final)
# Predictive power:
predPower <- powMCTInterim(
contMat = contMat,,
mu_0t = resultRepeated$mu0t,
S_0t = resultRepeated$S0t,
S_01 = S01,
alpha = 0.025,
type = "predictive",
alternative = "one.sided"
)
predPower[1] 0.9996943
attr(,"error")
[1] 1.504715e-07
attr(,"msg")
[1] "Normal Completion"
# Conditional power:
deltaAssumed <- pars$MeanMat[, "Week10"] - pars$bslMean
deltaAssumed <- deltaAssumed - deltaAssumed[1] # assumed treatment difference vs placebo
mu_assumed <- resultRepeated$mu0t[1] + deltaAssumed # to obtain group means add observed placebo response
condPower <- powMCTInterim(
contMat = contMat,
mu_0t = resultRepeated$mu0t,
S_0t = resultRepeated$S0t,
S_01 = S01,
mu_assumed = mu_assumed,
alpha = 0.025,
type = "conditional",
alternative = "one.sided"
)
condPower[1] 0.9978589
attr(,"error")
[1] 6.728031e-07
attr(,"msg")
[1] "Normal Completion"Simulation study
In practice, it will in most cases be important to run a simulation
study to evaluate the performance of the proposed combination of
generalized MCP-Mod with MMRM. In particular, in order to calibrate the
futility threshold for the interim analysis, it is important to
understand the power loss at the final analysis in relation to the
futility threshold. The following code may serve as the starting point
for such a simulation study, and can be adapted as needed. We include
here the argument documentation in roxygen2 format for
simplicity. The code can also be used to reproduce the simulation study
in Section 4.1 of Bornkamp et al. (2024).
##' @param times Vector assessment times
##' @param bslMean Mean at baseline
##' @param errSD Residual standard deviation (on absolute scale)
##' @param rho Correlation of measurements over time (within patient)
##' @param maxEff Maximum effect achieved at last time-point for the highest dose
##' @param RecrDist Recruitment distribution (see function datGen)
##' @param LPFV.t Time for the last patient first visit
##' @param models DoseFinding::Mods object
##' @param ia_timing Information time when IAs are performed (% of patients having last visit)
##' @param N Overall sample size
##' @param doses doses to use
##' @param alRatio allocation ratios to the different doses
##' @param eval.time Name of final
##' @param alpha Type 1 error for test
##' @param nSim Number of simulations
##' @param delta_assumed Vector of treatment effects assumed at the different doses
##' @return Data frame with all results
sim_one_trial <- function(times, bslMean, errSD, rho, maxEff,
RecrDist, LPFV.t, models, ia_timing = seq(0.3, 0.8, 0.05), N = 531,
doses = c(0, 0.5, 1, 2, 4, 8), alRatio = c(2, 1, 1, 1, 2, 2),
eval.time = "Week12", alpha = 0.025, nSim = 1000, delta_assumed) {
contMat <- optContr(models = models, w = alRatio / sum(alRatio))$contMat
pars <- parGen(times, doses,
maxEff = maxEff, bslMean = bslMean, errSD = errSD,
rho = rho
)
mydat <- datGen(N, pars, doses, alRatio, LPFV.t = LPFV.t, RecrDist = RecrDist, RandRecr = TRUE)
## fit final data
fit.final <- analyzeCompleters(mydat)
test <- MCTtest(
dose = doses, resp = fit.final$mu0t, S = fit.final$S0t, models = models,
alternative = "one.sided", type = "general", critV = TRUE, pVal = FALSE,
alpha = alpha
)
maxT <- max(test$tStat)
cVal <- test$critVal
n_final <- mydat |>
dplyr::filter(Visit == eval.time) |>
dplyr::group_by(Dose) |>
dplyr::count()
IAtm <- pp_long <- cp_long <- cp_long2 <- pp_compl <- cp_compl <- cp_compl2 <- inf_long <- inf_compl <- numeric()
IAdist <- matrix(nrow = length(times), ncol = length(ia_timing))
## fit interim data
for (ia in seq_along(ia_timing)) {
tp <- ia_timing[ia] # prop pts completes visit of eval.time
cutDat <- intDat(mydat, tp)
middat <- cutDat$dat
IAtm[ia] <- cutDat$int.time
IAdist[times %in% cutDat$dist$lastvis, ia] <- cutDat$dist$pct
## using all data from every patient (longitudinal MMRM)
fit_int_long <- analyzeRepeated(middat, eval.time = eval.time)
## only using the completers by visit eval time (cross-sectional linear model)
fit_int_compl <- analyzeCompleters(middat, eval.time = eval.time)
## The covariance matrix anticipated for the estimates at EoS
S_end <- diag(fit_int_long$sigma^2 / n_final$n)
# information fraction
inf_long[ia] <- det(diag(fit_int_long$sigma^2 / n_final$n))^(1 / length(doses)) /
det(fit_int_long$S0t)^(1 / length(doses))
inf_compl[ia] <- det(diag(fit_int_compl$sigma^2 / n_final$n))^(1 / length(doses)) /
det(fit_int_compl$S0t)^(1 / length(doses))
mu_assumed <- fit_int_long$mu0t[1] + delta_assumed # for conditional power use planned treatment effect
pp_long[ia] <- powMCTInterim(
contMat = contMat, S_0t = fit_int_long$S0t,
S_01 = S_end, mu_0t = fit_int_long$mu0t, type = "predictive", alpha = alpha
)
cp_long[ia] <- powMCTInterim(
contMat = contMat, S_0t = fit_int_long$S0t,
S_01 = S_end, mu_0t = fit_int_long$mu0t, type = "conditional", alpha = alpha,
mu_assumed = mu_assumed
)
cp_long2[ia] <- powMCTInterim(
contMat = contMat, S_0t = fit_int_long$S0t,
S_01 = S_end, mu_0t = fit_int_long$mu0t, type = "conditional", alpha = alpha,
mu_assumed = fit_int_long$mu0t
)
pp_compl[ia] <- powMCTInterim(
contMat = contMat, S_0t = fit_int_compl$S0t,
S_01 = S_end, mu_0t = fit_int_compl$mu0t, type = "predictive", alpha = alpha
)
cp_compl[ia] <- powMCTInterim(
contMat = contMat, S_0t = fit_int_compl$S0t,
S_01 = S_end, mu_0t = fit_int_compl$mu0t, type = "conditional", alpha = alpha,
mu_assumed = mu_assumed
)
cp_compl2[ia] <- powMCTInterim(
contMat = contMat, S_0t = fit_int_compl$S0t,
S_01 = S_end, mu_0t = fit_int_compl$mu0t, type = "conditional", alpha = alpha,
mu_assumed = fit_int_compl$mu0t
)
}
data.frame(
final_maxT = maxT, final_cVal = as.numeric(cVal),
ia_timing, pp_long, cp_long, cp_long2, pp_compl, cp_compl, cp_compl2, inf_long, inf_compl
)
}
## Function to run one scenario (arguments described in previous function)
run_scen <- function(n_sim, rho, maxEff,
LPFV.t, ia_timing, delta_assumed) {
## fixed parameters
doses <- c(0, 0.5, 1, 2, 4, 8) # doses
alRatio <- c(2, 1, 1, 1, 2, 2) # allocation ratio for doses
alpha <- 0.025
times <- c(0, 2, 4, 8, 12)
bslMean <- 0
errSD <- 0.56
RecrDist <- "Quad"
if (rho == 0.6 & LPFV.t == 50) N <- 820
if (rho == 0.6 & LPFV.t == 100) N <- 825
if (rho == 0.9 & LPFV.t == 50) N <- 248
if (rho == 0.9 & LPFV.t == 100) N <- 236
models <- Mods(
emax = c(0.5, 1, 2, 4), sigEmax = rbind(c(0.5, 3), c(1, 3), c(2, 3), c(4, 3)),
quadratic = -0.1, doses = doses
)
lst <- vector("list", n_sim)
for (i in 1:n_sim) {
res <- try(sim_one_trial(
times = times, bslMean = bslMean,
errSD = errSD, rho = rho, maxEff = maxEff, RecrDist = RecrDist,
LPFV.t = LPFV.t, models = models, ia_timing = ia_timing, N = N,
doses = doses, alRatio = alRatio, delta_assumed = delta_assumed,
alpha = alpha
), silent = FALSE)
if (!is.character(res)) {
lst[[i]] <- data.frame(sim = i, res)
}
}
out <- do.call("rbind", lst)
out$rho <- rho
out$maxEff <- maxEff
out$RecrDist <- RecrDist
out$N <- N
out$LPFV.t <- LPFV.t
out
}These functions can be used as follows:
doses <- c(0, 0.5, 1, 2, 4, 8) # doses
alRatio <- c(2, 1, 1, 1, 2, 2) # allocation ratio for doses
alpha <- 0.025 # for MCTtest
times <- c(0, 2, 4, 8, 12) # observation times (weeks)
bslMean <- 0 # mean at baseline
errSD <- 0.56 # SD for outcome on absolute scale
RecrDist <- "Quad" # recruitment distribution
rho <- 0.9 # correlation across time for outcome on absolute scale
LPFV.t <- 100 # time of last patient first visit
maxEff <- 0.12 # effect assumed for the highest dose at the last time point
delta_assumed <- calcMean(doses, times, maxEff)[, length(times)]
N <- 236
ia_timing <- 0.5
models <- Mods(
emax = c(0.5, 1, 2, 4), sigEmax = rbind(c(0.5, 3), c(1, 3), c(2, 3), c(4, 3)),
quadratic = -0.1, doses = doses
)
oneSim <- sim_one_trial(
times = times, bslMean = bslMean,
errSD = errSD, rho = rho, maxEff = maxEff, RecrDist = RecrDist,
LPFV.t = LPFV.t, models = models, ia_timing = ia_timing, N = N,
doses = doses, alRatio = alRatio, delta_assumed = delta_assumed,
alpha = alpha
)
scenRes <- run_scen(
n_sim = 10, rho = 0.9, maxEff = 0.12,
LPFV.t = 50, ia_timing = seq(0.3, 0.8, 0.05), delta_assumed = delta_assumed
)References
Bornkamp, B., Zhou, J., Xi, D., and Cao, W. (2024), “Futility analyses for the
MCP-mod methodology based on longitudinal models.”