Genentech / jmpost

https://genentech.github.io/jmpost/
17 stars 4 forks source link

Fix Random Slope Model Link Parameter #85

Closed gowerc closed 1 year ago

gowerc commented 1 year ago

Currently it looks like the link parameter to the random slope model isn't working. As in its credible intervals are massive even with reasonably large amounts of data. This may just be the nature of the model but it may also indicate a bug in the code. Needs to be investigated and resolved.

example code:

devtools::load_all(export_all = FALSE)

# Fully specified model
jm <- JointModel(
    longitudinal_model = LongitudinalRandomSlope(),
    link = LinkRandomSlope(),
    survival_model = SurvivalWeibullPH()
)

jlist <- simulate_joint_data(
    n = 1000,
    max_time = 2000,
    lambda_cen = 1 / 9000,
    beta_cat = c(
        "A" = 0,
        "B" = -0.1,
        "C" = 0.5
    ),
    beta_cont = 0.3,
    lm_fun = sim_lm_random_slope(
        intercept = 30,
        sigma = 3,
        slope = 2,
        z_sigma = 0.2,
        phi = 0.02
    ),
    os_fun = sim_os_weibull(
        lambda = 1 / 300,
        gamma = 0.97
    )
)

dat_os <- jlist$os
dat_lm <- jlist$lm |>
    dplyr::filter(time %in% c(1, 50, 100, 150, 200, 250, 300)) |>
    dplyr::arrange(time, pt)

stan_data <- as_stan_data(dat_os, dat_lm, ~ cov_cat + cov_cont)

dir.create(path = file.path("local"), showWarnings = FALSE)
mp <- sampleStanModel(
    jm,
    data = stan_data,
    iter_sampling = 1000,
    iter_warmup = 1000,
    chains = 1,
    parallel_chains = 1,
    exe_file = file.path("local", "full")
)
gowerc commented 1 year ago

Recommendations from @danielinteractive were to check the following:

gowerc commented 1 year ago

Phi with true value of 0.1

for 800 / 1600 / 4000 patients respectively:

variable       mean   median      sd     mad       q5      q95  rhat
link_lm_phi  0.132    0.133   3.49e-2 3.46e-2  0.0756   0.188   1.00
link_lm_phi  0.0881   0.0886  2.48e-2 2.46e-2  0.0458   0.127   1.01 
link_lm_phi  0.117    0.117   1.57e-2 1.49e-2  0.0898   0.142   1.00

Seems to be working exactly as expected now will close