Figure 3 represents the changes in CCI during 10 years following AS, RP, and RT for men with prostate cancer aged 55 and 65, with initial CCI = 0. It also shows absolute differences (with 95 % confidence intervals) in the proportion of men remaining on CCI = 0 by treatment regime (RT/AS and RP). Follow-up time started at treatment initiation. There was a rapid increase in CCI for men treated with RP or RT. This is likely explained by the fact that in Sweden all medical conditions of a patient are included in the discharge diagnoses following an in-hospital episode. Therefore, conditions that were previously not recorded will occur for the first time after hospitalisation despite the fact the patient may have had this condition for a long time. Furthermore, there was a more rapid increase in CCI for men treated with RT than AS. For instance, at 10 years after treatment a higher proportion of men age 65 had died or changed CCI after RT than after RP, 7.3 % (95 % CI: 5.1–9.5). The corresponding value for AS versus RP was 4.1 % (95 % CI: 2.4–5.9).

The following section further clarifies the model specifications and simplifications used in our approach to model temporal changes in comorbidity. Our choice of four-week time steps was a trade-off between precision of dates when CCI changes occurred and computer resources in terms of memory and CPU. A smaller time step would make the long format dataset very large. The implications of our choice of time steps can also be evaluated by looking at time intervals with several CCI changes. For instance, our time steps of four weeks resulted in 2 % of the CCI changes (updated on a daily basis) being pooled with other changes in the same time interval (data not shown). A smaller time step would have alleviated this problem, but at the expense of increased computer load. Moreover, larger time steps would also require a thorough examination of the assumption that changes in CCI follow a Poisson distribution.

It is also important to note that the Poisson distribution, as observed for the size of CCI changes in our example, was purely based on investigation of our data and not justified by any probability theory arguments. If all the CCI changes would have been of size one, a theoretical motivation for the Poisson distribution would have been justifiable. Nevertheless, the weights used to calculate the CCI (1, 2, 3, or 6) make the theoretical basis for usage of the Poisson-distribution limited. Thus, this simplification may need further investigation in other research situations as it is not a given that a Poisson distribution will fit any type of data and, as suggested above, this may also depend on the size of the time steps.

Further justification of the validity of our transition stage models follows from a comparison between simulated and observed data on CCI changes (Fig. 4). Simulated data were obtained using the approach described above, whereas observed data were taken from the traditional Kaplan-Meier estimates with the endpoint of an increase in CCI of at least one unit. Figure 4 shows a comparison for men aged 65 with CCI = 1 at time of prostate cancer diagnosis. The close agreement between the curves for observed and simulated data supports the choice of the models defined above.

Finally, our transition stage model required less model specifications than the traditional time-to-event approach. It captured the complete change in CCI over time using four regression models (three logistic and one Poisson model). This problem would have been difficult to address using traditional survival analysis such as Cox proportional hazards regression. These models only deal with time to a specific event. In the current setting such an event could, for example, be a transition from CCI = 1 to the CCI = 3. More specifically, each arrow in Fig. 1 could correspond to such an event, so that each arrow in Fig. 1 would correspond to a separate Cox model. In practice the cumulated CCI rarely exceeds 15. If we restrict the problem to CCI ≤ 15, the number of models needed to fit would be 16 + 15 + … + 1 = 136. When limiting the problem to CCI ≤ 2, still six Cox models would be needed, three for CCI and three for death. Even when using a Cox model for each size of CCI-change, several different Cox models would be required. Thus, the number of regression analyses needed when applying our proposed method is comparatively small.

Moreover, it is not well-defined how to perform simulations like the ones proposed, based on a set of Cox models. To simulate single event times from a Cox model one could use the estimated hazard function. However, when calculating confidence intervals, a measure of the uncertainty for the hazard function estimate is required. Such an estimate is possible to retrieve using the continuous-time stochastic processes described in [11]. Although possible, this problem is complicated even when CCI ≤ 2, because the variance-covariance matrix is practically infinite dimensional. In our complex simulation, where the results from several different models need to be combined, the mathematical difficulties would be overwhelming.

The simulation in our framework is much simpler since we repeatedly take a small step forward in time. This approach also makes it possible to estimate the uncertainty by using the unscented transform because the number of regression coefficients is limited and the estimated variance-covariance matrix is readily available for each time step.