- Email: [email protected]
Determining disease prevalence from incidence and survival using simulation techniques
Cancer Epidemiology 38 (2014) 193–199
Contents lists available at ScienceDirect
Cancer Epidemiology The International Journal of Cancer Epidemiology, Detection, and Prevention journal homepage: www.cancerepidemiology.net
Determining disease prevalence from incidence and survival using simulation techniques Simon Crouch *, Alex Smith, Dan Painter, Jinlei Li, Eve Roman Epidemiology & Cancer Statistics Group, Department of Health Sciences, University of York, YO10 5DD, UK
A R T I C L E I N F O
A B S T R A C T
Article history: Received 21 August 2013 Received in revised form 14 January 2014 Accepted 17 February 2014 Available online 18 March 2014
Objectives: We present a new method for determining prevalence estimates together with estimates of their precision, from incidence and survival data using Monte-Carlo simulation techniques. The algorithm also provides for the incidence process to be marked with the values of subject level covariates, facilitating calculation of the distribution of these variables in prevalent cases. Methods: Disease incidence is modelled as a marked stochastic process and simulations are made from this process. For each simulated incident case, the probability of remaining in the prevalent subpopulation is calculated from bootstrapped survival curves. This algorithm is used to determine the distribution of prevalence estimates and of the ancillary data associated with the marks of the incidence process. This is then used to determine prevalence estimates and estimates of the precision of these estimates, together with estimates of the distribution of ancillary variables in the prevalent subpopulation. This technique is illustrated by determining the prevalence of acute myeloid leukaemia from data held in the Haematological Malignancy Research Network (HMRN). In addition, the precision of these estimates is determined and the age distribution of prevalent cases diagnosed within twenty years of the prevalence index date is calculated. Conclusion: Determining prevalence estimates by using Monte-Carlo simulation techniques provides a means of calculation more flexible that traditional techniques. In addition to automatically providing precision estimates for the prevalence estimates, the distribution of any measured subject level variables can be calculated for the prevalent sub-population. Temporal changes in incidence and in survival offer no difficulties for the method. ß 2014 Elsevier Ltd. All rights reserved.
Keywords: Prevalence Incidence Survival Simulation Monte-Carlo Prevalence distribution
1. Introduction Estimation of disease prevalence is of fundamental interest in epidemiology [1]. As observed by Gigli et al. [2], three main methods of estimation of prevalence are commonly employed: cross-sectional population survey; direct count of cases in a disease register; and mathematical modelling based on incidence and survival rates. Gigli et al. [2] illustrate a method combining the latter two approaches in three steps: step 1 counts surviving cases at an index date from incident cases in a registry; step 2 estimates the number of prevalent cases lost-to-follow-up from the registry count; and step 3 estimates the number of prevalent cases at the index date that were incident before the start of the registry. Steps
* Corresponding author at: Epidemiology & Cancer Statistics Group, Department of Health Sciences, Seebohm-Rowntree Building, University of York, YO10 5DD, UK. Tel.: +44 01904 321938. E-mail address: [email protected] (S. Crouch). http://dx.doi.org/10.1016/j.canep.2014.02.005 1877-7821/ß 2014 Elsevier Ltd. All rights reserved.
1 and 2 together are often referred to as the ‘‘counting method’’ and step 3 as the ‘‘completeness index method’’. Previous approaches to steps 2 and 3 of this schema have focussed on various analytic techniques of estimation [3–6], themselves based on the relationships between the various measurable quantities [7–9], or by direct modelling [10]. These techniques have found wide application in the literature [11–13]. Consideration of the precision of prevalence estimates has focussed on the variation implied by considering the incidence process as Poisson [2,14]. In this paper we will consider techniques of estimation for steps 2 and 3 based entirely on simulation. We will illustrate our techniques using data drawn from a population based cohort of patients diagnosed with haematological malignancies; in particular we will provide prevalence estimates for acute myeloid leukaemia (AML). We first define what we mean by ‘‘prevalence’’. Broadly speaking the prevalence of a disease in a population is the number or proportion of the population alive at some index date,
194
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
Fig. 1. Types of prevalence. Complete prevalence includes diagnosed cases from any time before the index date still in the prevalent population. n-year prevalence includes cases diagnosed within the last n years. nth year prevalence includes cases diagnosed during the single year between (n 1) and n years before the index date.
previously diagnosed with the disease and not removed from the prevalent disease sub-population between diagnosis and the index date (by death or complete cure, for example); this is referred to as ‘‘complete prevalence’’. We also define ‘‘n-year prevalence’’ and ‘‘nth-year prevalence’’ to refer to those in the prevalent subpopulation at the index date having received a diagnosis of the disease in the previous n years or during the nth year before the index date respectively. Therefore n-year prevalence is the sum of kth-year prevalence for k between 1 and n. So, for example, 3-year prevalence refers to all those in the prevalent sub-population on the index date diagnosed in the three years before the index date; 3rd-year prevalence refers to all those in the prevalent population on the index date diagnosed during the third year before the index date (Fig. 1). For simplicity of presentation in this paper, we assume that the only removal mechanism from the prevalent subpopulation is death. The relevance of n-year and nth-year prevalence for different values of n depends upon the use made of the estimates and upon the disease under consideration. n-year and nth-year prevalence estimates for small values of n will typically correspond to periods of intense treatment for acute diseases; for nth-year prevalence, larger values of n will correspond to periods of long term monitoring. For chronic diseases, all values of n are typically of interest. In this paper we present a new method of determining prevalence based on the computationally intensive method of simulation. The advantages of this new method over existing methodology are that it naturally allows for the estimation of the precision of prevalence estimates and also allows for the estimation of ancillary information about the prevalent subpopulation (for example, it allows for the estimation of the age distribution of the prevalent sub-population). In addition more complex modelling of incidence and survival functions than is usually allowed for in current techniques provides no additional obstacle to simulation techniques. 2. Materials The determination by simulation of prevalence from incidence and survival estimates derived from a patient cohort is illustrated with data on patients diagnosed with acute myeloid leukaemia (AML) drawn from the UK’s population-based Haematological Malignancy Research Network (HMRN) [15]. Initiated in 2004, and covering a population of 3.6 million, this unique patient cohort was established to provide robust generalizable data to inform clinical practice and research. Comprehensive information about HMRN is available elsewhere [15] but briefly, all patients newly diagnosed with a haematological malignancy residing in the HMRN region (>2000 patients a year) have full-treatment, response and outcome data collected to clinical trial standards. HMRN has Section 251 support under the NHS Act 2006; enabling the Health and Social Care Information Centre (HSCIC) to routinely link to and
release nationwide information on deaths, subsequent cancer registrations, and Hospital Episode Statistics (HES). Loss-to-followup rates are very low in this registry, thanks to this comprehensive data linkage. In fact, for the small number of subjects that are lostto-follow-up (by emigration from the UK, for example), the actual date of loss is known with precision in this registry. The demographic structure of the region is similar to the demographic structure of the UK as a whole, allowing for reliable generalization from this population to the population of the UK. Incidence data on patients, 18 years and older, diagnosed with AML was available from the HMRN registry for seven years from 01/09/2005 to 31/08/2012. The index date for the calculations of prevalence was taken to be 31/08/2011 and years are taken to run from the first of September to the thirty-first of August. Survival outcome data was available up until 26/03/2013 for patients diagnosed between 01/09/2005 and 31/08/2011. Characteristics of these patients are shown in Table 1. 2.1. Methods We can estimate the number of prevalent cases of a disease at a particular index date by combining information on incidence and survival. An incident case at time t before the index date, characterized by a vector of explanatory variables for survival x, will contribute S(x, t) to the expected number of prevalent cases at the index date, where S(x, t) is the survival function conditional on explanatory variables x. Therefore, if there are n cases incident at times {t1, t2, . . ., tn} each with corresponding survival explanatory variables {x1, x2, . . ., xn}, then the expected number of prevalent cases at the index date T0 is given by
P¼
n X Sðxi ; T 0 t i Þ i¼1
In this paper we take prevalent cases to be those that have ever been diagnosed with the disease under consideration. This can easily be generalized so that prevalence refers to subjects that have not been removed from the prevalent sub-population by other means (such as cure) by taking the end-point for the survival function S to be time to removal from prevalent population rather than simply time to death. Complete prevalence takes the sum over all time before the index date; this generalizes to estimation of nyear and nth-year prevalence by restricting the sum to cover incident cases from the corresponding time period. The value of P can be calculated by simulation. If the times and associated survival explanatory variables of incident cases can be appropriately modelled, and if survival conditional on the explanatory variables can be estimated, then simulation from the incidence model, together with the survival function, will provide an estimate of P. What is more, sources of variation can be taken into account, so that calculations of P from repeated random draws from the incidence model will provide an estimate of the Table 1 Characteristics of the AML patient cohort.
Incidence dates Number of subjects Male Female Age range (years) Median age (IQR) Maximum follow-up (days) Median survival (95% CI) Total follow-up (years) IQR, inter quartile range.
Incidence dataset
Survival dataset
01/09/2005–31/08/2012 1079 592 (55%) 487 (45%) 18.7–97.8 71.9 (60.0–79.8) N/A N/A N/A
01/09/2005–31/08/2011 934 517 (55%) 417 (45%) 19.0–97.8 71.7 (60.0–79.3) 2720 132 (109–159) 1390
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
precision of the estimate of P. Additionally, the distribution of survival explanatory variables, and other ancillary data, in the prevalent sub-population can also be estimated by ‘‘marking’’ each simulated incident case with the values of these data so that the distribution of these data corresponds to their distribution in the observed incident cases. The simulation process identifies those simulated incident cases considered to be alive at the index date and the distribution of the variables of interest for such survivors are calculated. For example, suppose we are interested in calculating the prevalence, and the age and sex distributions of prevalent cases, for a non-communicable disease in which survival depends on age but not sex. Suppose, in addition, that incidence data, by year, for this disease is available from a registry for k years before a fixed index date. Incidence times for such a disease might reasonably be modelled as arising from a Poisson process. If time is broken up by year before the index date then we can simulate the number of incident cases in any particular year by making a single draw from a Poisson distribution with appropriate rate parameter. This rate parameter is itself the value of a random draw from a distribution that represents the uncertainty in the estimation of the overall incidence rate per year. For a Poisson process this distribution is approximately normal with mean given by the observed incidence rate R and variance R/k. Suppose that Ri cases are simulated in year i then, as the distribution of incidence times for a Poisson process is the uniform distribution [20], the incident case times can be simulated by making a random draw from a uniform distribution over the year with mean Ri. Now, for each simulated incident case, it can be ‘‘marked’’ with values for age and sex by making a single random draw from the joint age and sex distribution of the observed cases. (Alternatively, the joint age and sex distribution can be modelled from the observed distribution and a draw made from this model.) As the incidence time of each simulated incident case is known, the probability of survival for this incident case can be calculated from a survival model calculated from the observed registry. The binary survival (yes/no) for each simulated incident case can be simulated as a draw from a Bernoulli distribution with probability equal to the probability of survival. Therefore the survivors from the simulated incident cases can themselves be simulated and their age and sex distributions calculated.
195
As survival is estimated from a finite quantity of data, the simulation should take into account the variability in the survival function estimate. This can be done by using a bootstrap sampling strategy [21]. In the simple bootstrap a sample of size equal to the number of cases is drawn with replacement and the survival function is calculated for this sample. This is then repeated many times; the variation in the repeated estimates of the survival function captures the uncertainty in the estimate of the survival function. For the prevalence simulation algorithm each simulated incident case is paired at random with one of the bootstrap survival function estimates in order to calculate its probability of survival at the index date. One particular difficulty in the estimation of prevalence (especially if we are interested in surviving subjects who have ever had the disease) is that long term survival may be poorly estimated. Longer registration periods are clearly most desirable but if these are not available then we suggest an approach to this problem that provides reasonable approximate upper and lower bounds for survival estimates. For the lower bound, survival is modelled from an observed cohort using parametric survival modelling, allowing for extrapolation beyond observed survival. For an upper bound it is assumed that the disease in question has a ‘‘cure time’’ after diagnosis, such that survival characteristics after the cure time revert to the general population survival. It is possible that the extrapolation in the lower bound may be too optimistic, and that a ‘‘healthy survivor’’ effect renders the upper bound too pessimistic, but in the absence of data beyond the registry, this is the best that can be hoped for. So, if we let SO denote estimated survival from the cohort and SP denote population survival then the survival function S for the ‘‘cure’’ model with cure time TC is
Sðx; tÞ ¼
SO ðx; tÞ SO ðx; T C Þ Sp ðx0 ; tÞ
t < TC t > TC
where x0 denotes that the covariates (which may include, for example, age) are evaluated at time TC. The effect of this equation is to ‘‘glue’’ the population survival curve onto the survival curve estimated from the cohort (Fig. 2).
Fig. 2. Modelled survival by age. Blue lines show predicted survival from a Weibull survival model according to the ages shown. Red lines show the reversion to population survival after the cure time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
196 Table 2 a Measured incidenceb, prevalenceb and model fitc. 2005–2006
2006–2007
2007–2008
2008–2009
2009–2010
2010–2011
2011–2012
70 (4.05) 61 (3.31)
80 (4.62) 76 (4.13)
78 (4.51) 66 (3.58)
82 (4.74) 56 (3.04)
95 (5.49) 70 (3.80)
90 (5.20) 77 (4.18)
97 (5.61) 81 (4.40)
Measured prevalence Male 14 (0.81) Female 8 (0.43)
12 (0.69) 16 (0.87)
12 (0.69) 11 (0.60)
23 (1.33) 14 (0.76)
27 (1.56) 17 (0.92)
42 (2.43) 33 (1.79)
N/A N/A
Predicted (no cure) Male Female
8.4(0.49) 8.1(0.44)
11.1(0.64) 11.6(0.63)
12.9(0.75) 11.8(0.64)
17.0(0.98) 12.3(0.67)
26.8(1.55) 20.5(1.11)
45.5(2.63) 39.2(2.13)
N/A N/A
Predicted (5 year cure) Male 8.9(0.51) Female 8.6(0.47)
11.1(0.64) 11.6(0.63)
12.9(0.75) 11.8(0.64)
17.0(0.98) 12.3(0.67)
26.8(1.55) 20.5(1.11)
45.5(2.63) 39.2(2.13)
N/A N/A
Predicted (3 year cure) Male 12.3(0.71) Female 11.6(0.63)
14.3(0.83) 14.6(0.79)
14.2(0.82) 12.9(0.70)
17.0(0.98) 12.3(0.67)
26.8(1.55) 20.5(1.11)
45.5(2.63) 39.2(2.63)
N/A N/A
d
Incidence Male Female
a
Years run from 1st September to 31st August. Figures given are: number observed or predicted (number per 100,000). c Model fit: no cure model chi-squared = 7.6 on 11 d.f. (p = 0.75); 5 year cure model chi-squared = 7.3 on 11 d.f. (p = 0.78); 3 year cure model chi-squared = 7.1 on 11 d.f. (p = 0.79). d Incidence cases amongst males = 84.7 per year (=4.90 per 100,000 per year); females = 69.6 (=3.78 per 100,000 per year). Incidence rates are consistent with a homogeneous Poisson process with these rates (p = 0.34 for males, p = 0.30 for females; p values by simulation). b
These two approaches may be considered to give us upper and lower bounds on survival and therefore prevalence estimates. We note that this simulation approach deals in terms of numbers of cases rather than population proportions. It is straightforward to generalize from numbers estimated from a particular sample to proportions for populations. Results are shown in both forms in Tables 2 and 3. We summarize this description in the following algorithm and then illustrate the procedure for the cohort of AML patients. 2.2. Algorithm Preliminary stage (a). Using the registry data, determine an appropriate model for the incidence process for the disease. Estimate the parameters required to specify this model so that simulations of incident case times can be made from it. For example, a Poisson model for the number of yearly incident cases, together with a uniform distribution for the incident case times, models the incidence process as a Poisson process. Preliminary stage (b). Use survival data from the registry to model disease survival using an appropriate model with
explanatory variables suitable for the disease. Calculate survival functions for simple bootstrap samples and store these for use in the prevalence simulation. Step 1. Divide the time before the index date into k suitable time periods, such as years, counting time backwards from the index date. For each time period, simulate the number and times of incident cases during that time period using the model determined in Preliminary stage (a). For each simulated incident case, simulate the associated survival explanatory variables and other ancillary variables by making a random draw from the observed joint distribution of these variables or from a joint model determined from the observed data. Step 2. For each simulated incident case m, calculate its expected probability of survival pm at the index date, using a random survival function (with appropriate values of explanatory variables) from the bootstrap distribution of survival curves prepared in Preliminary stage (b). Step 3. For each simulated incident case m, make a random draw from a Bernoulli distribution with probability pm as parameter. This gives a binary variable Bm with value 1
Table 3 Predicted prevalence.a 5 year prevalence
10 year prevalence
15 year prevalence
20 year prevalence
No cure model Male numbersb Male proportionsb Female numbers Female proportions
116 (95, 137) 6.71 (5.46, 7.95) 91 (70, 112) 4.94 (3.81, 6.08)
161 (137, 185) 9.30 (7.90, 10.7) 128 (104, 152) 6.96 (5.66, 8.27)
188 (162, 215) 10.9 (9.36, 12.4) 155 (128, 181) 8.41 (6.97, 9.84)
209 (181, 237) 12.1 (10.4, 13.7) 175 (147, 203) 9.52 (8.00, 11.0)
5 year cure model Male numbers Male proportions Female numbers Female proportions
116 (95, 137) 6.71 (5.49, 7.92) 91 (70, 112) 4.94 (3.78, 6.10)
172 (147, 198) 9.95 (8.47, 11.4) 138 (112, 163) 7.47 (6.10, 8.84)
222 (192, 251) 12.8 (11.1, 14.5) 183 (156, 211) 9.96 (8.47, 11.5)
266 (234, 298) 15.4 (13.5, 17.2) 226 (195, 257) 12.3 (10.6, 14.0)
3 year cure model Male numbers Male Proportions Female numbers Female proportions
116 (94, 138) 6.71 (5.45, 7.96) 91 (68, 114) 4.94 (3.72, 6.16)
188 (160, 215) 10.9 (9.26, 12.4) 151 (123, 179) 8.18 (6.66, 9.70)
254 (222, 285) 14.7 (12.9, 16.5) 211 (179, 242) 11.4 (9.74, 13.1)
312 (276, 347) 18.0 (16.0, 20.0) 265 (231, 299) 14.4 (12.5, 16.3)
a b
Index date is 31st August 2011. Figures given are: estimate (95% confidence interval). Proportions are given as per 100,000 of the population.
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
corresponding to the simulated case m being a prevalent case at the index date and value 0 corresponding to a case that has died before the index date. Step 4. Sum the values of Bm over all the simulated incident cases m for each of the k time periods to give the values bi (1 i k). Each of the bi corresponds to the number of simulated prevalent cases at the index date that were incident during the time period i. So if years were chosen as the time periods in step 1, then each bn estimates nth-year prevalence and the sum of the bi (1 i n) estimates nyear prevalence. Step 5. For each incident case m with Bm = 1, collect the values of the survival explanatory variables and other ancillary variables. These can either simply be collected variable by variable over all the k time periods or they can be collected variable by variable for each of the k time periods (if estimates of the distributions of these variables are required for each time period before the index date). So far, this algorithm has calculated a single estimate for each of nth-year and n-year prevalence (if the k time periods were chosen as years) and of the distribution of the explanatory and ancillary variables. The algorithm is now repeated N_sim times in order to obtain the sampling distribution of these estimates. For example, in the case of 20-year prevalence P20 we would have calculated N_sim estimates of P20. The mean of the distribution of these estimates provides an overall estimate of P20 and the 2.5% and 97.5% quantiles of the distribution provide a 95% confidence interval for the estimate of P20. Estimates and confidence intervals for estimates of the other quantities calculated by this algorithm can be determined in similar fashion. It is clear that this algorithm facilitates step 3 of the prevalence calculation schema of Gigli et al. [2]: estimates of incidence and survival are made from the registry data and these are used to generate prevalence estimates for the cases incident before the start of the registry. What is perhaps less clear is that the algorithm also provides estimates for the years covered by the disease registry that are unbiased with respect to loss to follow-up, provided that survival estimates are based on survival times that are censored at the time of loss to follow-up. Lost-to-follow-up subjects are simulated using an analogue of the above algorithm and the resulting estimates of prevalent lost-to-follow-up cases are added to the observed totals.
197
We have calculated prevalence estimates for AML under the assumptions of no cure, 3-year cure and 5-year cure. In the simulation, the age distribution of the simulated incident cases was simulated by drawing, for each year, a sample with size equal to the Poisson rate for that year from the actual age distribution of cases. For each prevalent simulated case at the index date, age was recorded. Finally, means of the number of male and female prevalent cases were taken together with appropriate quantiles to estimate 95% confidence intervals. For this exercise, N_sim was set to 1000. In general N_sim may be chosen to ensure stable repeatable results to a chosen degree of precision. All calculations were performed in R version 3.0.1 [17] using the ‘‘survival’’ [18] and ‘‘rms’’ [19] libraries. 3. Results Incidence (see Table 2) amongst males was estimated at 84.7 cases per year (4.90 per 100,000 per year); amongst females at 69.6 cases per year (3.78 per 100,000 per year). The observed incidence by year was consistent with being drawn from a homogeneous Poisson process with these rates (p = 0.34 for males, p = 0.30 for females; p values by simulation). Of the 934 subjects, 58 died on the same day as or before final diagnosis was made. These subjects were allotted a nominal 1-day survival time. Overall survival, estimated by the Kaplan–Meier method, is shown in the upper figure of Fig. 3. Weibull survival modelling showed age at diagnosis to be a highly significant predictor (p < 1010). Modelling the effect of age on survival using spline techniques showed a slight (but statistically significant) non-linear effect; however, the magnitude of this non-linearity
2.3. AML example For the AML data, incidence was modelled as a homogeneous Poisson process, consistent with the observed incidence in the registry (Table 2, note c). Observed incidence rates were calculated for both sexes separately for each of the seven years for which complete data were available and the mean of these yearly rates was taken as the yearly Poisson rate. For each year before the index date, the number of incident cases was simulated by a draw from a normal distribution with mean equal to the Poisson rate and variance given by the Poisson rate divided by seven. The times of these incident cases in each year were then drawn from a uniform distribution, corresponding to the incidence times of a homogeneous Poisson process. AML cohort survival was modelled using a parametric Weibull model with age and sex as explanatory variables. Coefficients from fitting N_sim bootstrap samples were pre-calculated and stored for the main simulation. Population survival by age and sex was estimated using the LSHTM Tables[16] for the same region as the AML cohort, averaged over the final five years of data (2005–2009).
Fig. 3. (Upper) Overall survival in the sample. (Lower) Modelled survival for 70year-old men (blue lines) and women (red lines) compared with observed survival (green solid line, dashed line 95% confidence interval). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
198
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
Fig. 4. Age distribution of incident cases. Black dots represent observed number of incident cases in the registry by year of age at diagnosis; the blue line is a smoothing fit to this distribution. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
was small and had little effect on survival predictions; therefore age was modelled linearly in the final procedure. Sex, on the other hand was not associated with survival (p = 0.23). However, we retain sex in the model because for long term survival deriving from population data, sex is a significant explanatory variable for mortality. Survival model fit was assessed by inspection of predicted survival for ages 20, 30, 40, 50, 60, 70, 80, 90 ages compared with the Kaplan Meier estimate for all ages three years either side of each reference age. An example (for 70-year olds) is shown in the lower figure of Fig. 3. The distribution of age at diagnosis among the incident cases is shown in Fig. 4 together with a smooth of the plot. Fig. 2 shows predicted survival from the Weibull model, together with the predicted survival from the 5year cure model (in red). The overall simulation model was internally validated (see Table 2) by predicting the prevalent cases arising from the known number of incident cases for the six years 2005–2011 and by comparing this prediction with the actual number of prevalent cases at the index date. Model fit gave no cause for concern: no cure model chi-squared = 7.6 on 11 d.f. (p = 0.75); 5 year cure model chi-squared = 7.3 on 11 d.f. (p = 0.78); 3 year cure model chisquared = 7.1 on 11 d.f. (p = 0.79). Prevalence estimates from the model, together with 95% confidence intervals, are shown in Table 3. Prevalence for the first six years before the index date was simply counted from the observed data. Prevalence from before six years prior to the index date was estimated from the simulation model. Fig. 5 plots the nth year (upper figure) and the n-year (lower figure) prevalence by year before index date. As an illustration of the method, Fig. 6 shows the age distribution of incident cases as well as the age distribution of 20-year prevalent cases (assuming the 5-year cure model) at the index date. It is clear from Fig. 6 that the 20-year prevalent population is younger (median age 56.0) than the incident population (median age 71.7).
Fig. 5. Estimated nth year prevalence (upper) and cumulative n-year prevalence (lower) by year before index date. Male subjects are represented by blue dots and line; females by red dots and lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
estimation of the precision of these estimates. In addition to the conceptual simplicity of the approach, simulation allows for the modelling of more complex incidence and survival characteristics. For example, seasonal variation in incidence might be modelled by a periodic non-homogenous Poisson process. Similarly, temporal changes in the incidence process (due, for example, to changes in disease classification or to new preventative measures) or in the survival process (due, for example, to improvements in treatment)
4. Discussion We have presented an approach to the computation of estimates of disease prevalence based on the three-step schema of Gigli et al. [2], substituting simulation based techniques for analytical techniques in steps 2 and 3. The increasing availability of powerful computing resources has meant that epidemiological calculations are increasingly amenable to Monte-Carlo simulation techniques. Estimation of the prevalence of disease in the population from incidence and survival characteristics is ideally suited to a simulation approach, as is the
Fig. 6. Incidence (red) and prevalence (green) age distributions for 20 year prevalence (5 year cure model). Distribution overlap shows in darker green. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
can easily be incorporated into the simulation, provided that they can be modelled from observed data. Simulation also facilitates the calculation of other quantities that are not so easily calculated by conventional means. We have given the example of the prevalence age distribution compared to the incidence age distribution above; provided that the distribution of subject characteristics can be measured and sampled from, this example can be generalized to these other characteristics. We also note that the method allows for the calculation of n-year and nth year prevalence for any n, or indeed for the calculation of prevalence at an index date arising from incidence cases from any time interval before that date; similarly future projections of prevalence are easily calculated, provided that future incidence and survival can be modelled (using age-period-cohort techniques or discrete time models [22], for example). As we have noted, the simulation approach very naturally allows for the estimation of the precision of prevalence estimates, taking into account all the possible sources of variation: the precision of the incidence rate estimate; the stochastic nature of the incidence process; and the precision of the survival estimate. Previous authors [2,14] have noted that the precision of the prevalence estimate will be approximately governed by the variance of the incidence process; our simulations have confirmed this approximation, but have allowed for more precise estimates of this precision. We have illustrated our new technique with the example of acute myeloid leukaemia (AML). From the point of view of prevalence calculations, AML is a disease for which overall survival is poor (3 year survival 18%, Fig. 3 upper figure) but for which cure is possible with increasing probability with decreasing age (Fig. 2). In addition, we observe that the age distribution of incident cases is heavily skewed towards old age (Fig. 4). The consequence of these facts is that the age distribution of prevalent cases is markedly different from that of incident cases (Fig. 6). Although we have illustrated this effect for 20-year prevalence, the same remains true for n-year prevalence estimates for small n (data not shown); survival is so poor for elderly patients that the form of the long term prevalence age distribution is reached rapidly. In conclusion, the new technique for the calculation of disease prevalence that we have presented in this paper provides a means for prevalence estimation from incidence and survival data that is more flexible than the previous techniques that have appeared in the literature. Although we have illustrated the technique with the example of AML, an acute disease with poor overall survival, the technique not limited to diseases with such characteristics. Provided that incidence and survival can be modelled for the population under study, the simulation technique can be applied. Source of funding The authors are supported by a programme grant from Lymphoma and Leukaemia Research (http://leukaemialymphomaresearch.org.uk/).
199
Conflicts of interest No conflicts of interest declared.
References [1] Greenland S, Rothman KJ. Measures of occurrence. In: Rothman KJ, Greenland S, Lash TL, eds. Modern epidemiology. 3rd ed., Philadelphia: Lippincott Williams & Wilkins, 2012: 32–50. [2] Gigli A, Mariotti A, Clegg L, Tavilla A, Corazziari I, Capocaccia R, et al. Estimating the variance of cancer prevalence from population-based registries. Stat Methods Med Res 2006;15:235–53. [3] Verdecchia A, Capocaccia R, Egidi V, Golini A. A method for the estimation of chronic disease morbidity and trends from mortality and survival data. Stat Med 1989;8:201–6. [4] Capocaccia R, De Angelis R. Estimating the completeness of prevalence based on cancer registry data. Stat Med 1997;16:425–40. [5] Gail M, Kessler L, Midthune D, Scoppa S. Two approaches for estimating disease prevalence from population-based registries of incidence and total mortality. Biometrics 1999;55:1137–44. [6] Verdecchia A, De Angelis G, Capocaccia R. Estimation and projections of cancer prevalence from cancer registry data. Stat Med 2002;21:3511–26. [7] Preston P. Relations among standard epidemiological measures in a population. Am J Epidemiol 1987;126:336–45. [8] Keiding N. Age-specific incidence and prevalence: a statistical perspective. J R Stat Soc A 1991;154:371–412. [9] Capocaccia R. Relationships between incidence and mortality in nonreversible diseases. Stat Med 1993;12:2395–415. [10] Maddams J, Brewster D, Gavin A, Steward J, Elliott J, Utley M, et al. Cancer prevalence in the United Kingdom: estimates for 2008. Br J Cancer 2009;101: 541–7. [11] Merrill R, Capocaccia R, Feuer E, Mariotti A. Cancer prevalence estimates based on tumour registry data in the Surveillance, Epidemiology, and End Result (SEER) program. Int J Epidemiol 2000;29:197–207. [12] Pisani P, Bray F, Parkin D. Estimates of the world-wide prevalence of cancer for 25 sites in the adult population. Int J Cancer 2002;97:72–81. [13] Maddams J, Utley M, Mo¨ller H. Projections of cancer prevalence in the United Kingdom, 2010–2040. Br J Cancer 2012;107:1195–202. [14] Capocaccia R, Colonna M, Corazziari I, De Angelis R, Francisci S, Micheli A, et al. Measuring cancer prevalence in Europe: the EUROPREVAL project. Ann Oncol 2002;13:831–9. [15] Smith A, Roman E, Howell D, Jones R, Patmore R, Jack A. The Haematological Malignancy Research Network (HMRN): a new information strategy for population based epidemiology and health service research. Br J Haematol 2009;148:739–53. [16] Cancer Research UK Cancer Survival Group. strel computer program and life tables for cancer survival analysis. UK: Department of Non-Communicable Disease Epidemiology, London School of Hygiene & Tropical Medicine, 2013, Available from www.lshtm.ac.uk/ncde/cancersurvival/tools (accessed 30.07.13).. [17] R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing, 2013, http://www.R-project.org/. [18] Therneau TM, Grambsch PM. Modeling survival data: extending the Cox model. New York: Springer, 2000. [19] Harrell Jr FE. rms: Regression Modeling Strategies. R package version 3. 6-3; 2013, http://CRAN.R-project.org/package=rms. [20] Rizzo NL. Statistical computing with R. London: Chapman & Hall, 2008. [21] Davison AC, Hinkley DV. Bootstrap methods and their application. Cambridge: Cambridge University Press, 1997. [22] Fiorentino F, Maddams J, Mo¨ller H, Utley M. Modelling to estimate future trends in cancer prevalence. Health Care Manage Sci 2011;14:262–6.
Contents lists available at ScienceDirect
Cancer Epidemiology The International Journal of Cancer Epidemiology, Detection, and Prevention journal homepage: www.cancerepidemiology.net
Determining disease prevalence from incidence and survival using simulation techniques Simon Crouch *, Alex Smith, Dan Painter, Jinlei Li, Eve Roman Epidemiology & Cancer Statistics Group, Department of Health Sciences, University of York, YO10 5DD, UK
A R T I C L E I N F O
A B S T R A C T
Article history: Received 21 August 2013 Received in revised form 14 January 2014 Accepted 17 February 2014 Available online 18 March 2014
Objectives: We present a new method for determining prevalence estimates together with estimates of their precision, from incidence and survival data using Monte-Carlo simulation techniques. The algorithm also provides for the incidence process to be marked with the values of subject level covariates, facilitating calculation of the distribution of these variables in prevalent cases. Methods: Disease incidence is modelled as a marked stochastic process and simulations are made from this process. For each simulated incident case, the probability of remaining in the prevalent subpopulation is calculated from bootstrapped survival curves. This algorithm is used to determine the distribution of prevalence estimates and of the ancillary data associated with the marks of the incidence process. This is then used to determine prevalence estimates and estimates of the precision of these estimates, together with estimates of the distribution of ancillary variables in the prevalent subpopulation. This technique is illustrated by determining the prevalence of acute myeloid leukaemia from data held in the Haematological Malignancy Research Network (HMRN). In addition, the precision of these estimates is determined and the age distribution of prevalent cases diagnosed within twenty years of the prevalence index date is calculated. Conclusion: Determining prevalence estimates by using Monte-Carlo simulation techniques provides a means of calculation more flexible that traditional techniques. In addition to automatically providing precision estimates for the prevalence estimates, the distribution of any measured subject level variables can be calculated for the prevalent sub-population. Temporal changes in incidence and in survival offer no difficulties for the method. ß 2014 Elsevier Ltd. All rights reserved.
Keywords: Prevalence Incidence Survival Simulation Monte-Carlo Prevalence distribution
1. Introduction Estimation of disease prevalence is of fundamental interest in epidemiology [1]. As observed by Gigli et al. [2], three main methods of estimation of prevalence are commonly employed: cross-sectional population survey; direct count of cases in a disease register; and mathematical modelling based on incidence and survival rates. Gigli et al. [2] illustrate a method combining the latter two approaches in three steps: step 1 counts surviving cases at an index date from incident cases in a registry; step 2 estimates the number of prevalent cases lost-to-follow-up from the registry count; and step 3 estimates the number of prevalent cases at the index date that were incident before the start of the registry. Steps
* Corresponding author at: Epidemiology & Cancer Statistics Group, Department of Health Sciences, Seebohm-Rowntree Building, University of York, YO10 5DD, UK. Tel.: +44 01904 321938. E-mail address: [email protected] (S. Crouch). http://dx.doi.org/10.1016/j.canep.2014.02.005 1877-7821/ß 2014 Elsevier Ltd. All rights reserved.
1 and 2 together are often referred to as the ‘‘counting method’’ and step 3 as the ‘‘completeness index method’’. Previous approaches to steps 2 and 3 of this schema have focussed on various analytic techniques of estimation [3–6], themselves based on the relationships between the various measurable quantities [7–9], or by direct modelling [10]. These techniques have found wide application in the literature [11–13]. Consideration of the precision of prevalence estimates has focussed on the variation implied by considering the incidence process as Poisson [2,14]. In this paper we will consider techniques of estimation for steps 2 and 3 based entirely on simulation. We will illustrate our techniques using data drawn from a population based cohort of patients diagnosed with haematological malignancies; in particular we will provide prevalence estimates for acute myeloid leukaemia (AML). We first define what we mean by ‘‘prevalence’’. Broadly speaking the prevalence of a disease in a population is the number or proportion of the population alive at some index date,
194
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
Fig. 1. Types of prevalence. Complete prevalence includes diagnosed cases from any time before the index date still in the prevalent population. n-year prevalence includes cases diagnosed within the last n years. nth year prevalence includes cases diagnosed during the single year between (n 1) and n years before the index date.
previously diagnosed with the disease and not removed from the prevalent disease sub-population between diagnosis and the index date (by death or complete cure, for example); this is referred to as ‘‘complete prevalence’’. We also define ‘‘n-year prevalence’’ and ‘‘nth-year prevalence’’ to refer to those in the prevalent subpopulation at the index date having received a diagnosis of the disease in the previous n years or during the nth year before the index date respectively. Therefore n-year prevalence is the sum of kth-year prevalence for k between 1 and n. So, for example, 3-year prevalence refers to all those in the prevalent sub-population on the index date diagnosed in the three years before the index date; 3rd-year prevalence refers to all those in the prevalent population on the index date diagnosed during the third year before the index date (Fig. 1). For simplicity of presentation in this paper, we assume that the only removal mechanism from the prevalent subpopulation is death. The relevance of n-year and nth-year prevalence for different values of n depends upon the use made of the estimates and upon the disease under consideration. n-year and nth-year prevalence estimates for small values of n will typically correspond to periods of intense treatment for acute diseases; for nth-year prevalence, larger values of n will correspond to periods of long term monitoring. For chronic diseases, all values of n are typically of interest. In this paper we present a new method of determining prevalence based on the computationally intensive method of simulation. The advantages of this new method over existing methodology are that it naturally allows for the estimation of the precision of prevalence estimates and also allows for the estimation of ancillary information about the prevalent subpopulation (for example, it allows for the estimation of the age distribution of the prevalent sub-population). In addition more complex modelling of incidence and survival functions than is usually allowed for in current techniques provides no additional obstacle to simulation techniques. 2. Materials The determination by simulation of prevalence from incidence and survival estimates derived from a patient cohort is illustrated with data on patients diagnosed with acute myeloid leukaemia (AML) drawn from the UK’s population-based Haematological Malignancy Research Network (HMRN) [15]. Initiated in 2004, and covering a population of 3.6 million, this unique patient cohort was established to provide robust generalizable data to inform clinical practice and research. Comprehensive information about HMRN is available elsewhere [15] but briefly, all patients newly diagnosed with a haematological malignancy residing in the HMRN region (>2000 patients a year) have full-treatment, response and outcome data collected to clinical trial standards. HMRN has Section 251 support under the NHS Act 2006; enabling the Health and Social Care Information Centre (HSCIC) to routinely link to and
release nationwide information on deaths, subsequent cancer registrations, and Hospital Episode Statistics (HES). Loss-to-followup rates are very low in this registry, thanks to this comprehensive data linkage. In fact, for the small number of subjects that are lostto-follow-up (by emigration from the UK, for example), the actual date of loss is known with precision in this registry. The demographic structure of the region is similar to the demographic structure of the UK as a whole, allowing for reliable generalization from this population to the population of the UK. Incidence data on patients, 18 years and older, diagnosed with AML was available from the HMRN registry for seven years from 01/09/2005 to 31/08/2012. The index date for the calculations of prevalence was taken to be 31/08/2011 and years are taken to run from the first of September to the thirty-first of August. Survival outcome data was available up until 26/03/2013 for patients diagnosed between 01/09/2005 and 31/08/2011. Characteristics of these patients are shown in Table 1. 2.1. Methods We can estimate the number of prevalent cases of a disease at a particular index date by combining information on incidence and survival. An incident case at time t before the index date, characterized by a vector of explanatory variables for survival x, will contribute S(x, t) to the expected number of prevalent cases at the index date, where S(x, t) is the survival function conditional on explanatory variables x. Therefore, if there are n cases incident at times {t1, t2, . . ., tn} each with corresponding survival explanatory variables {x1, x2, . . ., xn}, then the expected number of prevalent cases at the index date T0 is given by
P¼
n X Sðxi ; T 0 t i Þ i¼1
In this paper we take prevalent cases to be those that have ever been diagnosed with the disease under consideration. This can easily be generalized so that prevalence refers to subjects that have not been removed from the prevalent sub-population by other means (such as cure) by taking the end-point for the survival function S to be time to removal from prevalent population rather than simply time to death. Complete prevalence takes the sum over all time before the index date; this generalizes to estimation of nyear and nth-year prevalence by restricting the sum to cover incident cases from the corresponding time period. The value of P can be calculated by simulation. If the times and associated survival explanatory variables of incident cases can be appropriately modelled, and if survival conditional on the explanatory variables can be estimated, then simulation from the incidence model, together with the survival function, will provide an estimate of P. What is more, sources of variation can be taken into account, so that calculations of P from repeated random draws from the incidence model will provide an estimate of the Table 1 Characteristics of the AML patient cohort.
Incidence dates Number of subjects Male Female Age range (years) Median age (IQR) Maximum follow-up (days) Median survival (95% CI) Total follow-up (years) IQR, inter quartile range.
Incidence dataset
Survival dataset
01/09/2005–31/08/2012 1079 592 (55%) 487 (45%) 18.7–97.8 71.9 (60.0–79.8) N/A N/A N/A
01/09/2005–31/08/2011 934 517 (55%) 417 (45%) 19.0–97.8 71.7 (60.0–79.3) 2720 132 (109–159) 1390
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
precision of the estimate of P. Additionally, the distribution of survival explanatory variables, and other ancillary data, in the prevalent sub-population can also be estimated by ‘‘marking’’ each simulated incident case with the values of these data so that the distribution of these data corresponds to their distribution in the observed incident cases. The simulation process identifies those simulated incident cases considered to be alive at the index date and the distribution of the variables of interest for such survivors are calculated. For example, suppose we are interested in calculating the prevalence, and the age and sex distributions of prevalent cases, for a non-communicable disease in which survival depends on age but not sex. Suppose, in addition, that incidence data, by year, for this disease is available from a registry for k years before a fixed index date. Incidence times for such a disease might reasonably be modelled as arising from a Poisson process. If time is broken up by year before the index date then we can simulate the number of incident cases in any particular year by making a single draw from a Poisson distribution with appropriate rate parameter. This rate parameter is itself the value of a random draw from a distribution that represents the uncertainty in the estimation of the overall incidence rate per year. For a Poisson process this distribution is approximately normal with mean given by the observed incidence rate R and variance R/k. Suppose that Ri cases are simulated in year i then, as the distribution of incidence times for a Poisson process is the uniform distribution [20], the incident case times can be simulated by making a random draw from a uniform distribution over the year with mean Ri. Now, for each simulated incident case, it can be ‘‘marked’’ with values for age and sex by making a single random draw from the joint age and sex distribution of the observed cases. (Alternatively, the joint age and sex distribution can be modelled from the observed distribution and a draw made from this model.) As the incidence time of each simulated incident case is known, the probability of survival for this incident case can be calculated from a survival model calculated from the observed registry. The binary survival (yes/no) for each simulated incident case can be simulated as a draw from a Bernoulli distribution with probability equal to the probability of survival. Therefore the survivors from the simulated incident cases can themselves be simulated and their age and sex distributions calculated.
195
As survival is estimated from a finite quantity of data, the simulation should take into account the variability in the survival function estimate. This can be done by using a bootstrap sampling strategy [21]. In the simple bootstrap a sample of size equal to the number of cases is drawn with replacement and the survival function is calculated for this sample. This is then repeated many times; the variation in the repeated estimates of the survival function captures the uncertainty in the estimate of the survival function. For the prevalence simulation algorithm each simulated incident case is paired at random with one of the bootstrap survival function estimates in order to calculate its probability of survival at the index date. One particular difficulty in the estimation of prevalence (especially if we are interested in surviving subjects who have ever had the disease) is that long term survival may be poorly estimated. Longer registration periods are clearly most desirable but if these are not available then we suggest an approach to this problem that provides reasonable approximate upper and lower bounds for survival estimates. For the lower bound, survival is modelled from an observed cohort using parametric survival modelling, allowing for extrapolation beyond observed survival. For an upper bound it is assumed that the disease in question has a ‘‘cure time’’ after diagnosis, such that survival characteristics after the cure time revert to the general population survival. It is possible that the extrapolation in the lower bound may be too optimistic, and that a ‘‘healthy survivor’’ effect renders the upper bound too pessimistic, but in the absence of data beyond the registry, this is the best that can be hoped for. So, if we let SO denote estimated survival from the cohort and SP denote population survival then the survival function S for the ‘‘cure’’ model with cure time TC is
Sðx; tÞ ¼
SO ðx; tÞ SO ðx; T C Þ Sp ðx0 ; tÞ
t < TC t > TC
where x0 denotes that the covariates (which may include, for example, age) are evaluated at time TC. The effect of this equation is to ‘‘glue’’ the population survival curve onto the survival curve estimated from the cohort (Fig. 2).
Fig. 2. Modelled survival by age. Blue lines show predicted survival from a Weibull survival model according to the ages shown. Red lines show the reversion to population survival after the cure time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
196 Table 2 a Measured incidenceb, prevalenceb and model fitc. 2005–2006
2006–2007
2007–2008
2008–2009
2009–2010
2010–2011
2011–2012
70 (4.05) 61 (3.31)
80 (4.62) 76 (4.13)
78 (4.51) 66 (3.58)
82 (4.74) 56 (3.04)
95 (5.49) 70 (3.80)
90 (5.20) 77 (4.18)
97 (5.61) 81 (4.40)
Measured prevalence Male 14 (0.81) Female 8 (0.43)
12 (0.69) 16 (0.87)
12 (0.69) 11 (0.60)
23 (1.33) 14 (0.76)
27 (1.56) 17 (0.92)
42 (2.43) 33 (1.79)
N/A N/A
Predicted (no cure) Male Female
8.4(0.49) 8.1(0.44)
11.1(0.64) 11.6(0.63)
12.9(0.75) 11.8(0.64)
17.0(0.98) 12.3(0.67)
26.8(1.55) 20.5(1.11)
45.5(2.63) 39.2(2.13)
N/A N/A
Predicted (5 year cure) Male 8.9(0.51) Female 8.6(0.47)
11.1(0.64) 11.6(0.63)
12.9(0.75) 11.8(0.64)
17.0(0.98) 12.3(0.67)
26.8(1.55) 20.5(1.11)
45.5(2.63) 39.2(2.13)
N/A N/A
Predicted (3 year cure) Male 12.3(0.71) Female 11.6(0.63)
14.3(0.83) 14.6(0.79)
14.2(0.82) 12.9(0.70)
17.0(0.98) 12.3(0.67)
26.8(1.55) 20.5(1.11)
45.5(2.63) 39.2(2.63)
N/A N/A
d
Incidence Male Female
a
Years run from 1st September to 31st August. Figures given are: number observed or predicted (number per 100,000). c Model fit: no cure model chi-squared = 7.6 on 11 d.f. (p = 0.75); 5 year cure model chi-squared = 7.3 on 11 d.f. (p = 0.78); 3 year cure model chi-squared = 7.1 on 11 d.f. (p = 0.79). d Incidence cases amongst males = 84.7 per year (=4.90 per 100,000 per year); females = 69.6 (=3.78 per 100,000 per year). Incidence rates are consistent with a homogeneous Poisson process with these rates (p = 0.34 for males, p = 0.30 for females; p values by simulation). b
These two approaches may be considered to give us upper and lower bounds on survival and therefore prevalence estimates. We note that this simulation approach deals in terms of numbers of cases rather than population proportions. It is straightforward to generalize from numbers estimated from a particular sample to proportions for populations. Results are shown in both forms in Tables 2 and 3. We summarize this description in the following algorithm and then illustrate the procedure for the cohort of AML patients. 2.2. Algorithm Preliminary stage (a). Using the registry data, determine an appropriate model for the incidence process for the disease. Estimate the parameters required to specify this model so that simulations of incident case times can be made from it. For example, a Poisson model for the number of yearly incident cases, together with a uniform distribution for the incident case times, models the incidence process as a Poisson process. Preliminary stage (b). Use survival data from the registry to model disease survival using an appropriate model with
explanatory variables suitable for the disease. Calculate survival functions for simple bootstrap samples and store these for use in the prevalence simulation. Step 1. Divide the time before the index date into k suitable time periods, such as years, counting time backwards from the index date. For each time period, simulate the number and times of incident cases during that time period using the model determined in Preliminary stage (a). For each simulated incident case, simulate the associated survival explanatory variables and other ancillary variables by making a random draw from the observed joint distribution of these variables or from a joint model determined from the observed data. Step 2. For each simulated incident case m, calculate its expected probability of survival pm at the index date, using a random survival function (with appropriate values of explanatory variables) from the bootstrap distribution of survival curves prepared in Preliminary stage (b). Step 3. For each simulated incident case m, make a random draw from a Bernoulli distribution with probability pm as parameter. This gives a binary variable Bm with value 1
Table 3 Predicted prevalence.a 5 year prevalence
10 year prevalence
15 year prevalence
20 year prevalence
No cure model Male numbersb Male proportionsb Female numbers Female proportions
116 (95, 137) 6.71 (5.46, 7.95) 91 (70, 112) 4.94 (3.81, 6.08)
161 (137, 185) 9.30 (7.90, 10.7) 128 (104, 152) 6.96 (5.66, 8.27)
188 (162, 215) 10.9 (9.36, 12.4) 155 (128, 181) 8.41 (6.97, 9.84)
209 (181, 237) 12.1 (10.4, 13.7) 175 (147, 203) 9.52 (8.00, 11.0)
5 year cure model Male numbers Male proportions Female numbers Female proportions
116 (95, 137) 6.71 (5.49, 7.92) 91 (70, 112) 4.94 (3.78, 6.10)
172 (147, 198) 9.95 (8.47, 11.4) 138 (112, 163) 7.47 (6.10, 8.84)
222 (192, 251) 12.8 (11.1, 14.5) 183 (156, 211) 9.96 (8.47, 11.5)
266 (234, 298) 15.4 (13.5, 17.2) 226 (195, 257) 12.3 (10.6, 14.0)
3 year cure model Male numbers Male Proportions Female numbers Female proportions
116 (94, 138) 6.71 (5.45, 7.96) 91 (68, 114) 4.94 (3.72, 6.16)
188 (160, 215) 10.9 (9.26, 12.4) 151 (123, 179) 8.18 (6.66, 9.70)
254 (222, 285) 14.7 (12.9, 16.5) 211 (179, 242) 11.4 (9.74, 13.1)
312 (276, 347) 18.0 (16.0, 20.0) 265 (231, 299) 14.4 (12.5, 16.3)
a b
Index date is 31st August 2011. Figures given are: estimate (95% confidence interval). Proportions are given as per 100,000 of the population.
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
corresponding to the simulated case m being a prevalent case at the index date and value 0 corresponding to a case that has died before the index date. Step 4. Sum the values of Bm over all the simulated incident cases m for each of the k time periods to give the values bi (1 i k). Each of the bi corresponds to the number of simulated prevalent cases at the index date that were incident during the time period i. So if years were chosen as the time periods in step 1, then each bn estimates nth-year prevalence and the sum of the bi (1 i n) estimates nyear prevalence. Step 5. For each incident case m with Bm = 1, collect the values of the survival explanatory variables and other ancillary variables. These can either simply be collected variable by variable over all the k time periods or they can be collected variable by variable for each of the k time periods (if estimates of the distributions of these variables are required for each time period before the index date). So far, this algorithm has calculated a single estimate for each of nth-year and n-year prevalence (if the k time periods were chosen as years) and of the distribution of the explanatory and ancillary variables. The algorithm is now repeated N_sim times in order to obtain the sampling distribution of these estimates. For example, in the case of 20-year prevalence P20 we would have calculated N_sim estimates of P20. The mean of the distribution of these estimates provides an overall estimate of P20 and the 2.5% and 97.5% quantiles of the distribution provide a 95% confidence interval for the estimate of P20. Estimates and confidence intervals for estimates of the other quantities calculated by this algorithm can be determined in similar fashion. It is clear that this algorithm facilitates step 3 of the prevalence calculation schema of Gigli et al. [2]: estimates of incidence and survival are made from the registry data and these are used to generate prevalence estimates for the cases incident before the start of the registry. What is perhaps less clear is that the algorithm also provides estimates for the years covered by the disease registry that are unbiased with respect to loss to follow-up, provided that survival estimates are based on survival times that are censored at the time of loss to follow-up. Lost-to-follow-up subjects are simulated using an analogue of the above algorithm and the resulting estimates of prevalent lost-to-follow-up cases are added to the observed totals.
197
We have calculated prevalence estimates for AML under the assumptions of no cure, 3-year cure and 5-year cure. In the simulation, the age distribution of the simulated incident cases was simulated by drawing, for each year, a sample with size equal to the Poisson rate for that year from the actual age distribution of cases. For each prevalent simulated case at the index date, age was recorded. Finally, means of the number of male and female prevalent cases were taken together with appropriate quantiles to estimate 95% confidence intervals. For this exercise, N_sim was set to 1000. In general N_sim may be chosen to ensure stable repeatable results to a chosen degree of precision. All calculations were performed in R version 3.0.1 [17] using the ‘‘survival’’ [18] and ‘‘rms’’ [19] libraries. 3. Results Incidence (see Table 2) amongst males was estimated at 84.7 cases per year (4.90 per 100,000 per year); amongst females at 69.6 cases per year (3.78 per 100,000 per year). The observed incidence by year was consistent with being drawn from a homogeneous Poisson process with these rates (p = 0.34 for males, p = 0.30 for females; p values by simulation). Of the 934 subjects, 58 died on the same day as or before final diagnosis was made. These subjects were allotted a nominal 1-day survival time. Overall survival, estimated by the Kaplan–Meier method, is shown in the upper figure of Fig. 3. Weibull survival modelling showed age at diagnosis to be a highly significant predictor (p < 1010). Modelling the effect of age on survival using spline techniques showed a slight (but statistically significant) non-linear effect; however, the magnitude of this non-linearity
2.3. AML example For the AML data, incidence was modelled as a homogeneous Poisson process, consistent with the observed incidence in the registry (Table 2, note c). Observed incidence rates were calculated for both sexes separately for each of the seven years for which complete data were available and the mean of these yearly rates was taken as the yearly Poisson rate. For each year before the index date, the number of incident cases was simulated by a draw from a normal distribution with mean equal to the Poisson rate and variance given by the Poisson rate divided by seven. The times of these incident cases in each year were then drawn from a uniform distribution, corresponding to the incidence times of a homogeneous Poisson process. AML cohort survival was modelled using a parametric Weibull model with age and sex as explanatory variables. Coefficients from fitting N_sim bootstrap samples were pre-calculated and stored for the main simulation. Population survival by age and sex was estimated using the LSHTM Tables[16] for the same region as the AML cohort, averaged over the final five years of data (2005–2009).
Fig. 3. (Upper) Overall survival in the sample. (Lower) Modelled survival for 70year-old men (blue lines) and women (red lines) compared with observed survival (green solid line, dashed line 95% confidence interval). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
198
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
Fig. 4. Age distribution of incident cases. Black dots represent observed number of incident cases in the registry by year of age at diagnosis; the blue line is a smoothing fit to this distribution. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
was small and had little effect on survival predictions; therefore age was modelled linearly in the final procedure. Sex, on the other hand was not associated with survival (p = 0.23). However, we retain sex in the model because for long term survival deriving from population data, sex is a significant explanatory variable for mortality. Survival model fit was assessed by inspection of predicted survival for ages 20, 30, 40, 50, 60, 70, 80, 90 ages compared with the Kaplan Meier estimate for all ages three years either side of each reference age. An example (for 70-year olds) is shown in the lower figure of Fig. 3. The distribution of age at diagnosis among the incident cases is shown in Fig. 4 together with a smooth of the plot. Fig. 2 shows predicted survival from the Weibull model, together with the predicted survival from the 5year cure model (in red). The overall simulation model was internally validated (see Table 2) by predicting the prevalent cases arising from the known number of incident cases for the six years 2005–2011 and by comparing this prediction with the actual number of prevalent cases at the index date. Model fit gave no cause for concern: no cure model chi-squared = 7.6 on 11 d.f. (p = 0.75); 5 year cure model chi-squared = 7.3 on 11 d.f. (p = 0.78); 3 year cure model chisquared = 7.1 on 11 d.f. (p = 0.79). Prevalence estimates from the model, together with 95% confidence intervals, are shown in Table 3. Prevalence for the first six years before the index date was simply counted from the observed data. Prevalence from before six years prior to the index date was estimated from the simulation model. Fig. 5 plots the nth year (upper figure) and the n-year (lower figure) prevalence by year before index date. As an illustration of the method, Fig. 6 shows the age distribution of incident cases as well as the age distribution of 20-year prevalent cases (assuming the 5-year cure model) at the index date. It is clear from Fig. 6 that the 20-year prevalent population is younger (median age 56.0) than the incident population (median age 71.7).
Fig. 5. Estimated nth year prevalence (upper) and cumulative n-year prevalence (lower) by year before index date. Male subjects are represented by blue dots and line; females by red dots and lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
estimation of the precision of these estimates. In addition to the conceptual simplicity of the approach, simulation allows for the modelling of more complex incidence and survival characteristics. For example, seasonal variation in incidence might be modelled by a periodic non-homogenous Poisson process. Similarly, temporal changes in the incidence process (due, for example, to changes in disease classification or to new preventative measures) or in the survival process (due, for example, to improvements in treatment)
4. Discussion We have presented an approach to the computation of estimates of disease prevalence based on the three-step schema of Gigli et al. [2], substituting simulation based techniques for analytical techniques in steps 2 and 3. The increasing availability of powerful computing resources has meant that epidemiological calculations are increasingly amenable to Monte-Carlo simulation techniques. Estimation of the prevalence of disease in the population from incidence and survival characteristics is ideally suited to a simulation approach, as is the
Fig. 6. Incidence (red) and prevalence (green) age distributions for 20 year prevalence (5 year cure model). Distribution overlap shows in darker green. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
S. Crouch et al. / Cancer Epidemiology 38 (2014) 193–199
can easily be incorporated into the simulation, provided that they can be modelled from observed data. Simulation also facilitates the calculation of other quantities that are not so easily calculated by conventional means. We have given the example of the prevalence age distribution compared to the incidence age distribution above; provided that the distribution of subject characteristics can be measured and sampled from, this example can be generalized to these other characteristics. We also note that the method allows for the calculation of n-year and nth year prevalence for any n, or indeed for the calculation of prevalence at an index date arising from incidence cases from any time interval before that date; similarly future projections of prevalence are easily calculated, provided that future incidence and survival can be modelled (using age-period-cohort techniques or discrete time models [22], for example). As we have noted, the simulation approach very naturally allows for the estimation of the precision of prevalence estimates, taking into account all the possible sources of variation: the precision of the incidence rate estimate; the stochastic nature of the incidence process; and the precision of the survival estimate. Previous authors [2,14] have noted that the precision of the prevalence estimate will be approximately governed by the variance of the incidence process; our simulations have confirmed this approximation, but have allowed for more precise estimates of this precision. We have illustrated our new technique with the example of acute myeloid leukaemia (AML). From the point of view of prevalence calculations, AML is a disease for which overall survival is poor (3 year survival 18%, Fig. 3 upper figure) but for which cure is possible with increasing probability with decreasing age (Fig. 2). In addition, we observe that the age distribution of incident cases is heavily skewed towards old age (Fig. 4). The consequence of these facts is that the age distribution of prevalent cases is markedly different from that of incident cases (Fig. 6). Although we have illustrated this effect for 20-year prevalence, the same remains true for n-year prevalence estimates for small n (data not shown); survival is so poor for elderly patients that the form of the long term prevalence age distribution is reached rapidly. In conclusion, the new technique for the calculation of disease prevalence that we have presented in this paper provides a means for prevalence estimation from incidence and survival data that is more flexible than the previous techniques that have appeared in the literature. Although we have illustrated the technique with the example of AML, an acute disease with poor overall survival, the technique not limited to diseases with such characteristics. Provided that incidence and survival can be modelled for the population under study, the simulation technique can be applied. Source of funding The authors are supported by a programme grant from Lymphoma and Leukaemia Research (http://leukaemialymphomaresearch.org.uk/).
199
Conflicts of interest No conflicts of interest declared.
References [1] Greenland S, Rothman KJ. Measures of occurrence. In: Rothman KJ, Greenland S, Lash TL, eds. Modern epidemiology. 3rd ed., Philadelphia: Lippincott Williams & Wilkins, 2012: 32–50. [2] Gigli A, Mariotti A, Clegg L, Tavilla A, Corazziari I, Capocaccia R, et al. Estimating the variance of cancer prevalence from population-based registries. Stat Methods Med Res 2006;15:235–53. [3] Verdecchia A, Capocaccia R, Egidi V, Golini A. A method for the estimation of chronic disease morbidity and trends from mortality and survival data. Stat Med 1989;8:201–6. [4] Capocaccia R, De Angelis R. Estimating the completeness of prevalence based on cancer registry data. Stat Med 1997;16:425–40. [5] Gail M, Kessler L, Midthune D, Scoppa S. Two approaches for estimating disease prevalence from population-based registries of incidence and total mortality. Biometrics 1999;55:1137–44. [6] Verdecchia A, De Angelis G, Capocaccia R. Estimation and projections of cancer prevalence from cancer registry data. Stat Med 2002;21:3511–26. [7] Preston P. Relations among standard epidemiological measures in a population. Am J Epidemiol 1987;126:336–45. [8] Keiding N. Age-specific incidence and prevalence: a statistical perspective. J R Stat Soc A 1991;154:371–412. [9] Capocaccia R. Relationships between incidence and mortality in nonreversible diseases. Stat Med 1993;12:2395–415. [10] Maddams J, Brewster D, Gavin A, Steward J, Elliott J, Utley M, et al. Cancer prevalence in the United Kingdom: estimates for 2008. Br J Cancer 2009;101: 541–7. [11] Merrill R, Capocaccia R, Feuer E, Mariotti A. Cancer prevalence estimates based on tumour registry data in the Surveillance, Epidemiology, and End Result (SEER) program. Int J Epidemiol 2000;29:197–207. [12] Pisani P, Bray F, Parkin D. Estimates of the world-wide prevalence of cancer for 25 sites in the adult population. Int J Cancer 2002;97:72–81. [13] Maddams J, Utley M, Mo¨ller H. Projections of cancer prevalence in the United Kingdom, 2010–2040. Br J Cancer 2012;107:1195–202. [14] Capocaccia R, Colonna M, Corazziari I, De Angelis R, Francisci S, Micheli A, et al. Measuring cancer prevalence in Europe: the EUROPREVAL project. Ann Oncol 2002;13:831–9. [15] Smith A, Roman E, Howell D, Jones R, Patmore R, Jack A. The Haematological Malignancy Research Network (HMRN): a new information strategy for population based epidemiology and health service research. Br J Haematol 2009;148:739–53. [16] Cancer Research UK Cancer Survival Group. strel computer program and life tables for cancer survival analysis. UK: Department of Non-Communicable Disease Epidemiology, London School of Hygiene & Tropical Medicine, 2013, Available from www.lshtm.ac.uk/ncde/cancersurvival/tools (accessed 30.07.13).. [17] R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing, 2013, http://www.R-project.org/. [18] Therneau TM, Grambsch PM. Modeling survival data: extending the Cox model. New York: Springer, 2000. [19] Harrell Jr FE. rms: Regression Modeling Strategies. R package version 3. 6-3; 2013, http://CRAN.R-project.org/package=rms. [20] Rizzo NL. Statistical computing with R. London: Chapman & Hall, 2008. [21] Davison AC, Hinkley DV. Bootstrap methods and their application. Cambridge: Cambridge University Press, 1997. [22] Fiorentino F, Maddams J, Mo¨ller H, Utley M. Modelling to estimate future trends in cancer prevalence. Health Care Manage Sci 2011;14:262–6.