Closed prockenschaub closed 3 months ago
Thanks for the detailed example. I will study it and see what changes will need to be made to the underlying code of predict()
.
I changed the underlying code in the current development version on GitHub. If you have some time, could you check your examples to see if the results are better now?
I reran my examples and to output looks, agreement between the linear predictor and the predictions is high again.
Problem
I am looking to use the models obtained via JMbayes2 for risk prediction. However, like others before me (#72, #81, #90), I noticed that the AUC of the joint model was considerably worse than that of a simple Cox regression using baseline variables. The explanations that had been given so far were that a) different definitions of AUC/c-index might have been used (#72, #90) and b) the longitudinal covariate may not be a good predictor (#81).
Additional experiments
In my case, the performance of the joint model was so much worse that I decided to dig a little deeper. In particular, I decided to look into the correlation of linear predictor at
t=0
and the cumulative incidence provided bypredict
. Even if the longitudinal predictor is a bad predictor -- i.e., it does not add additional information -- we would expect these to be at least moderately correlated. I adapted the prediction tutorial as follows:Counter to my expectations, there was practically no correlation ($\rho \approx 0.03)$ between linear predictor and cumulative incidence. This seemed odd.
Potential bug
I started debugging the above code to see what was going on under the hood. I stumbled across the
SurvData_HazardModel()
, which to the best of my (limited) understanding expands the prediction frame with the GK points that are used to calculate the hazard in between prediction times. The GK points are given in thetimes_to_fill
.There might be a bug when trying to do so. Consider predicting with the joint model from before, but for readability we will set
GK_k = 7
and only predict for the first two patients:The first call to
SurvData_HazardModel()
occurs inprepare_Data_preds()
, and there was no problem here in my example. We can just continue to debug by typingc
. It the second call toSurvData_HazardModel()
occurs inprepare_DataE_preds()
, something seems to go awry.length(times_to_fill)
is a multiple ofnrow(data)
and the dataset needs to be broadcast so that each patient is matched with their corresponding CK points. However, in our example, we end up with the following:Patient 1 ends up with all zeros, whereas Patient 2 gets the same 7 points between times 0-5 twice, whereas each patient should get one set of each. A
rep(data, times=2)
situation seems to be happening instead of the correctrep(data, each=2)
. This gives a clear reason why the correlation between linear predictor and cumulative incidence was so low above.Potential solution
I haven't fully understood all the parts of JMbayes2 in which
SurvData_HazardModel()
is used. I also don't know if the only time when the bug arises is during prediction. I therefore was only able to try the following quick and dirty fix ofSurvData_HazardModel
and its corresponding call in
prepare_DataE_preds()
This only fixes the replication issue during prediction, leaving any other call in JMbayes2 untouched. By repeating my above example with this code, I get a very high correlation ($\rho \approx 1)$.
Open questions
I haven't been able to figure out yet why the example in the tutorial is able to achive such a high AUC despite this error. Nevertheless, when applying my fix, I at least see a moderate increase from the $AUC = 0.8088$ reported in the tutorial to $AUC \approx 0.83$.