search

drizopoulos / JMbayes2

Extended Joint Models for Longitudinal and Survival Data
https://drizopoulos.github.io/JMbayes2/
78 stars 22 forks source link

JMBayes2 results don't seem to match JMBayes results #1

Closed jpkrooney closed 3 years ago

jpkrooney commented 3 years ago

Hi Dimirtris, Great to see this new package. I was just testing it out and I found that results from JMBayes and JMBayes2 don't match when the same sub-models are used as arguments.

A reprex: library(JMbayes) lmeFit.aids <- lme(CD4 ~ obstime * drug, random = ~ obstime | patient, data = aids) survFit.aids <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE) JMB1_model <- jointModelBayes(lmeFit.aids, survFit.aids, timeVar = "obstime") library(JMbayes2) JMB2_model <- jm(survFit.aids, lmeFit.aids, time_var = "obstime")

Looking at the results:

Screenshot 2021-01-14 at 13 29 10 Screenshot 2021-01-14 at 13 33 47

As you can see the results are quite different. I note that some that there are less parameters in JMBayes2, so perhaps I should not expect the same results ? However in my own data the AssocT parameter changed from a significant negative number to a significant positive number so I thought this issue was worth raising. Many thanks!

harryparr commented 3 years ago

Hi - I'm also having similar issues that my alpha value (and other event) parameters are in opposing directions between JMBayes2 and JMBayes::jointModelBayes & mvjointModelBayes

drizopoulos commented 3 years ago

JMbayes2 implements a better MCMC algorithm than JMbayes. In particular, it uses the hierarchical centering parameterization and also center covariates in both submodels. In addition, JMbayes2 uses three chains compared to the single chain of JMbayes. In that way, we can assess convergence better using, e.g., the German-Rubin diagnostic (i.e., Rhat values in the results). In this example, the Rhat values indicate successful convergence. Hence, all in all, the results from JMbayes2 are more trustworthy.

drizopoulos commented 3 years ago

But in any case, thanks for reporting. We will investigate this further...

jpkrooney commented 3 years ago

Thanks Dimitris that is great. The updates in JMB2 are great and so far it seems much faster to me!

Maybe it helpful to give the issue some real world context ?:

In ALS most (2/3) die from respiratory failure and in this paper (https://hrbopenresearch.org/articles/2-23 ) we are using a longitudinal metric of respiratory function such that lower scores => worse respiratory function. I reran one of the models from the paper in original JMBayes and in JMBayes2:

lme_ref <- lme(SNIP.occ ~ ns(onset2snip_mnths, 2) ,
               random = ~ ns(onset2snip_mnths, 2)  | ID,
               data = df.snip_72m)

coxC9 <- coxph(Surv(time = dx_delay, time2 = surv_t2_months, fail) ~ age_on +
                   dx_delay2 + simp_site + C9,
               data = df.clin_72m, x = TRUE, model = TRUE)

JMB1_mod <- jointModelBayes(lme_ref, coxC9, timeVar = "onset2snip_mnths",
                             verbose = FALSE)
library(JMbayes2)
JMB2_mod <- jm(coxC9, lme_ref, time_var = "onset2snip_mnths")

Here are the Event process results from both:

Screenshot 2021-01-14 at 17 45 07

As you can see the AssocT term has flipped sign in the JMB2 model (in addition to other changes). So this would then imply that decreased respiratory function is associated with improved survival -> which doesn't make sense given background knowledge.

drizopoulos commented 3 years ago

Thanks for the additional info. Would it be possible to try fitting with jm() again by altering the prior for the tau_bs_gammas. Namely using the following

jm(..., priors = list(A_tau_bs_gammas = 1, B_tau_bs_gammas = 0.01))

—— Dimitris Rizopoulos Professor of Biostatistics Erasmus University Medical Center The Netherlands

From: jpkrooney notifications@github.com Sent: Thursday, January 14, 2021 5:50:03 PM To: drizopoulos/JMbayes2 JMbayes2@noreply.github.com Cc: D. Rizopoulos d.rizopoulos@erasmusmc.nl; Comment comment@noreply.github.com Subject: Re: [drizopoulos/JMbayes2] JMBayes2 results don't seem to match JMBayes results (#1)

Thanks Dimitris that is great. The updates in JMB2 are great and so far it seems much faster to me!

Maybe it helpful to give the issue some real world context ?:

In ALS most (2/3) die from respiratory failure and in this paper (https://hrbopenresearch.org/articles/2-23https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fhrbopenresearch.org%2Farticles%2F2-23&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cc4afac1c59984532994308d8b8ac712d%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462398058107199%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=ono0G%2FToTjiMmxsW5uocV7BH2EtMNPYNJmTzNyz0TNQ%3D&reserved=0 ) we are using a longitudinal metric of respiratory function such that lower scores => worse respiratory function. I reran one of the models from the paper in original JMBayes and in JMBayes2:

lme_ref <- lme(SNIP.occ ~ ns(onset2snip_mnths, 2) , random = ~ ns(onset2snip_mnths, 2) | ID, data = df.snip_72m)

coxC9 <- coxph(Surv(time = dx_delay, time2 = surv_t2_months, fail) ~ age_on + dx_delay2 + simp_site + C9, data = df.clin_72m, x = TRUE, model = TRUE)

JMB1_mod <- jointModelBayes(lme_ref, coxC9, timeVar = "onset2snip_mnths", verbose = FALSE) library(JMbayes2) JMB2_mod <- jm(coxC9, lme_ref, time_var = "onset2snip_mnths")

Here are the Event process results from both: [Screenshot 2021-01-14 at 17 45 07]https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F41968663%2F104621417-61249980-5690-11eb-81ec-2f7054506bb8.png&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cc4afac1c59984532994308d8b8ac712d%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462398058117155%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=jCOyqdwN7EiwZZ6489DQWqQN1qCltUpo4ebyy8A4yDw%3D&reserved=0

As you can see the AssocT term has flipped sign in the JMB2 model. So this would then imply that decreased respiratory function is associated with improved survival -> which doesn't make sense given background knowledge.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdrizopoulos%2FJMbayes2%2Fissues%2F1%23issuecomment-760321037&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cc4afac1c59984532994308d8b8ac712d%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462398058117155%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=xFvXtrwhVTgfGYdu73rKMN2cEjF7e75k4p16BztGdGw%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FADE7TT4BNFHMRRYQZ3YMPT3SZ4ODXANCNFSM4WCL2N4Q&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cc4afac1c59984532994308d8b8ac712d%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462398058117155%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=adaEXQOVsiXwg3D28Vl69Ow%2F4xQQ%2F4g9tEBJaHU9cPU%3D&reserved=0.

jpkrooney commented 3 years ago

Yes no problem! I just fit this:

JMB2_mod2 <- jm(coxC9, lme_ref, time_var = "onset2snip_mnths",
                priors = list(A_tau_bs_gammas = 1, B_tau_bs_gammas = 0.01))

And the results are not much different than the default priors:

Screenshot 2021-01-14 at 19 12 40
drizopoulos commented 3 years ago

I just pushed a new version (JMbayes2_0.1-3) on GitHub - could you try with this one.

From: jpkrooney notifications@github.com Sent: Thursday, January 14, 2021 7:13 PM To: drizopoulos/JMbayes2 JMbayes2@noreply.github.com Cc: D. Rizopoulos d.rizopoulos@erasmusmc.nl; Comment comment@noreply.github.com Subject: Re: [drizopoulos/JMbayes2] JMBayes2 results don't seem to match JMBayes results (#1)

Yes no problem! I just fit this:

JMB2_mod2 <- jm(coxC9, lme_ref, time_var = "onset2snip_mnths",

            priors = list(A_tau_bs_gammas = 1, B_tau_bs_gammas = 0.01))

And the results are not much different than the default priors:

[Screenshot 2021-01-14 at 19 12 40]https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F41968663%2F104631403-7dc6ce80-569c-11eb-92a6-e403097d46f8.png&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cceec210393284c9b169508d8b8b81487%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462448048527371%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=E4qe%2B7K0DxX9%2FEVOKs0KxXvV7FxTC9NaZy0ybiRNBaU%3D&reserved=0

- You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdrizopoulos%2FJMbayes2%2Fissues%2F1%23issuecomment-760371530&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cceec210393284c9b169508d8b8b81487%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462448048527371%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=uq%2F5xJ7HMoxS0N0PfFpmQXELKZRGCBHWm5sj5pp7jBU%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FADE7TT2MTP63BTKXMKWGIYDSZ4X4DANCNFSM4WCL2N4Q&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cceec210393284c9b169508d8b8b81487%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462448048537327%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=oCxaAs07ta1cZHfZmYRfSaaxdpaKPT5nJQUVfFU5CxU%3D&reserved=0.

jpkrooney commented 3 years ago

Hi Dimitris - I just ran the models again after updating. They are in much closer agreement now between JMBayes original, JMBayes2 and JMBayes2 with the tweaked priors:

Screenshot 2021-01-15 at 00 46 56

If I understand your code tweaks properly, this was something to do with the position of spline knots for the baseline hazards ?

drizopoulos commented 3 years ago

Yes, indeed. I had changed the default position of the knots to improve convergence, but apparently this was not a great idea. Thanks for the feedback.

—— Dimitris Rizopoulos Professor of Biostatistics Erasmus University Medical Center The Netherlands

From: jpkrooney notifications@github.com Sent: Friday, January 15, 2021 12:52:30 AM To: drizopoulos/JMbayes2 JMbayes2@noreply.github.com Cc: D. Rizopoulos d.rizopoulos@erasmusmc.nl; Comment comment@noreply.github.com Subject: Re: [drizopoulos/JMbayes2] JMBayes2 results don't seem to match JMBayes results (#1)

Hi Dimitris - I just ran the models again after updating. They are in much closer agreement now between JMBayes original, JMBayes2 and JMBayes2 with the tweaked priors:

[Screenshot 2021-01-15 at 00 46 56]https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fuser-images.githubusercontent.com%2F41968663%2F104662701-afa35980-56cb-11eb-89b7-360cd3aed4b0.png&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C6237a7705db44e9954af08d8b8e77573%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462651535010910%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=8wrukbjmxYiaCylnAjz3p4bTqO2Noivffx0cTLqbBq4%3D&reserved=0

If I understand your code tweaks properly, this was something to do with the position of spline knots for the baseline hazards ?

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdrizopoulos%2FJMbayes2%2Fissues%2F1%23issuecomment-760551000&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C6237a7705db44e9954af08d8b8e77573%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462651535020866%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=8kd%2Fb1O%2FWItuRK82PATr5VhNZyFKxbgFrEj0QtMjgO8%3D&reserved=0, or unsubscribehttps://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FADE7TTZ6HYUWPBX35FJS4XDSZ57T5ANCNFSM4WCL2N4Q&data=04%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C6237a7705db44e9954af08d8b8e77573%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637462651535020866%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=rUz8AzUCE4AFK%2BnTogkgxj7LW0LDZwN2FnCmFFLAnEY%3D&reserved=0.

jpkrooney commented 3 years ago

Thats good to know - thanks for looking into this so quickly!