In pavlakrotka/NCC: Simulation and Analysis of Platform Trials with Non-Concurrent Controls

library(knitr)
opts_chunk$set(echo = TRUE,
               message = FALSE,
               error = FALSE,
               warning = FALSE,
               comment = "",
               fig.align = "center",
               out.width = "100%")

options(mc.cores=3)

library(NCC)
library(ggplot2)

Preparing scenarios

To perform a simulation study with the NCC package, we first create a data frame with the desired scenarios that contains all the parameters needed for data generation and analysis.

In this simple example, we consider a platform trial with 4 experimental treatment arms entering sequentially, where the null hypothesis holds for all experimental arms. We vary the strength and pattern of the time trend (parameters lambda and trend) in order to investigate their impact on the type I error, bias and MSE of the treatment effect estimates.

sim_scenarios <- data.frame(num_arms = 4, 
                            n_arm = 250, 
                            d1 = 250*0,
                            d2 = 250*1,
                            d3 = 250*2,
                            d4 = 250*3,
                            period_blocks = 2, 
                            mu0 = 0,
                            sigma = 1,
                            theta1 = 0,
                            theta2 = 0,
                            theta3 = 0,
                            theta4 = 0,
                            lambda0 = rep(seq(-0.15, 0.15, length.out = 9), 2),
                            lambda1 = rep(seq(-0.15, 0.15, length.out = 9), 2),
                            lambda2 = rep(seq(-0.15, 0.15, length.out = 9), 2),
                            lambda3 = rep(seq(-0.15, 0.15, length.out = 9), 2),
                            lambda4 = rep(seq(-0.15, 0.15, length.out = 9), 2),
                            trend = c(rep("linear", 9), rep("stepwise_2", 9)),
                            alpha = 0.025,
                            ncc = TRUE)

head(sim_scenarios)

Running simulations

We use the function sim_study_par() to perform a simulation study with the created scenarios. Here we evaluate the 4th experimental treatment arm using the regression model with period adjustment, as well as the separate and pooled analyses. Each scenario will be replicated 1000 times.

set.seed(1234)
sim_results <- sim_study_par(nsim = 1000, scenarios = sim_scenarios, arms = 4, 
                             models = c("fixmodel", "sepmodel", "poolmodel"), 
                             endpoint = "cont")

[1] "Starting the simulations. 18 scenarios will be simulated. Starting time: 2023-02-19 14:24:09"  
[1] "Scenario 1/18 done. Time: 2023-02-19 14:24:21"
[1] "Scenario 2/18 done. Time: 2023-02-19 14:24:26"  
[1] "Scenario 3/18 done. Time: 2023-02-19 14:24:32"  
[1] "Scenario 4/18 done. Time: 2023-02-19 14:24:37"
[1] "Scenario 5/18 done. Time: 2023-02-19 14:24:42"
[1] "Scenario 6/18 done. Time: 2023-02-19 14:24:47"
[1] "Scenario 7/18 done. Time: 2023-02-19 14:24:53"
[1] "Scenario 8/18 done. Time: 2023-02-19 14:24:58"
[1] "Scenario 9/18 done. Time: 2023-02-19 14:25:03"
[1] "Scenario 10/18 done. Time: 2023-02-19 14:25:08"
[1] "Scenario 11/18 done. Time: 2023-02-19 14:25:13"
[1] "Scenario 12/18 done. Time: 2023-02-19 14:25:19"
[1] "Scenario 13/18 done. Time: 2023-02-19 14:25:24"
[1] "Scenario 14/18 done. Time: 2023-02-19 14:25:30"
[1] "Scenario 15/18 done. Time: 2023-02-19 14:25:36"
[1] "Scenario 16/18 done. Time: 2023-02-19 14:25:41"
[1] "Scenario 17/18 done. Time: 2023-02-19 14:25:47"
[1] "Scenario 18/18 done. Time: 2023-02-19 14:25:52"

The function reports the system time after each scenario finishes in order to track the progress of the simulations.

Simulation results

# Add results manually because of long runtime of the simulation (CRAN check)

sim_results <- sim_scenarios[rep(seq_len(nrow(sim_scenarios)), each = 3), ]
rownames(sim_results) <- c(1:nrow(sim_results))

sim_results$study_arm <- 4
sim_results$model <- rep(c("fixmodel", "poolmodel", "sepmodel"), 18)

sim_results$reject_h0 <- c(0.030, 0.007, 0.027, 0.020, 0.008, 0.020, 0.025, 0.009, 0.018, 0.023, 0.018, 0.032, 0.031, 0.023, 0.027, 0.029, 0.033, 0.031, 0.031, 0.049, 0.032, 0.024, 0.053, 0.022, 0.027, 0.093, 0.025, 0.031, 0.000, 0.020, 0.020, 0.000, 0.022, 0.034, 0.002, 0.035, 0.013, 0.003, 0.018, 0.021, 0.020, 0.015, 0.025, 0.082, 0.021, 0.024, 0.199, 0.029, 0.028, 0.380, 0.029, 0.021, 0.617, 0.019)

sim_results$bias <- c(-0.0005506988, -0.0438708399, 0.0005236461, -0.0001912505, -0.0332907139, -0.0002745666, 0.0014866433, -0.0224244517, 0.0014689339, 0.0019243764, -0.0132531632, 0.0012354448, -0.0008100829, -0.0006305263, 0.0005062764, 0.0007936728, 0.0105594167, 0.0007603154, -0.0029095179, 0.0177743581, -0.0020748971, -0.0004732032, 0.0320234632, -0.0003444024, -0.0018542583, 0.0441438285, -0.0020824008, 0.0041050797, -0.1684501438, 0.0035396805, -0.0052914700, -0.1338118193, -0.0061579896, 0.0007975244, -0.0861302426, -0.0017007487, -0.0008201267, -0.0447944941, -0.0021701199, -0.0021257738, -0.0007531575, -0.0026113360, -0.0021726726,  0.0429857395, -0.0039479999, 0.0011433884, 0.0888110442, 0.0028729606, -0.0011313254, 0.1288540403, -0.0005165626, -0.0003261784, 0.1735776879, 0.0009341703)

sim_results$MSE <- c(0.007717052, 0.008220890, 0.008487647, 0.006601056, 0.006741297, 0.007570587, 0.007080945, 0.006059043, 0.007976836, 0.007338159, 0.006113367, 0.008393631, 0.007527355, 0.006224911, 0.008203350, 0.007192991, 0.006072723, 0.008228014, 0.008187574, 0.006498952, 0.009307924, 0.007103126, 0.006636311, 0.007794670, 0.007697291, 0.008650210, 0.008328635, 0.007325503, 0.034469812, 0.007933315, 0.007306962, 0.023967797, 0.008274504, 0.007618799, 0.013630920, 0.008346264, 0.006848033, 0.007418390, 0.007645439, 0.006918091, 0.005676356, 0.007640708, 0.007514734, 0.007862916, 0.008116238, 0.007187846, 0.013994993, 0.008139475, 0.006881168, 0.022324478, 0.007473595, 0.006596491, 0.035255110, 0.007170597)

sim_results$failed <- 0
sim_results$nsim <- 1000

The resulting data frame contains the considered scenarios, as well as simulation results - probability to reject the null hypothesis, bias and MSE of the treatment effect estimates.

head(sim_results)

We can now visualize the performance of the considered analysis methods with respect to the strength and pattern of the time trend.

Type I error

ggplot(sim_results, aes(x=lambda0, y=reject_h0, color=model)) +
  geom_point() +
  geom_line() +
  facet_grid(~ trend) +
  geom_hline(aes(yintercept = 0.025), linetype = "dotted") +
  labs(x="Strength of time trend", y="Type I error", color="Analysis approach") +
  theme_bw()

Bias

ggplot(sim_results, aes(x=lambda0, y=bias, color=model)) +
  geom_point() +
  geom_line() +
  facet_grid(~ trend) +
  geom_hline(aes(yintercept = 0), linetype = "dotted") +
  labs(x="Strength of time trend", y="Bias", color="Analysis approach") +
  theme_bw()

MSE

ggplot(sim_results, aes(x=lambda0, y=MSE, color=model)) +
  geom_point() +
  geom_line() +
  facet_grid(~ trend) +
  labs(x="Strength of time trend", y="MSE", color="Analysis approach") +
  theme_bw()

pavlakrotka/NCC documentation built on April 17, 2025, 3:11 a.m.