The jmpost package includes data simulation functionality for the included joint models. The data simulation is based on specifying the longitudinal model and the survival model including link parameters.
Example
set.seed(129)
sim_data <- SimJointData(
design = list(
SimGroup(50, "Arm-A", "Study-X"),
SimGroup(50, "Arm-B", "Study-X")
),
longitudinal = SimLongitudinalSteinFojo(
times = c(-2, 1, 30, 90, 150, 210),
mu_g = log(c(0.005, 0.005)),
mu_s = log(c(0.06, 0.007)),
omega_g = 0.3,
link_dsld = 0.1,
),
survival = SimSurvivalWeibullPH(
lambda = 1 / 300,
gamma = 0.97,
time_max = 2000,
time_step = 1,
lambda_cen = 1 / 9000,
beta_cat = c(
"A" = 0,
"B" = -0.1,
"C" = 0.5
),
beta_cont = 0.3
)
)This object has survival and longitudinal components.
head(sim_data@survival)
#> # A tibble: 6 × 7
#> subject study arm time cov_cont cov_cat event
#> <chr> <fct> <fct> <dbl> <dbl> <fct> <dbl>
#> 1 subject_001 Study-X Arm-A 42 -1.12 B 1
#> 2 subject_002 Study-X Arm-A 20 -0.990 C 1
#> 3 subject_003 Study-X Arm-A 832 -1.37 C 1
#> 4 subject_004 Study-X Arm-A 112 -1.36 C 1
#> 5 subject_005 Study-X Arm-A 13 2.00 B 1
#> 6 subject_006 Study-X Arm-A 135 0.696 B 1
head(sim_data@longitudinal)
#> # A tibble: 6 × 6
#> subject arm study time sld observed
#> <chr> <fct> <fct> <dbl> <dbl> <lgl>
#> 1 subject_001 Arm-A Study-X -2 64.0 TRUE
#> 2 subject_001 Arm-A Study-X 1 61.9 TRUE
#> 3 subject_001 Arm-A Study-X 30 27.7 TRUE
#> 4 subject_001 Arm-A Study-X 90 58.9 FALSE
#> 5 subject_001 Arm-A Study-X 150 124. FALSE
#> 6 subject_001 Arm-A Study-X 210 226. FALSEWe can see the trajectory of the tumour size.
library(ggplot2)
ggplot(
sim_data@longitudinal,
aes(x = time, y = sld, group = subject)
) +
geom_line(alpha = 0.1) +
geom_point() +
facet_wrap(~arm)
We can also visualise the Kaplan-Meier survival curves.
Administrative Cut Off
We may wish to have a fixed maximum study duration and remove any simulated values after that time. This can be specified as a single fixed value or as a vector if we additionally wish to consider enrolment time.
Progression
We might also be interested in time to progression defined as an increase in tumour size over a threshold. The threshold is based on the minimum observed SLD up to the given time and the relative and absolute growth. For example we might require the SLD to increase at least 20% and at least 5mm. We don’t include any observations before time 0 in the calculation of the minimum and we don’t observe any SLD values after progression.
sim_data_pd <- add_pfs(
sim_data_end,
relative_threshold = 1.2,
absolute_threshold = 5,
from_time = 0,
observed_after = FALSE
)Let’s look at the plots again up to progression.
ggplot(
sim_data_pd@longitudinal |> dplyr::filter(observed), # here we only include observed values
aes(x = time, y = sld, group = subject)
) +
geom_line(alpha = 0.1) +
geom_point() +
facet_wrap(~arm)
Model Fitting from Simulated Data
As already descried in the quick start vignette, we can use this data to fit a joint model.
os_data <- sim_data_end@survival
long_data <- sim_data_end@longitudinal
joint_data <- DataJoint(
subject = DataSubject(
data = os_data,
subject = "subject",
arm = "arm",
study = "study"
),
survival = DataSurvival(
data = os_data,
formula = Surv(time, event) ~ cov_cat + cov_cont
),
longitudinal = DataLongitudinal(
data = long_data,
formula = sld ~ time,
threshold = 5
)
)
sf_model <- JointModel(
longitudinal = LongitudinalSteinFojo(
mu_bsld = prior_normal(log(60), 0.2),
mu_ks = prior_normal(-3, 0.4),
mu_kg = prior_normal(log(0.005), 0.3),
centred = TRUE
),
survival = SurvivalWeibullPH(),
link = linkDSLD()
)
set.seed(202671)
mcmc_results <- sampleStanModel(
sf_model,
data = joint_data,
iter_sampling = 1000,
iter_warmup = 500,
chains = 4,
parallel_chains = 4,
step_size = 0.01
)
knitr::kable(
mcmc_results@results$summary(
c("lm_sf_mu_bsld", "lm_sf_mu_ks", "lm_sf_mu_kg", "lm_sf_omega_bsld", "lm_sf_omega_ks", "lm_sf_omega_ks")
),
digits = 3
)| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| lm_sf_mu_bsld[1] | 4.079 | 4.079 | 0.020 | 0.020 | 4.046 | 4.113 | 1.003 | 654.772 | 1048.022 |
| lm_sf_mu_ks[1] | -2.866 | -2.867 | 0.026 | 0.026 | -2.909 | -2.823 | 1.005 | 538.208 | 1092.617 |
| lm_sf_mu_ks[2] | -4.976 | -4.977 | 0.030 | 0.030 | -5.024 | -4.924 | 1.005 | 508.038 | 674.176 |
| lm_sf_mu_kg[1] | -5.375 | -5.374 | 0.045 | 0.044 | -5.449 | -5.302 | 1.028 | 339.787 | 455.300 |
| lm_sf_mu_kg[2] | -5.344 | -5.344 | 0.044 | 0.043 | -5.417 | -5.271 | 1.008 | 494.226 | 758.850 |
| lm_sf_omega_bsld[1] | 0.208 | 0.208 | 0.015 | 0.014 | 0.186 | 0.234 | 1.005 | 418.249 | 833.851 |
| lm_sf_omega_ks[1] | 0.183 | 0.181 | 0.019 | 0.017 | 0.155 | 0.218 | 1.011 | 306.336 | 505.912 |
| lm_sf_omega_ks[2] | 0.211 | 0.209 | 0.022 | 0.020 | 0.179 | 0.248 | 1.012 | 296.813 | 476.726 |
Simulating from a Fitted Model
The package also have the functionality to generate new data based on the fitted model. We use the [simulate.JointModelSamples] function.
For this simulation we set
- the times we wish to have longitudinal observations
- variance for jitter around those observation times
- time_max and time_step for numerical integration
- lambda parameter for exponential censoring times
- scaled_variance which should be set to the same as the fitted model. This is described more in the statistical specifications vignette.
Now we can inspect the new simulated data.
head(new_model_data@longitudinal)
#> subject arm study time sld observed
#> 1 subject_001 Arm-A Study-X -2.000000 51.55855 TRUE
#> 2 subject_001 Arm-A Study-X 1.092538 49.50418 TRUE
#> 3 subject_001 Arm-A Study-X 31.020617 17.69610 TRUE
#> 4 subject_001 Arm-A Study-X 60.795321 14.73317 TRUE
#> 5 subject_001 Arm-A Study-X 91.002436 19.34256 TRUE
#> 6 subject_001 Arm-A Study-X 120.535991 26.21532 TRUE
head(new_model_data@survival)
#> # A tibble: 6 × 9
#> subject study arm time cov_cat.x cov_cont.x event cov_cat.y cov_cont.y
#> <fct> <fct> <fct> <dbl> <fct> <dbl> <dbl> <fct> <dbl>
#> 1 subject_001 Study… Arm-A 545 B -1.12 1 B -1.12
#> 2 subject_002 Study… Arm-A 316 C -0.990 1 C -0.990
#> 3 subject_003 Study… Arm-A 220 C -1.37 1 C -1.37
#> 4 subject_004 Study… Arm-A 574 C -1.36 1 C -1.36
#> 5 subject_005 Study… Arm-A 25 B 2.00 1 B 2.00
#> 6 subject_006 Study… Arm-A 102 B 0.696 1 B 0.696
ggplot(
dplyr::bind_rows(
"model sim" = new_model_data@longitudinal,
"original data" = sim_data@longitudinal,
.id = "sim"
),
aes(x = time, y = sld, group = subject)
) +
geom_line(alpha = 0.1) +
geom_point() +
facet_wrap(~ arm + sim) +
coord_cartesian(ylim = c(0, 400)) +
ggtitle("Original Data versus Model Simulated Data")
plot(
survfit(Surv(time, event) ~ arm, data = new_model_data@survival),
col = 1:2,
lwd = 2,
main = "Overall Survival (simulated from model)"
)
lines(
survfit(Surv(time, event) ~ arm, data = sim_data@survival),
col = 1:2,
lwd = 2,
lty = 2
)
legend("topright",
col = c(1, 2, 1, 2), lwd = 2, lty = c(1, 1, 2, 2),
legend = c("Arm-A from model", "Arm-B from model", "Arm-A original", "Arm-B original")
)